12. Multicollinearity, variable selection, and validation

Measuring body fat well is a nuisance. The accurate way is to weigh a person underwater and work out their density, which needs a tank, a technician, and a cooperative subject. The cheap way is a tape measure. So a natural question follows: can a handful of tape measurements (waist, chest, thigh, wrist, and the rest) predict the percent body fat that the underwater weighing gives? A 1985 study recorded both for 252 men, and that dataset, fat.csv, is the one this chapter works with from start to finish. (This is the larger 252-man study with thirteen measurements, not the three-measurement, twenty-man body-fat table you met for extra sums of squares in 8.5 Added-variable plots; the two share a “thigh” column but are different data.)

There is an obvious plan: regress body fat on all thirteen measurements and read off the coefficients. When you do, the model explains about 75 percent of the variation, which sounds like a success. Then you look at the coefficients and something is wrong. The coefficient on weight is negative, as if heavier men had less fat once you hold the other measures fixed. Several predictors that everyone knows matter come out statistically insignificant. Drop one measurement and the others lurch. The problem is that the tape measurements are all measuring roughly the same thing, body size, so they are tangled together. Figure 1 shows just how tangled.

This chapter is about living with too many predictors, most of them redundant. Three questions organize it. First, how do we diagnose and describe the damage that correlated predictors do, using the variance inflation factor? Second, when several models are on the table, how do we choose, using criteria like Mallows $C_p$ and cross-validation that reward prediction rather than mere fit? Third, what modern tools (the lasso and ridge regression) sidestep the whole selection problem by shrinking coefficients instead of deleting them? Along the way you will meet the most important warning in applied statistics: a model chosen to fit your data will always look better on that data than it deserves.

12.1 Multicollinearity and the variance inflation factor¶

Intuition¶

Chapters 7 and 8 built a multiple regression and read each coefficient as a partial slope, the effect of one predictor with the others held fixed (8.5 Added-variable plots). That reading quietly assumed the predictors carried separate information. The vignette showed the damage when they do not; now we name the cause, multicollinearity (Definition 12.1).

A little of it is harmless and normal. A lot is trouble, because the model cannot tell the correlated predictors apart. Imagine two predictors that move almost in lockstep, like weight and hip circumference. The data show what happens to body fat when both rise together, but they show almost nothing about what happens when one rises and the other stays fixed, because that combination barely occurs. So the model has no firm basis for splitting the shared effect between them. It can give weight a big positive coefficient and hip a big negative one, or the reverse, and fit the data about equally well either way. Figure 2 draws the picture.

The variance inflation factor (Definition 12.2) measures exactly how much the variance is inflated.

The full model, and why its coefficients misbehave¶

Fit the full model first, so the symptoms are concrete. The response is brozek, percent body fat by Brozek’s formula, and the thirteen predictors are the tape and scale measurements. Recall from 4.2 Correlation and the regression slope that for a single predictor, $R^2$ is just the squared correlation; with many predictors $R^2$ keeps climbing as you add columns, whether or not they help, so a high $R^2$ alone proves nothing.

You saw a small version of this in 8.5 Added-variable plots, where a twenty-man body-fat table had a triceps and a thigh measurement so correlated that their coefficients turned unstable. Here the same disease appears at full scale, with thirteen tangled predictors instead of two. A quick look at the correlations among the size measurements confirms the diagnosis.

round(cor(fat[, c("weight", "chest", "abdom", "hip", "thigh")]), 3)

       weight chest abdom   hip thigh
weight  1.000 0.894 0.888 0.941 0.869
chest   0.894 1.000 0.916 0.829 0.730
abdom   0.888 0.916 1.000 0.874 0.767
hip     0.941 0.829 0.874 1.000 0.896
thigh   0.869 0.730 0.767 0.896 1.000

print(fat[["weight", "chest", "abdom", "hip", "thigh"]].corr().round(3))

        weight  chest  abdom    hip  thigh
weight   1.000  0.894  0.888  0.941  0.869
chest    0.894  1.000  0.916  0.829  0.730
abdom    0.888  0.916  1.000  0.874  0.767
hip      0.941  0.829  0.874  1.000  0.896
thigh    0.869  0.730  0.767  0.896  1.000

Every pair correlates above 0.7, and weight with hip reaches 0.941. These columns are near-duplicates of one another.

Formula¶

In words: regress predictor $k$ on the remaining predictors, see how well they predict it, and $\mathrm{VIF}_k$ is one over the leftover fraction. If $X_k$ is unrelated to the others, $R_k^2 = 0$ and $\mathrm{VIF}_k = 1$. If the others explain 90 percent of it, $R_k^2 = 0.9$ and $\mathrm{VIF}_k = 10$. A common rule of thumb flags $\mathrm{VIF}_k > 10$ (equivalently $R_k^2 > 0.9$) as serious, though the cutoff is a convention, not a law.

The shape of that formula is worth a look, because it explains why a little collinearity is harmless and a lot is a disaster. Figure 3 plots the VIF against $R_k^2$. The curve is nearly flat while a predictor is mostly its own, then bends sharply upward as the others close in on reproducing it. Weight, which the other twelve measurements explain 97 percent of, sits far up the steep part.

