Problems: High-dimensional Regression

Problems: High-dimensional Regression #

Get this problem set

Download everything (Topic2.5-High_dimensional_regression_Problems.zip) — the notebook and solvent_ir_calibration.csv, in a folder that is ready to run as-is.

The Download badge at the top of the page will also give you the notebook on its own. On Vocareum everything is already set up for you.

Before you start

This problem set accompanies High-dimensional Regression. It is worth 100 points:

Part A — Skill Checks (30 pts) — short answers, auto-graded, resubmit as often as you like until they pass.
Part B — Visualization (35 pts) — plots plus written interpretation, peer graded.
Part C — Open Ended (35 pts) — one synthesis problem, peer graded.

The parts build on each other: Part A works out the syntax you need for Part B, and Part B produces the evidence you argue from in Part C. Do them in order.

Setup #

Same solvent-recovery column as Topic 2.4, same FTIR on the overhead vapor, same four analytes: acetone, toluene, n-hexane and 2-propanol. A channel is again one 4 cm⁻¹ wavenumber bin, and there are 880 of them from 450 to 3966 cm⁻¹. Topic 2.4 asked what the design matrix looks like. This set asks the question the instrument was installed for: predict the acetone mole fraction from the spectrum.

The file here adds one thing. It holds two batches:

Batch A, 120 scans, the calibration run, and byte for byte the same rows as Topic 2.4’s file.
Batch B, 60 scans, the same column three months later, after the sample cell was replaced with a longer one and the nitrogen purge started letting more water through.

Everything in Parts A and B uses batch A only. Batch B stays sealed until Part C, which is the point of it: a calibration is only worth anything if it still works on spectra recorded after something about the instrument changed.

What you need to know about the spectra

No infrared spectroscopy is assumed. Two facts about where things absorb, both read off the reference spectra rather than from a textbook:

2850–3000 cm⁻¹ is the C–H stretching region, and all four analytes absorb there, since all four are organic molecules with C–H bonds. It is the largest absorption in the data, and it says how much organic vapor is present rather than which one.
1700–1780 cm⁻¹ is the C=O (carbonyl) stretching region, and of the four analytes only acetone absorbs there. Water also absorbs in that region, which is why the purge matters.

And the two things that differ about batch B, in terms of what they do to a spectrum:

A longer sample cell means the light travels further through the same gas. Absorbance is proportional to path length, so a longer cell multiplies the whole spectrum by a common factor — every band grows together, and the shape is unchanged.
A leakier purge means more water vapor and CO₂ sitting in the beam. Those two have their own bands (water near 1500–1600 and 3700–3900 cm⁻¹, CO₂ near 2310–2360 cm⁻¹), and how much of them is present has nothing to do with the solvent composition, so they add structure that carries no information about acetone.

Part C task 1 asks you to work out which of the two is responsible for most of the difference between the batches.

%matplotlib inline
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.decomposition import PCA
from sklearn.cross_decomposition import PLSRegression
from sklearn.pipeline import make_pipeline
from sklearn.model_selection import train_test_split
try:
    plt.style.use('../settings/plot_style.mplstyle')   # available inside the book
except OSError:
    pass                                              # downloaded notebook: use defaults

df = pd.read_csv('data/solvent_ir_calibration.csv')

composition_cols = [c for c in df.columns if c.startswith('x_')]
channel_cols = [c for c in df.columns
                if c not in composition_cols + ['sample_id', 'batch']]
wavenumber = np.array([float(c) for c in channel_cols])   # cm^-1, 880 of them

is_A = df['batch'] == 'A'
X_A = df.loc[is_A, channel_cols].values
y_A = df.loc[is_A, 'x_acetone'].values
X_B = df.loc[~is_A, channel_cols].values          # Part C only
y_B = df.loc[~is_A, 'x_acetone'].values

# Parts A and B work inside batch A. One split, used by every question.
X_train, X_test, y_train, y_test = train_test_split(
    X_A, y_A, test_size=0.3, random_state=0)

print(f"batch A: {X_A.shape[0]} scans x {X_A.shape[1]} channels")
print(f"batch B: {X_B.shape[0]} scans (untouched until Part C)")
print(f"train / test inside batch A: {X_train.shape[0]} / {X_test.shape[0]}")

