Class 3/Lab1: Subset Selection and Shrinkage

Live in-browser R — edit and re-run any cell

Published

September 11, 2026

Introduction

In this lab, we will:

Fit a linear regression baseline and interpret its coefficients.
Compare best subset, forward stepwise, and backward stepwise selection.
Fit ridge, lasso, and elastic net, and see how they treat correlated predictors differently.
Use partial least squares as a dimension reduction alternative.

Every code cell below runs live in your browser via webR — no local R install needed. Edit any cell and re-run it. Packages and data load once when the page opens; give it a few seconds on first load.

We use the caret package throughout. Every model below is fit with the same train() call and the same cross-validation settings, so you can compare methods by changing one argument instead of learning four different interfaces.

The data: prostate cancer

The dataset has 97 men who were about to undergo radical prostatectomy. The outcome is logPSA, the log of prostate-specific antigen — the blood marker used to monitor prostate cancer.

The clinical question is one a biostatistician is often asked: which tumor characteristics explain PSA, and can we predict it from a smaller set of measurements? This matters in practice because the predictors differ in cost. Age is free. Cancer volume and Gleason score require a biopsy.

Variable	Clinical meaning
`lcavol`	log cancer volume — how much tumor is present
`lweight`	log prostate weight
`age`	patient age in years
`lbph`	log benign prostatic hyperplasia — non-cancerous enlargement
`svi`	seminal vesicle invasion (0/1) — tumor has spread beyond the prostate
`lcp`	log capsular penetration — tumor breaching the prostate capsule
`gleason`	Gleason score (6–9) — biopsy grade, how aggressive the cells look
`pgg45`	percent of biopsy tissue that is Gleason grade 4 or 5

Note that svi is binary — only 21 of 97 men had seminal vesicle invasion. Keep that in mind when we interpret its coefficient.

One shared cross-validation setup

Every model in this lab uses the same resampling scheme, defined once:

On the choice of 10-fold CV

With only 97 observations, a single 10-fold split is noisy — each held-out fold has about 10 patients, so the estimated error moves around depending on the random split. Results can shift if you change the seed.

For a more stable estimate, use repeated cross-validation:

ctrl_repeated <- trainControl(method = "repeatedcv", number = 10, repeats = 5)

This runs the whole 10-fold procedure 5 times with different splits and averages the results. We use plain 10-fold here to keep the lab fast in the browser. Try swapping in ctrl_repeated at the end and see how much the selected models change.

Part 1: Linear regression baseline

Start with the single strongest predictor. Cancer volume has the highest correlation with PSA, which is what a clinician would expect — more tumor produces more PSA.

Now fit the full model with all eight predictors and compare.

Adding seven predictors improves cross-validated RMSE from about 0.78 to about 0.74. The gain is small compared to the number of parameters we spent to get it. That gap is what the rest of the lab is about.

A warning sign in the coefficients

Look at lcp (capsular penetration). Its coefficient is negative. Now check its correlation with the outcome on its own:

Both variables flip sign. On its own, capsular penetration is strongly and positively associated with PSA (r = 0.55). In the full model its coefficient is negative.

What you would tell the clinician: capsular penetration does raise PSA. But it travels closely with cancer volume (r = 0.68) and seminal vesicle invasion (r = 0.67). Once the model already knows those two things, lcp has nothing left to explain, and its coefficient is determined by noise. The sign flip is a symptom of correlated predictors, not a claim that capsular penetration protects against high PSA.

This is exactly the instability that shrinkage methods are built to control. Note that gleason and pgg45 are correlated at 0.75 — we return to that pair in Part 3.

Part 2: Subset selection

Rather than keeping all eight predictors, choose a subset. We compare three search strategies, all through the same train() interface:

leapSeq — best subset (exhaustive search)
leapForward — forward stepwise
leapBackward — backward stepwise

The tuning parameter nvmax is the number of predictors allowed in the model.

All three strategies agree: about three predictors. For a small predictor set like this one the three methods often agree, but that is not guaranteed. Forward and backward selection are greedy — they commit to one variable at a time and can miss the subset that exhaustive search would find.

Which variables enter at each size?

caret reports cross-validated error, but it does not show the entry order. For that we call regsubsets() directly.

Each row is a model size; each * marks an included predictor. Read down the first column: lcavol enters first and never leaves. Then lweight, then svi. These are the same three predictors the cross-validation chose.

Choosing a model size by information criteria

Cross-validation is one way to pick the size. The classical alternative is a penalized criterion computed on the training fit.

Interpretation: RSS always decreases as you add predictors, so it always picks the largest model and cannot be used to choose a size. That is why the other three criteria exist — each penalizes complexity differently.

The criteria disagree here: BIC picks 3, Cp picks 5, adjusted \(R^2\) picks 7. That disagreement is the lesson. BIC penalizes complexity most heavily, so it gives the smallest model. There is no single correct answer; the choice depends on whether you want a parsimonious model to explain to a clinician or the best possible predictive fit.

Part 3: Ridge, lasso, and elastic net