The name is literal. The VIF is the factor by which the variance of $b_k$ is multiplied, compared with the variance it would have if $X_k$ were uncorrelated with the other predictors:

In words: the variance of a coefficient is the usual $\sigma^2 / S_{kk}$ (small when the predictor is spread wide) times the inflation factor. Collinearity enters only through that last factor.

Derivation (the VIF from the auxiliary regression)¶

Proof. Fix a predictor $X_k$. By the added-variable construction of 8.5 Added-variable plots, the least-squares coefficient $b_k$ equals the slope from a simple regression of $Y$ on the part of $X_k$ that the other predictors do not explain. Make that part concrete: run the auxiliary regression of $X_k$ on all the other predictors, and call its residuals $\tilde{X}_{ik}$. These residuals are $X_k$ with the other predictors’ information removed.

The added-variable result says $b_k = \sum_i \tilde{X}_{ik} Y_i / \sum_i \tilde{X}_{ik}^2$, a simple-regression slope with the $\tilde{X}_{ik}$ as the predictor. From 2.5 Sampling behavior and the Gauss-Markov theorem, a simple-regression slope $\sum k_i Y_i$ with weights $k_i = \tilde{X}_{ik} / \sum_j \tilde{X}_{jk}^2$ has variance $\sigma^2 \sum_i k_i^2 = \sigma^2 / \sum_i \tilde{X}_{ik}^2$. So

Now identify the denominator. The auxiliary regression has total sum of squares $\sum_i (X_{ik} - \bar{X}_k)^2 = S_{kk}$ and error sum of squares $\sum_i \tilde{X}_{ik}^2$. By the definition of $R^2$ applied to that auxiliary regression,

Substitute:

The first factor $\sigma^2 / S_{kk}$ is what the variance would be if $X_k$ were free of the others ($R_k^2 = 0$). The second factor $1/(1 - R_k^2) = \mathrm{VIF}_k$ is the price of collinearity. When the other predictors nearly reproduce $X_k$, $R_k^2 \to 1$ and the variance explodes. $\blacksquare$

R and Python¶

R’s car package computes every VIF at once. The formula above lets you check any one of them by hand from a single auxiliary regression.

Remedies¶

There are four honest responses to high VIFs. First, do nothing, if you only care about prediction inside the range of the data: collinearity does not bias $\hat{Y}$ or widen its interval much, so a predictive model can carry redundant predictors safely. Second, drop predictors: if weight, hip, and abdomen say nearly the same thing, keep the one you can measure best and discard the rest. Third, combine them into a single index (an average, or a principal component) that captures the shared size signal. Fourth, shrink the coefficients with ridge regression (Section 12.5), which tolerates collinearity by design. What you must not do is read a collinear coefficient as a real causal effect and act on its sign.

Reading about inflated variance is one thing; here is the same disease with a dial on it, so you can decide for yourself how much correlation is too much.

12.2 Choosing among models: selection criteria¶

Intuition¶

The full model has redundant predictors; a smaller model might predict just as well and be easier to trust. But which smaller model? With thirteen predictors there are $2^{13} = 8192$ possible subsets. We need a score that ranks them. The obvious score, $R^2$, is useless for this, because $R^2$ never decreases when you add a predictor: the biggest model always wins, even if the added predictors are noise. Recall from 4.2 Correlation and the regression slope that $R^2$ is squared correlation, a measure of fit, not of prediction. A good selection criterion must instead reward fit while charging a penalty for each extra parameter, so that a predictor earns its place only if it pays for itself.

Four criteria do this, each in its own currency. Adjusted $R^2$ discounts $R^2$ by the degrees of freedom spent. Mallows $C_p$ estimates the total prediction error and compares it to the number of parameters. AIC and BIC trade goodness of fit against a size penalty from information theory. PRESS predicts each case from a model fit without it. They usually agree on the broad message (drop the deadweight) while disagreeing on the exact size, which is itself a lesson: there is no single correct model.

Formula¶

Let a candidate model have $p$ parameters (predictors plus intercept), error sum of squares $\mathrm{SSE}_p = \sum_i (Y_i - \hat{Y}_i)^2$, and $\mathrm{SSTO} = \sum_i (Y_i - \bar{Y})^2$. Let $\hat{\sigma}^2$ be the mean square error of the largest (full) model, our best estimate of the noise variance $\sigma^2$. The four criteria are

• $R^2_{a,p}$ (adjusted $R^2$): like $R^2$ but with $\mathrm{SSE}$ and $\mathrm{SSTO}$ divided by their degrees of freedom. Adding a useless predictor lowers it. Bigger is better.

• $C_p$ (Mallows): small is better; a model with no bias has $C_p \approx p$, so you look for models with low $C_p$ that also sit near the line $C_p = p$.