batch A: 120 scans x 880 channels
batch B: 60 scans (untouched until Part C)
train / test inside batch A: 84 / 36

Note the shape of the training set: 84 scans, 880 channels. There are ten times more unknowns than equations, and ordinary least squares will still return an answer.

Center, do not standardize

Every model in this set is fit to mean-centered channels, never to standardized ones. PCA centers the columns for you, so make_pipeline(PCA(n_components=k), LinearRegression()) is already doing the right thing and needs no StandardScaler.

PLSRegression needs one change: its constructor defaults to scale=True, which standardizes every column before fitting. Write PLSRegression(n_components=k, scale=False) everywhere in this set.

This is the direct consequence of Topic 2.4 Part B: all 880 channels are the same physical quantity in the same units, so the differences in their variances are real information about where the chemistry is. Standardizing promotes 504 near-empty channels to the same weight as the strongest absorption band. Part B task 4 measures what that costs.

Run this once. It defines the check helper used after each Part A answer. It will not run inside the book — use the downloaded notebook.

# Run this once. It defines `check`, used after each Part A answer below.
import hashlib

_EXPECTED = {
    'q1': (6, ['270108127578', '0a560ed2dbfe', 'e9f92987e2fa', 'ad12040e6410', '39160b503ccd', 'b74bab9eb184', '7730e9c7ca96', '9423f6cb0076', '441a47639ff4', '2174deb09aa0', 'ce5b193fdfcb', '02988fc45243', 'b005eb751a57', '5932b0414b46', '35922f90a104', '9dc8bdea34d8', '2de2230357fe', '140d2ea6df10', '8d84192aaffd', '076af9d5c838', '0707986c0a68']),
    'q2': (5, ['3226cd6eba1d', '07091e00cfeb', '759b30632dab', 'f2e7f5a3da3c', 'fba5f6b169d9', '1100d555290d', '6906688bb0ae', '967480214611', 'c956e3e3d526', '0145995b1889', '4bd1bece3c82', '57c8d0e75d82', '0fd24460bac8']),
    'q3': (5, ['48511fb12846', '7ee13570280e', '4296a012dee1', 'ae286948b388', 'c085824963ee', '77b64d081a6f', '06054b143ffd', '6a6614a247e6', '7c152dfcbb61', '629b226f752b', '50ee5f7f430d', '2b6738331c48', '3c814c21bb42', 'c941a02957e3', 'c288bd9645f8', '1a19f2a980a9', '4168112c04e1', '00aedf423627', '5d089fbbd047', 'dba7b4e39640', '292d7743e0ee']),
}

def check(qid, value):
    """Compare an answer against a stored digest. Practice only -- nothing is recorded."""
    places, digests = _EXPECTED[qid]
    if value is Ellipsis:
        print(f"{qid}: not answered yet")
        return
    got = hashlib.sha256(f"{float(value):.{places}f}".encode()).hexdigest()[:12]
    print(f"{qid}: {'correct' if got in digests else 'not the expected value'}")

Part A — Skill Checks (30 pts)#

Three questions, 10 points each. All three fit on X_train/y_train and score on X_test/y_test, so the three numbers are directly comparable. A1 is the model with no dimensionality reduction at all, and A2 and A3 are the two ways of reducing it to the same two components.

A1. Least squares with more unknowns than equations (10 pts)#

Exercise 1

Fit LinearRegression() on all 880 channels of X_train and score it on X_test.

There are 880 coefficients and 84 training scans, so infinitely many coefficient vectors fit the training data exactly. LinearRegression does not fail. It returns the one with the smallest norm. Print the training \(r^2\) as well as the test \(r^2\); you will need to have seen both for Part B.

Assign the test \(r^2\) to r2_ols_test.

# YOUR CODE HERE
r2_ols_test = ...

check("q1", r2_ols_test)

A2. Two principal components (10 pts)#

Exercise 2

Build a principal component regression with two components:

make_pipeline(PCA(n_components=2), LinearRegression())

Fit it on X_train/y_train and score on X_test/y_test. PCA centers the columns itself, so do not add a StandardScaler.

Assign the test \(r^2\) to r2_pcr2_test.

# YOUR CODE HERE
r2_pcr2_test = ...

check("q2", r2_pcr2_test)

A3. Two PLS components (10 pts)#

Exercise 3

Now the supervised reduction. Fit

PLSRegression(n_components=2, scale=False)