Subset selection makes a hard in-or-out decision for each predictor. Shrinkage methods instead pull all coefficients toward zero continuously, trading a small amount of bias for a large reduction in variance.

All three methods are the same glmnet model with a different mixing parameter \(\alpha\):

\(\alpha = 0\) — ridge, L2 penalty, shrinks coefficients but keeps all of them
\(\alpha = 1\) — lasso, L1 penalty, can set coefficients to exactly zero
\(0 < \alpha < 1\) — elastic net, a blend of both

caret tunes \(\alpha\) and \(\lambda\) together on one grid. We standardize the predictors with preProcess so the penalty applies fairly to variables measured on different scales.

How the penalty changes what is kept

Fit each \(\alpha\) separately, letting each one choose its own best \(\lambda\), and compare the coefficients side by side.

Interpretation. Ridge (\(\alpha = 0\)) keeps all eight predictors. As \(\alpha\) increases, more coefficients are driven to exactly zero. This is the practical difference between the two penalties: the L1 penalty has corners on the coordinate axes, so the solution can land exactly on zero. The L2 penalty is smooth and shrinks coefficients toward zero without reaching it.

Note that lcavol, lweight, and svi keep large coefficients at every value of \(\alpha\). These are the same three predictors subset selection chose. Tumor burden and spread beyond the prostate drive PSA, and every method agrees.

The grouping effect: `gleason` and `pgg45`

Look at the gleason row of the coefficient table above as \(\alpha\) increases.

These two variables measure the same clinical concept. Both come from the same biopsy and both grade how aggressive the tumor cells look — gleason is the summary score, pgg45 is the percentage of high-grade tissue. They correlate at 0.75.

What you would tell the clinician. Ridge keeps both markers and splits the credit between them (0.050 and 0.085). As \(\alpha\) rises, gleason shrinks steadily — 0.022, then 0.010, then 0.0001 — until the lasso sets it to zero, while pgg45 stays near 0.06 throughout. If the urologist asks why the model discarded the Gleason score, the answer is that pgg45 already carries the same information, and the lasso was forced to choose one of the two. Tumor grade still matters for PSA — the model simply measures it through the other variable.

The important caution is that this choice is driven by noise. With a different sample of 97 patients, the lasso might keep gleason and drop pgg45 instead. When you read a paper reporting a lasso-selected biomarker panel, this is the reason to be careful about which specific markers made the list.

Elastic net sits between the two extremes. It shrinks correlated predictors toward each other instead of forcing a hard choice, so groups of related markers tend to stay in or drop out together. This is the behavior that makes it a sensible default when predictors come in clinically related clusters.

lcp is set to zero by every method with \(\alpha > 0\). The unstable, sign-flipped coefficient from the OLS fit is removed by the penalty. This is shrinkage doing the job it was designed for.

Part 4: Partial least squares

Subset selection and shrinkage both work on the original predictors. Dimension reduction takes a different route: build a small number of new variables as combinations of the originals, then regress on those.

Partial least squares constructs components that explain variance in the predictors and correlate with the outcome.

Interpretation: most of the predictive signal is captured by about three components. Adding more components does not improve cross-validated error, because the later components describe directions in the predictor space that carry little information about PSA.

Three components, three predictors from subset selection, three large coefficients under the lasso. All three approaches point at the same conclusion — this problem has roughly three dimensions of useful signal.

The tradeoff is interpretability. A PLS component is a weighted combination of all eight predictors, so you cannot tell the clinician “drop the Gleason score from the panel.” You can only say “this combination predicts well.”

Reading the first component: it loads heavily on lcavol, lcp, and svi — all measures of tumor burden and spread. The component has a clinical meaning even though it is not a single measurement.

Part 5: Putting it together

Compare every method on the same cross-validation folds.

Interpretation. The methods land close together on this dataset. That is a realistic result — with 97 observations and 8 predictors, no method has much room to separate itself from the others, and the differences are within the noise of the cross-validation.

What differs is the kind of model each produces:

Method	Predictors used	Can explain to a clinician?
OLS	all 8, two with unstable signs	coefficients are hard to defend
Best subset	3	yes — a short, concrete list
Lasso	6 — drops `lcp` and `gleason`	yes, but the choice among correlated markers is arbitrary
Elastic net	7 — `lcp` dropped, `gleason` shrunk toward zero	yes, with grouped markers kept together
PLS	all 8, combined into 3 components	predicts well, individual predictors are not separable

The counts above come from the coefficient table in Part 3. Note that the lasso here drops two predictors rather than reducing to the three that subset selection chose — at its cross-validated \(\lambda\), the penalty is mild enough to retain small coefficients on age, lbph, and pgg45.

Choosing among them is a decision about what the model is for, not only about which has the lowest RMSE.

A caveat worth stating

This dataset has 97 patients. The three-predictor consensus is a finding within this sample, and the cross-validated differences between methods are small relative to their variability. Treat this lab as an illustration of how the methods behave, not as clinical evidence about PSA.

Try re-running the shrinkage section with set.seed(1) instead of set.seed(2155) and see which conclusions hold and which move.