• $\mathrm{AIC}_p$, $\mathrm{BIC}_p$: fit term plus penalty; smaller is better. BIC’s penalty $p \ln n$ is heavier than AIC’s $2p$ whenever $n > 7$, so BIC prefers smaller models.

That last claim is easy to see in a picture. Figure 6 draws the two size penalties as straight lines. With $n = 252$ men, BIC charges about $\ln 252 \approx 5.5$ for each extra parameter while AIC charges only 2, so the BIC line climbs almost three times as steeply. A predictor has to improve the fit more to pay off its heavier BIC toll, which is exactly why BIC lands on smaller models.

The fifth criterion, PRESS, gets its own subsection because it measures something different: honest out-of-sample error.

where $\hat{Y}_{(i)}$ is the prediction for case $i$ from the model refit with case $i$ removed. In words: hide each observation, predict it from the rest, and add up the squared misses. Smaller is better, and unlike $\mathrm{SSE}$, PRESS cannot be driven down just by adding predictors, because a junk predictor helps the fit but hurts the left-out prediction.

Derivation (Mallows $C_p$)¶

Proof. We want a criterion that is small when a model’s fitted values $\hat{Y}_i$ are close to the true means $\mu_i = E\{Y_i\}$, across all $n$ cases. Define the target

the total mean squared error of the fitted values, standardized by $\sigma^2$. Split each term into variance and squared bias with the identity $E\{W^2\} = \operatorname{Var}\{W\} + (E\{W\})^2$ applied to $W = \hat{Y}_i - \mu_i$:

For the variance term, the fitted values are $\hat{\mathbf{Y}} = \mathbf{H}\mathbf{Y}$ with the candidate model’s hat matrix $\mathbf{H}$, and $\sum_i \operatorname{Var}\{\hat{Y}_i\} = \sigma^2\,\mathrm{tr}(\mathbf{H}) = p\sigma^2$, using $\mathrm{tr}(\mathbf{H}) = p$ from 7.3 The hat matrix. So

Now connect $B$ to something we can compute, the expected error sum of squares. Write $\mathbf{Y} = \boldsymbol{\mu} + \boldsymbol{\varepsilon}$ and $\mathrm{SSE}_p = \|(\mathbf{I} - \mathbf{H})\mathbf{Y}\|^2$. Then

because $(\mathbf{I}-\mathbf{H})\boldsymbol{\mu}$ has $i$-th entry $\mu_i - E\{\hat{Y}_i\}$ so its squared length is $B$, and $\mathbf{I}-\mathbf{H}$ is idempotent with trace $n - p$ (again 7.3 The hat matrix). Solve for the bias: $B = E\{\mathrm{SSE}_p\} - (n - p)\sigma^2$. Substitute into $\Gamma_p$:

Replace the unknown $E\{\mathrm{SSE}_p\}$ by the observed $\mathrm{SSE}_p$ and the unknown $\sigma^2$ by the full model’s $\hat{\sigma}^2$, and you have Mallows’ statistic:

Two consequences fall out. If the candidate model has no bias ($B = 0$, meaning it contains all the predictors that matter), then $E\{\mathrm{SSE}_p\} = (n - p)\sigma^2$ and $\Gamma_p = (n - p) - (n - 2p) = p$. So unbiased models have $C_p \approx p$: you look for models where $C_p$ is both small and close to $p$. A model that omits an important predictor has $B > 0$, which pushes $C_p$ well above $p$. $\blacksquare$

R and Python¶

R’s leaps package searches all 2¹³ subsets and reports the best model of each size.

Deriving and checking PRESS¶

PRESS looks expensive: refitting the model $n$ times, once per deleted case, is 252 fits per model here. A beautiful identity makes it free. You can get every deleted prediction from a single fit of the full data, using the ordinary residual and the leverage value.

Proof. Write the design matrix $\mathbf{X}$ with $i$-th row $\mathbf{x}_i'$, the full-data estimate $\mathbf{b} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{Y}$, the fitted value $\hat{Y}_i = \mathbf{x}_i'\mathbf{b}$, the residual $e_i = Y_i - \hat{Y}_i$, and the leverage value $h_{ii} = \mathbf{x}_i'(\mathbf{X}'\mathbf{X})^{-1}\mathbf{x}_i$ from 7.3 The hat matrix (the diagonal of the hat matrix). Let $\mathbf{b}_{(i)}$ be the estimate from the data with case $i$ removed, and $\hat{Y}_{(i)} = \mathbf{x}_i'\mathbf{b}_{(i)}$ the deleted prediction. We show

Deleting case $i$ changes the cross-product matrix and vector to $\mathbf{X}'\mathbf{X} - \mathbf{x}_i\mathbf{x}_i'$ and $\mathbf{X}'\mathbf{Y} - \mathbf{x}_i Y_i$. The Sherman-Morrison inverse-update identity (derived and checked in 9.3 Influence: which points actually change the fit) gives