on the same X_train/y_train and score on the same X_test/y_test.

scale=False is not the default and it matters here. See the warning in the setup.

Assign the test \(r^2\) to r2_pls2_test.

# YOUR CODE HERE
r2_pls2_test = ...

check("q3", r2_pls2_test)

Part B — Visualization (35 pts)#

Part A fit three models to one split: every channel, two unsupervised components, and two supervised ones. This part is about the gap between the last two.

Exercise 4

1. The component sweep. For k in range(1, 21), fit four models on X_train and score each on X_test:

make_pipeline(PCA(n_components=k), LinearRegression())                    # PCR, centered
make_pipeline(StandardScaler(), PCA(n_components=k), LinearRegression())  # PCR, standardized
PLSRegression(n_components=k, scale=False)                                # PLS, centered
PLSRegression(n_components=k)                                             # PLS, standardized

Plot all four test \(r^2\) curves on one axes, with the centered pair solid and the standardized pair dashed. Clip the \(y\) axis so the interesting region is readable; the first PCR point is below zero.

Report the smallest k at which each of the four first exceeds 0.99.

In two sentences: say which method is ahead at small k when both are centered, and say whether standardizing helped or hurt.

2. Why PCA orders its components the way it does. Fit PCA(n_components=6) on X_train. For each of the six components, compute its explained variance ratio and the absolute correlation between that component’s score and y_train.

Plot both against component number on one axes — two series, 1 to 6 on the \(x\) axis. Use a second \(y\) axis or normalize them so both are readable together, and mark which component leads on each.

In two sentences: name the component carrying the most variance and the one most strongly related to acetone, and say why an unsupervised method has no way to prefer the second.

3. What the directions look like. On one figure with two panels, plot against wavenumber:

the first PCA loading vector, pca.components_[0]
the first PLS weight vector, pls.x_weights_[:, 0], from a fit with n_components=2

Mark 1738 cm⁻¹ and 2978 cm⁻¹ on both, and report which wavenumber carries the largest absolute value in each vector.

In two sentences: say what each method’s first direction is responding to. Use the two regions from the setup note: 2850–3000 cm⁻¹ is C–H stretching, which all four analytes have, and 1700–1780 cm⁻¹ is the carbonyl, which among the analytes only acetone has.

# YOUR CODE HERE
k_range = range(1, 21)

Part C — Open Ended (35 pts)#

Batch B has been sitting untouched. It is 60 scans of the same column taken three months later, through a longer sample cell and a leakier purge. Nothing about the chemistry changed; the instrument did.

This is the question that decides whether the calibration was worth building. Every number in Parts A and B was measured on scans recorded in the same session as the training data, and none of them tell you what happens next.

Exercise 5

1. Measure what changed. Compare the two batches directly: plot the mean spectrum of batch A and of batch B on one axes, and plot their difference on a second. Report the mean absorbance over all channels for each batch, and the wavenumber at which the mean spectra differ most.

In one sentence: say which of the two stated instrument changes the difference spectrum is mostly showing.

2. Transfer the models you already have. Refit, on all 120 scans of batch A, the three models from Part A — OLS on every channel, PCR at k = 2, and PLS at k = 2 — and score each on batch B. Put the batch-A test scores from Part A beside them in one table.

In one sentence: say which model loses the least and which loses the most.

3. Sweep again on batch B, and find the ceiling. For k in range(1, 21), fit PCR and PLS on all of batch A and score on batch B. Plot these two curves together with the two centered batch-A curves from Part B task 1: four curves, one axes.

Report the k that maximizes each batch-B curve and the best score each reaches, then the best batch-B \(r^2\) achieved by any model you have fit, including OLS, beside the best batch-A test score from Part B.

In two sentences: say whether adding components closes the gap between the two batches, and what that implies about where the remaining error comes from.

4. Decide. You have to put one model into service on the column. State the method, the number of components, and what you would fit it on. Defend it in two sentences using your own numbers from tasks 2 and 3, and name one thing you would change about how the calibration data was collected so that the next version transfers better.

There is more than one defensible answer to task 4. The reasoning and the numbers are what is graded, not the model you land on.

# YOUR CODE HERE
# X_B and y_B were loaded in the setup cell.

Problems: High-dimensional Regression

Contents

Problems: High-dimensional Regression#

Setup#