Write $\mathbf{u} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{x}_i$, so that $\mathbf{x}_i'\mathbf{u} = h_{ii}$ and $\mathbf{x}_i'(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{Y} = \mathbf{x}_i'\mathbf{b} = \hat{Y}_i$. Multiply the updated inverse by $\mathbf{X}'\mathbf{Y} - \mathbf{x}_i Y_i$:

where the numerator collapsed to $\hat{Y}_i - Y_i = -e_i$. Now form the deleted prediction:

So each deleted residual is just the ordinary residual divided by $1 - h_{ii}$, and

High-leverage points (large $h_{ii}$) get their residuals magnified, which is right: those are the cases the model props itself up on, and it predicts them worst when they are removed. $\blacksquare$

Before you accept any one criterion’s verdict, walk the whole ladder of model sizes yourself and watch the four scores pull in different directions.

12.3 The dangers of automatic selection¶

Intuition¶

It is tempting to automate the search: let the computer add and drop predictors by their $p$-values until nothing improves. This is stepwise selection (Definition 12.7), and it is available in one line in every statistics package. It is also the source of more bad science than almost any other routine procedure, and you should understand exactly why before you ever use it.

Two things go wrong. The first is overfitting (Definition 12.8): a procedure that chases the best fit on your sample will seize on accidents of that sample, patterns that are noise and will not repeat. The second, subtler and worse, is post-selection inference (Definition 12.9): once you have chosen a model because it fit well, the $p$-values, confidence intervals, and $F$ tests that software prints for that model are wrong. They were derived assuming the model was fixed in advance, before you looked at the data. Choosing the model by looking at the data breaks that assumption, and the reported significance becomes fiction.

A demonstration with pure noise¶

The cleanest way to see the disaster is to run selection on data with no signal at all. Take a response and forty predictors that are all independent random noise, unrelated to the response and to each other. Honest testing should almost never find anything. Selection finds plenty.

set.seed(4210)
noise <- data.frame(y = rnorm(100), matrix(rnorm(100 * 40), 100, 40))
names(noise)[-1] <- paste0("x", 1:40)
pvals <- sapply(1:40, function(j)
  summary(lm(noise$y ~ noise[[paste0("x", j)]]))$coefficients[2, 4])
best5 <- order(pvals)[1:5]
selected <- lm(y ~ ., data = noise[, c("y", paste0("x", best5))])
fstat <- summary(selected)$fstatistic
c(R2 = summary(selected)$r.squared,
  overall_F_pvalue = as.numeric(pf(fstat[1], fstat[2], fstat[3], lower.tail = FALSE)))

              R2 overall_F_pvalue
      0.13497399       0.01663716

We picked the five predictors with the smallest individual $p$-values out of forty pure-noise columns, fit a model on just those five, and got an overall $F$-test $p$-value of 0.017. By the usual rule, that model is “significant at the 0.05 level,” and its $R^2 = 0.135$ looks respectable. Every bit of it is an illusion: there is no relationship in the data at all. The selection step manufactured the significance by cherry-picking, and the $F$ test, blind to the cherry-picking, reported it as real.

This is not a one-off fluke of the seed. Figure 10 repeats the whole experiment 2000 times. Under honest testing the overall $F$ $p$-value should be uniform on $[0, 1]$, so only 5 percent should fall below 0.05. After selection, 94 percent do.

What to do instead¶

None of this means variable selection is forbidden. It means three things. First, prefer the criterion-based methods of Section 12.2 (compare whole models by $C_p$, BIC, or cross-validated error) over $p$-value-driven stepwise search, and report that you searched. Second, never quote the $p$-values from a selected model as if the model had been fixed in advance; recall from 8.4 The general linear test that the $F$ test’s distribution assumes the comparison was chosen before seeing the data. Third, and most reliably, judge the final model on data it was not chosen with. That is validation, the subject of the next section, and it is the only honest referee for a model you built by searching.

12.4 Validation: train, test, and cross-validate¶

Intuition¶

The residuals from a fitted model tell you how well it fits the data it was built on, which is always flattering. To know how well a model predicts, you must test it on data it has never seen. The idea echoes the out-of-sample thinking behind the bootstrap in 5.4 The bootstrap for regression and the prediction intervals of 3.5 Prediction interval for a new observation: honest error estimates come from treating some data as genuinely new.

The simplest version splits the data once: a training set (Definition 12.10) to fit the model, and a test set, held back, to measure its error. The split must be made before any modeling, and the test set must never influence a single choice. A more efficient version, $k$-fold cross-validation (Definition 12.11), splits the data into $k$ equal folds, then rotates: each fold takes a turn as the test set while the other $k - 1$ folds train the model. Averaging the $k$ test errors uses every case for testing exactly once, so it wastes no data. Figure 11 shows the rotation.

We measure error with the root mean squared error (Definition 12.12), the typical size of a prediction miss, in the units of the response (percentage points of body fat). A smaller one means the model is closer, on average, to the truth for a new man.

The single train/test split¶

Split the 252 men with a fixed seed so the analysis reproduces, fit on the training set, and compare training error to test error.

$k$-fold cross-validation, written by hand¶

A single split is noisy: a lucky or unlucky test set can mislead. Cross-validation averages over $k$ splits. It is worth writing the loop yourself once, so the mechanism holds no mystery. Assign each case a fold number, then loop: train on the other folds, predict the held-out fold, collect the squared errors.

Example 12.5 reported the test RMSE from one split. This widget reruns that split a thousand times, so you can see all the other answers it could just as easily have given you.

12.5 A look at shrinkage: ridge and lasso¶

Intuition¶

Selection is a blunt instrument: a predictor is either in (full coefficient) or out (zero). Shrinkage (Definition 12.13) offers a gentler alternative. Instead of deleting predictors, it keeps them all but pulls their coefficients toward zero, trading a little bias for a large drop in variance. When predictors are collinear, ordinary least squares can produce enormous coefficients that cancel each other (the negative weight from Example 12.1). Shrinkage refuses to let coefficients grow that large, which stabilizes them.

Formula¶

Ordinary least squares minimizes the residual sum of squares. Ridge (Definition 12.14) and lasso (Definition 12.15) add a penalty on the size of the coefficients.

• The first term is the usual least-squares loss; it wants a good fit.

• The second term is the penalty; it wants small coefficients. The sum runs over the $q$ slope coefficients (not the intercept), and predictors are standardized first so the penalty treats them fairly.

• $\lambda \ge 0$ is the tuning parameter that sets the strength of the penalty. At $\lambda = 0$ both reduce to least squares; as $\lambda \to \infty$ all coefficients are crushed to zero. We choose $\lambda$ by cross-validation.

Ridge uses the squared penalty ($\sum \beta_k^2$, the $L_2$ norm); lasso uses the absolute-value penalty ($\sum |\beta_k|$, the $L_1$ norm). The difference sounds minor and is not. Ridge shrinks every coefficient toward zero but rarely to exactly zero. Lasso can set coefficients exactly to zero, so it does variable selection and shrinkage at once. The reason is geometric, shown in Figure 15: the lasso’s constraint region is a diamond with corners on the axes, and the best-fitting contour tends to touch a corner, where a coefficient is zero. Ridge’s constraint region is a smooth disk with no corners.

R and Python¶

The R package glmnet is the standard tool. Its cv.glmnet fits the whole path of $\lambda$ values and picks one by cross-validation, returning lambda.min (the error minimizer) and lambda.1se (the largest $\lambda$ within one standard error of the minimum, a deliberately sparser choice). Predictors are standardized automatically.

The plain summary of shrinkage¶

Ridge and lasso answer the collinearity and selection problems in one move. Ridge keeps every predictor but shrinks the coefficients so no two collinear predictors can blow up against each other; it is the natural tool when you believe many predictors each contribute a little. Lasso shrinks and selects at once, zeroing the coefficients it does not need, so it is the natural tool when you believe only a few predictors matter and you want the model to find them. Both replace the discrete, unstable question “which predictors are in?” with a smooth, tunable dial $\lambda$, set by cross-validation. That is why shrinkage, not stepwise selection, is the modern default when prediction is the goal.

12.6 Chapter summary¶

You can now handle a regression with more predictors than you can trust. You diagnosed multicollinearity with the variance inflation factor and derived it from an auxiliary regression; you compared candidate models with adjusted $R^2$, Mallows $C_p$, AIC, BIC, and PRESS, deriving $C_p$ and proving the deleted-residual shortcut that makes leave-one-out free; you saw, in a pure-noise demonstration, why stepwise selection overfits and why post-selection $p$-values lie; you wrote a $k$-fold cross-validation loop by hand and used a held-out test set to judge a model honestly; and you described and ran ridge and lasso, reading shrinkage as a bias-variance trade made by a tunable penalty. On the body-fat data the eight-predictor candidate beat the full model on every honest measure: lower $C_p$, lower PRESS, and lower cross-validated error.

Figure 17 puts the whole chapter on one page: diagnose the redundancy, let your goal decide the remedy, choose a model by a prediction-minded score, and confirm the winner on data it never saw.

Key results at a glance

Key terms. Multicollinearity, variance inflation factor, auxiliary regression, adjusted $R^2$, Mallows $C_p$, AIC, BIC, PRESS, deleted residual, stepwise selection, overfitting, post-selection inference, training and test sets, $k$-fold cross-validation, root mean squared error, bias-variance trade-off, shrinkage, ridge regression, lasso, tuning parameter $\lambda$.

You should now be able to

• Diagnose multicollinearity with the VIF and derive it from an auxiliary regression (Theorem 12.4).

• Explain what multicollinearity does and does not harm, and list the standard remedies.

• Compute and compare adjusted $R^2$, Mallows $C_p$, AIC, BIC, and PRESS.

• Derive Mallows $C_p$ (Theorem 12.5) and prove the PRESS deleted-residual shortcut (Theorem 12.6).

• Explain, bluntly, why stepwise selection and post-selection inference mislead.

• Implement $k$-fold cross-validation by hand in R and Python.

• Describe ridge and lasso shrinkage and demonstrate them on the body-fat data.

Where this fits. This chapter is the FIT and CHECK stages of The modeling workflow done together and repeatedly: you fit many models, check each for collinearity and prediction error, and iterate toward one you can defend. Chapters 8 through 11 gave you a single multiple regression and the tools to interpret and diagnose it; this chapter faced the harder situation of many competing models and taught you to choose among them by prediction rather than fit. The validation mindset carries straight into 13. Logistic regression, where the same held-out thinking becomes the classification metrics (accuracy, ROC) that judge a logistic model, and the shrinkage and criterion ideas reappear wherever a model has more parameters than the data can comfortably support.

12.7 Frequently asked questions¶

Q1. Does multicollinearity bias my coefficients? No. Least squares stays unbiased under collinearity: on average the coefficients are still correct. What collinearity does is inflate their variance, so any single sample’s estimates can be wildly off, with flipped signs and huge standard errors. The average is fine; the individual estimate is unreliable.

Q2. If collinearity does not hurt prediction, why care about it? It depends on your goal. If you only want to predict body fat for men like those in the data, a collinear model predicts fine and you can ignore the VIFs. If you want to interpret coefficients (to say which measurement matters and by how much), collinearity makes those coefficients meaningless, and you must fix it by dropping, combining, or shrinking predictors.

Q3. Which selection criterion should I use? For interpretation, BIC and its heavier penalty give you a lean, defensible model. For prediction, cross-validated error is the most direct answer, and Mallows $C_p$ is a fast, well-motivated stand-in. Do not agonize over small differences: report that the criteria pointed to a range of models (here, four to eight predictors) and pick within that range for a reason you can state.

Q4. Why is the $C_p$ of the full model always equal to the number of its parameters? By construction. $C_p = \mathrm{SSE}_p/\hat{\sigma}^2 - (n - 2p)$, and for the full model $\hat{\sigma}^2 = \mathrm{SSE}_{\text{full}}/(n - p_{\text{full}})$, so $\mathrm{SSE}_{\text{full}}/\hat{\sigma}^2 = n - p_{\text{full}}$ and $C_p = (n - p_{\text{full}})

• (n - 2p_{\text{full}}) = p_{\text{full}}$. The full model is the reference point, so its$C_p$lands exactly on the$C_p = p$ line.

Q5. Is a lower PRESS always better? As a comparison of models on the same data, yes: lower PRESS means better leave-one-out prediction. But PRESS is still an in-sample estimate in the sense that it uses the same 252 men to build and to test (just never the same man for both at once), so it can be optimistic if the whole dataset is unusual. A truly held-out test set, or $k$-fold cross-validation repeated with different seeds, is a stronger check.

Q6. Can I just trust the model my software’s step() function returns? No. Stepwise selection overfits and its reported $p$-values are invalid (Section 12.3). If you use it, treat its output as a candidate to be validated on held-out data, never as a finished model, and never quote its coefficients’ significance as if the model had been chosen in advance.

Q7. Lasso and stepwise both drop predictors. What is the difference? Stepwise makes hard, unstable in-or-out decisions by $p$-value and reports dishonest inference. Lasso makes the decision continuously through a single penalty $\lambda$ chosen by cross-validation, shrinking coefficients smoothly and turning them off only as the penalty grows. Lasso is more stable, has a coherent bias-variance justification, and is built for prediction. When in doubt, prefer it.

Q8. Why standardize predictors before ridge or lasso? The penalty $\sum \beta_k^2$ or $\sum |\beta_k|$ adds up coefficients measured in each predictor’s own units. Without standardizing, a predictor measured in small units (and so carrying a large coefficient) would be penalized more than an identical predictor in large units. Standardizing puts every predictor on the same scale so the penalty is fair. Both glmnet and the scaled scikit-learn fit here do this.

12.8 Practice problems¶

1. (A) In one sentence each, say what multicollinearity does to the bias of $b_k$, to the variance of $b_k$, and to predictions $\hat{Y}$ made inside the data range. Then name two of the standard remedies for a predictor with a high VIF.

2. (A) A predictor has $\mathrm{VIF} = 5$. What is the auxiliary $R_k^2$, and how many times larger is $\operatorname{Var}\{b_k\}$ than it would be if the predictor were uncorrelated with the others?

3. (A) Explain why $R^2$ cannot be used to choose between a model with 5 predictors and one with 8 predictors, but adjusted $R^2$ can.

4. (A) State the rule of thumb that connects a “good” model’s $C_p$ to its number of parameters $p$, and say what a $C_p$ far above $p$ indicates.

5. (A) Why does BIC tend to select smaller models than AIC? Point to the specific term that differs.

6. (A) In plain words, what does the PRESS statistic measure that ordinary $\mathrm{SSE}$ does not?

7. (A) A student selects predictors by stepwise regression, then reports the model’s $p$-values as evidence it is correct. Name the fallacy and explain it in two sentences.

8. (A) Describe the difference between what ridge and lasso do to a coefficient, and say which one performs variable selection.

9. (A) A model’s training RMSE is 3.5 and its 10-fold CV RMSE is 4.6. What does the gap tell you, and which number estimates performance on new data?

10. (A) Explain why predictors are standardized before applying a ridge or lasso penalty.

11. (B) Derive $\operatorname{Var}\{b_k\} = \sigma^2/(S_{kk}(1 - R_k^2))$ (Theorem 12.4) starting from the added-variable fact that $b_k$ is the slope of $Y$ on the residuals of $X_k$ regressed on the other predictors. Identify where $\mathrm{VIF}_k$ appears.

12. (B) Show that the auxiliary regression’s error sum of squares equals $S_{kk}(1 - R_k^2)$, and use it to explain why a predictor nearly determined by the others has a nearly zero denominator in its coefficient variance.

13. (B) Derive Mallows $C_p$ (Theorem 12.5) from the standardized total mean squared error $\Gamma_p$, including the steps $\sum_i \operatorname{Var}\{\hat{Y}_i\} = p\sigma^2$ and $E\{\mathrm{SSE}_p\} = B + (n - p)\sigma^2$. Conclude that an unbiased model has $C_p \approx p$.

14. (B) Prove the PRESS deleted-residual shortcut $Y_i - \hat{Y}_{(i)} = e_i/(1 - h_{ii})$ (Theorem 12.6), stating clearly where the Sherman-Morrison inverse update is used.

15. (B) Show that for the full model, $C_p$ equals its number of parameters exactly. (Start from the definition of $C_p$ and the fact that $\hat{\sigma}^2$ is the full model’s MSE.)

16. (B) A high-leverage case has $h_{ii} = 0.9$ and ordinary residual $e_i = 2$. Compute its deleted residual, and explain why leave-one-out prediction magnifies the residuals of high-leverage cases.

17. (B) Explain, using the bias-variance decomposition of prediction error, why a model can have both lower training error and higher test error than a smaller model. Which term of the decomposition rises as the model grows?

18. (B) Show that as $\lambda \to \infty$ the ridge coefficient of a single standardized predictor $\hat{\beta}_{\text{ridge}} = \sum x_i Y_i / (\sum x_i^2 + \lambda)$ tends to zero, and that at $\lambda = 0$ it equals the least-squares slope. (Use the one-predictor ridge solution, which you may derive in Problem 19.)

19. (B) Derive the one-predictor ridge estimator in Problem 18 by setting the derivative of $\sum (Y_i - \beta x_i)^2 + \lambda \beta^2$ to zero, and interpret the denominator as “signal plus penalty.”

20. (C) Fit the full model and reproduce the VIFs. Identify every predictor with $\mathrm{VIF}

    10$ and, for one of them, confirm the VIF by running its auxiliary regression by hand.

21. (C) Run the exhaustive best-subset search (R leaps, or the Python loop) and report the best model of each size with its $C_p$. State the sizes chosen by $C_p$, adjusted $R^2$, and BIC.

22. (C) For the four-predictor best subset (weight, abdomen, forearm, wrist), compute $C_p$, AIC, BIC, and PRESS. Compare each to the eight-predictor candidate and say which model you would report and why.

23. (C) Write a function that computes PRESS from a fitted model using the shortcut, and verify it against a brute-force leave-one-out loop on the candidate model. Report both values.

24. (C) Reproduce the pure-noise post-selection demonstration with a different seed. Report the selected model’s $R^2$ and overall $F$ $p$-value, and explain why they are misleading.

25. (C) Perform a seeded 70/30 train/test split. Fit the full and candidate models on the training set and report both test-set RMSEs. Which model predicts better, and does the order match the training-set fit?

26. (C) Write a 5-fold cross-validation loop by hand and compare the full model, the eight-predictor candidate, and the four-predictor model (weight, abdomen, forearm, wrist) by CV RMSE. Rank them.

27. (C) Fit a lasso with cv.glmnet (R) or LassoCV (Python) to the whole dataset. Report which predictors have nonzero coefficients at the cross-validated penalty, and compare that set to the eight-predictor $C_p$ model.

28. (C) Fit ridge regression across a grid of $\lambda$ values and plot the coefficient of weight against $\log\lambda$. Describe how the coefficient changes, and connect its behavior to the VIF of 33.5 you found for weight.

29. (C) Using the candidate model, plot the deleted residuals $e_i/(1 - h_{ii})$ against the ordinary residuals $e_i$. Which cases move the most, and what property of those cases explains it?

30. (C) Take the 25 percent of men with the largest leverages and refit the candidate model without them. Report how much the coefficients and PRESS change, and discuss whether the model is being propped up by a few cases.

12.9 Exam practice¶

These five questions are written in the style of the course exams: each asks you to explain your reasoning in full sentences, not just to produce a number. Where output is shown, it was produced on the course machine from fat.csv; read the numbers you need off the printout and say which you used. Work each one before opening its model answer.

EP 12.1 (concepts, evaluate a claim). The full body-fat model of Example 12.1 was refit and the weight coefficient and three variance inflation factors were printed.

full <- lm(brozek ~ age + weight + height + neck + chest + abdom + hip +
             thigh + knee + ankle + biceps + forearm + wrist, data = fat)
round(coef(summary(full))["weight", ], 4)
round(car::vif(full)[c("weight", "hip", "abdom")], 2)

  Estimate Std. Error    t value   Pr(>|t|)
   -0.0803     0.0496    -1.6198     0.1066
weight    hip  abdom
 33.51  14.80  11.77

A student reads the negative weight coefficient and concludes: “holding the other twelve measurements fixed, heavier men have less body fat, so we should tell clients that gaining weight lowers their body fat percentage.” Evaluate this claim. Use the variance inflation factor in your answer, and state clearly what the coefficient can and cannot support.

EP 12.2 (interpret output in context). An exhaustive best-subset search was run on a six-predictor version of the model, brozek ~ age + weight + height + abdom + hip + thigh, reporting the best subset of each size.

rs <- regsubsets(brozek ~ age + weight + height + abdom + hip + thigh,
                 data = fat, nvmax = 6)
ss <- summary(rs)
rbind(Cp = round(ss$cp, 2), adjR2 = round(ss$adjr2, 4), BIC = round(ss$bic, 1))
c(min_Cp = which.min(ss$cp), max_adjR2 = which.max(ss$adjr2), min_BIC = which.min(ss$bic))

        [,1]    [,2]    [,3]    [,4]    [,5]    [,6]
Cp     55.07    6.29    3.88    4.40    5.13    7.00
adjR2 0.6608  0.7165  0.7203  0.7208  0.7212  0.7202
BIC  -262.40 -303.10 -302.00 -298.00 -293.70 -288.30
   min_Cp max_adjR2   min_BIC
        3         5         2

The best size-2 model is weight + abdom, the best size-3 model is weight + abdom + thigh, and the best size-5 model is weight + height + abdom + hip + thigh. Interpret this output. State which model size each of the three criteria selects, explain in full sentences why Mallows $C_p$ and BIC land on different sizes, and say which model you would report to a colleague who wants to interpret the coefficients, with your reason.

EP 12.3 (what would change if). For the four-predictor model brozek ~ weight + abdom + forearm + wrist, the error sum of squares is $\mathrm{SSE}_p = 3994.31$ with $n = 252$ and $p = 5$ parameters. The full model’s estimate of the noise variance is $\hat{\sigma}^2 = 15.90$.

cand4 <- lm(brozek ~ weight + abdom + forearm + wrist, data = fat)
c(SSE = sum(resid(cand4)^2), sigma2_full = summary(full)$sigma^2)

       SSE sigma2_full
  3994.311      15.904

First compute $C_p$ for this model. Then answer: what would change if you estimated $\sigma^2$ not from the full model but from the abdomen-only model, whose residual mean square is 20.38? Compute the resulting $C_p$ and explain, in full sentences, why $C_p$ becomes meaningless with that choice.

EP 12.4 (evaluate a claim, interpret output). Three models were compared by their training RMSE and by 5-fold cross-validated RMSE on the full body-fat data: the thirteen-predictor full model, the eight-predictor candidate (age, weight, neck, abdomen, hip, thigh, forearm, wrist), and a two-predictor model (abdomen, weight).

set.seed(4210)
folds <- sample(rep(1:5, length.out = nrow(fat)))
# training RMSE and 5-fold CV RMSE for each model

full   train RMSE 3.876   5-fold CV RMSE 4.097
cand8  train RMSE 3.893   5-fold CV RMSE 4.003
two    train RMSE 4.103   5-fold CV RMSE 4.183

A student argues: “the full model has the lowest training RMSE (3.876), so it predicts best; report the full model.” Evaluate this argument. Say which number estimates how the model will do on a new man and why, explain why the ranking flips between the two columns, and state which model you would report.

EP 12.5 (what would change if). A lasso was fit to the whole body-fat data with the penalty $\lambda$ chosen by cross-validation, and the coefficients were read at two penalties, lambda.min (the error minimizer) and the sparser lambda.1se.

X <- as.matrix(fat[, preds]); Y <- fat$brozek
set.seed(4210)
cvl <- cv.glmnet(X, Y, alpha = 1)
round(cbind(min = coef(cvl, s = "lambda.min")[, 1],
            onese = coef(cvl, s = "lambda.1se")[, 1]), 3)

                min   onese
(Intercept) -13.683  -6.733
age           0.055   0.052
weight       -0.065   0.000
height       -0.081  -0.149
neck         -0.403  -0.137
chest         0.000   0.000
abdom         0.830   0.635
hip          -0.147   0.000
thigh         0.171   0.000
knee          0.000   0.000
ankle         0.086   0.000
biceps        0.114   0.000
forearm       0.387   0.112
wrist        -1.437  -1.265

At lambda.min eleven predictors have nonzero coefficients; at lambda.1se only six do. Explain what happens to the coefficients as $\lambda$ increases from lambda.min to lambda.1se, why the lasso sets some coefficients exactly to zero while ridge regression would not, and what would change if the predictors were not standardized before fitting.

Chapter game¶