ExpertIndependently reviewed, not yet spot-checked by a human

Shrinkage and penalised regression

When the sample is too small, maximum likelihood returns coefficients that are too extreme — the model has learned this particular sample's noise as if it were signal. This page fits one model two hundred times at each of seven development sample sizes to measure how large the overfitting is, and shows that Van Houwelingen's uniform shrinkage factor can almost predict in advance how far the calibration slope will fall. It also covers what separates ridge, lasso and elastic net, and why a penalised coefficient can no longer be read as an effect size.

Overfitting is not “the model is too complex”, it is “the model memorised the noise”

What maximum likelihood estimation does is this: find the coefficients that make the data in front of you as probable as possible. That data contains signal and noise, and the objective function does not distinguish between them — if fitting the noise raises the likelihood a little further, it will fit the noise.

The consequence is very concrete: the coefficients come out too large in absolute value. The model pushes high-risk people higher than they really are and presses low-risk people lower. So on new people the predictions are too extreme — which the calibration page will show up as a calibration slope below 1.

You can make this happen on demand. survival::rotterdam has 2982 patients, which is large enough to manufacture an overfitting scenario: draw a small subset to fit the model, then measure its performance on everyone left over. Those people come from the same cohort and the same measurements, so the gap can only be overfitting — no population difference is mixed in.

The same 8-variable Cox model, seven development sample sizes, 200 repetitions each:

Left: apparent performance, measured on the same patients the model was fitted to, is always better than held-out performance, and the gap widens as the sample shrinks. Right: Van Houwelingen's heuristic shrinkage factor, computable from the development data alone, against the calibration slope actually measured on the held-out patients. The two lines nearly coincide.Plotting script figures/scripts/B5-02-shrinkage.R

Development n	Events	EPV	Apparent C-index	Held-out C-index	Optimism	Held-out calibration slope
80	46	5.7	0.697	0.628	0.0688	0.527
120	69	8.6	0.683	0.639	0.0438	0.648
200	115	14.4	0.679	0.651	0.0283	0.779
350	201	25.1	0.674	0.657	0.0177	0.852
600	345	43.1	0.672	0.662	0.0102	0.909
1200	689	86.2	0.670	0.665	0.0048	0.960
2400	1379	172.3	0.670	0.668	0.0019	0.994

Shrinkage: pull the coefficients a little towards zero

If the problem is that the coefficients are too extreme, the most direct fix is to pull them towards zero. That is shrinkage.

By how much? Van Houwelingen and Le Cessie gave a heuristic that needs nothing but the development data:

\hat{s} = \frac{\chi^2_{\mathrm{LR}} - p}{\chi^2_{\mathrm{LR}}}

The numerator is the model’s likelihood ratio chi-square minus the number of parameters; the denominator is the chi-square itself. Multiply every coefficient by this uniform shrinkage factor and the model stops being so extreme.

Why it is reasonable: part of the chi-square is real signal and part of it is what estimating p extra parameters buys you for free, and the expected value of that second part is exactly p. Subtract it and take the ratio, and you have “how much of this is signal”.

Can that number tell you in advance how far the slope will actually fall? The right-hand panel of the figure above is checking exactly that — the heuristic shrinkage factor, which only looks at the development data, against the calibration slope actually measured on the held-out patients. Across all seven sample sizes they track each other closely. At the smallest sample size the heuristic predicts 0.591 and the measurement comes out at 0.527.

ridge, lasso, elastic net: one idea, three penalties

Penalised regression writes the shrinkage directly into the estimation. Instead of maximising the log likelihood alone, you maximise

\ell(\beta) - \lambda \left[ \frac{1-\alpha}{2}\sum_j \beta_j^2 + \alpha \sum_j |\beta_j| \right]

The bracket is the penalty; $\lambda$ decides how hard it bites and $\alpha$ decides its shape:

	Penalty	Can a coefficient become exactly 0?	Typical use
ridge	Sum of squared coefficients (alpha is 0)	No, everything simply gets smaller	You want to keep every candidate, just stop them being extreme
lasso	Sum of absolute coefficients (alpha is 1)	Yes, so it does variable selection at the same time	Many candidates, and you want a shorter model
elastic net	A mixture of the two	Yes, but more gently than lasso	Predictors that are highly correlated (lasso will pick one of a correlated group more or less arbitrarily)

The difference comes from the shape of the penalty at zero: the absolute value has a corner there, so the optimum gets caught on zero; the square is smooth at zero, so the solution approaches zero without ever arriving.

Take 200 patients out of rotterdam (125 events, EPV 15.6) as the development data and keep the remaining 2782 as the held-out sample:

The left panel is the lasso coefficient path: the horizontal axis is the log of the penalty strength and eight lines show the eight predictors' coefficients multiplied by their own standard deviations, being pressed one after another to exactly zero as the penalty grows; a vertical dashed line marks the penalty chosen by cross-validation. The right panel is a histogram of the penalty values chosen when the same cross-validation is rerun sixty times on the same data; the distribution is clearly right-skewed and the largest value is about three point seven times the smallest. — Left: the lasso path. Coefficients are multiplied by their predictors' standard deviations before being drawn on one axis, since variables measured in years, millimetres and 0/1 are otherwise not comparable. The dashed vertical line is the penalty chosen by cross-validation. Right: the same data and the same code, changing only the random cross-validation folds, rerun sixty times.Plotting script `figures/scripts/B5-02-shrinkage.R`

Method	Non-zero coefficients	Held-out C-index	Held-out calibration slope
Maximum likelihood (no penalty)	8	0.6419	0.807
Uniform shrinkage factor	8	0.6419	0.962
ridge (alpha = 0)	8	0.6532	1.268
elastic net	5	0.6517	1.010
lasso (alpha = 1)	5	0.6501	0.953

Read this table across the rows, not down the C-index column.

The uniform shrinkage row has exactly the same C-index as maximum likelihood — because multiplying every coefficient by one constant changes nobody’s ranking, and the C-index only sees rankings. What it changes is the calibration slope: from 0.807 up to 0.962.

Shrinkage fixes calibration, not discrimination. Expecting a penalty to raise the AUC will usually end in disappointment.

The ridge row is a cautionary example in the other direction: its C-index is the highest in the table, but its calibration slope is 1.268 — above 1, which means it shrank too far. The predictions are now too conservative and the risk differences too compressed. More penalty is not better, and the $\lambda$ that cross-validation picks is chosen to minimise some prediction error, not to make the calibration slope equal 1.

A penalised coefficient is not an effect size

On this dataset lasso pressed 3 variables to zero. Here are the coefficients under four approaches:

Variable	Maximum likelihood	After uniform shrinkage	ridge	lasso
age	0.019	0.016	0.005	0.007
menopausal	-0.183	-0.153	0.022	0
tumour size	0.009	0.008	0.009	0.009
grade 3	0.639	0.536	0.359	0.517
nodes	0.071	0.059	0.048	0.067
log PgR	0.006	0.005	-0.015	0
log ER	-0.115	-0.097	-0.054	-0.087
hormonal tx	-0.062	-0.052	0.008	0

Not one of the non-zero numbers in this table can be written up as an “adjusted log hazard ratio”.

There are three layers to the reason:

Penalised estimation is deliberately biased. It trades bias for variance, which makes the overall prediction better, but the expected value of any single coefficient is no longer the true value.
There is no usable standard error. glmnet will not give you p-values or confidence intervals, and that is a design decision rather than a missing feature — the sampling distribution after penalisation has no simple form, and forcing a normal approximation onto it produces intervals that are too narrow.
A lasso “zero” does not mean “unrelated”. It means that at this $\lambda$ , in this sample, the variable’s contribution was too small to be worth the price of the penalty. Different people, different $\lambda$ , and it may well come back.

So the output of a penalised model has exactly one legitimate use: treat the whole set of coefficients as a machine for computing risk. To say anything about one variable’s effect, go back and fit an unpenalised model, and be ready to answer the whole confounding question.

The penalty parameter is itself an estimate

$\lambda$ is usually chosen by cross-validation: split into k folds, compute the cross-validated deviance for a whole grid of candidate $\lambda$ values, and take the smallest.

The problem is that the split itself is random. Same data, same code, changing only the random folds, rerun 60 times:

	Minimum	Median	Maximum
lambda chosen by cross-validation	0.0166	0.0264	0.0610
Number of non-zero coefficients kept	5	5	5

The largest value is 3.68 times the smallest.

There is an honest surprise here. The penalty strength moved by a factor of 3.68, but the number of variables kept did not change once across the 60 reruns.

The reason is that a few variables in this dataset (node count, grade, tumour size) are far stronger than the rest and survive the entire $\lambda$ range, while the weak ones were already gone — nothing sits on the boundary in between. On a dataset where strong and weak are not so cleanly separated, the same wobble would change what the model is made of.

The practical advice is to rerun cv.glmnet a few dozen times and take the median $\lambda$ , or to use lambda.1se instead (a heavier penalty than the minimum, and more stable), rather than trusting a single run. Whichever you choose, that decision is also made from the data, so it has to go inside the resampling loop — the reasoning is on the internal validation page.

What shrinkage cannot rescue

A penalty is not a cure for an inadequate sample. The table above already hinted at this: shrinkage pulled the calibration slope back, but doing so cannot conjure discrimination out of nothing. Keep the two apart. Uniform shrinkage multiplies every coefficient by the same factor, changes nobody’s relative ranking, and gives a held-out C-index identical to maximum likelihood (both 0.6419) — that is a mathematical necessity, not a measured difference. ridge shrinks different coefficients by different amounts, so it does move the ranking: on this single split the point estimate went from 0.6419 to 0.6532, a gap of 0.0113 and the largest in this table; but it is a point estimate from one split, and this page does not estimate its uncertainty.

What shrinkage does is admit “I know less than I appeared to”. It makes the model’s predictions more honest; it does not create signal the data never contained. When the sample is severely inadequate, the methodological literature’s advice is not to develop the model at all, rather than to prop it up with a penalty. The tool for deciding whether the sample is adequate is on the next page.

Run it yourself

library(survival); library(glmnet)
data(cancer, package = "survival")

rot <- rotterdam
rot$rfs_time  <- pmin(rot$rtime, rot$dtime)
rot$rfs_event <- as.integer(rot$recur == 1 | rot$death == 1)
rot$size_mm   <- c("<=20" = 15, "20-50" = 35, ">50" = 60)[as.character(rot$size)]
rot$grade3    <- as.integer(rot$grade >= 3)
rot$log_pgr   <- log1p(rot$pgr); rot$log_er <- log1p(rot$er)
cand <- c("age", "meno", "size_mm", "grade3", "nodes", "log_pgr", "log_er", "hormon")

set.seed(1)
i     <- sample(nrow(rot), 200)          # manufacture an inadequate-sample scenario
small <- rot[i, ]; rest <- rot[-i, ]

fit <- coxph(as.formula(paste("Surv(rfs_time, rfs_event) ~",
                              paste(cand, collapse = " + "))), data = small)

# Uniform shrinkage factor: a fitted model is all you need
lr <- 2 * diff(fit$loglik); p <- length(coef(fit))
(shrink <- (lr - p) / lr)
coef(fit) * shrink                       # the shrunken coefficients

# Penalised regression
x <- as.matrix(small[, cand]); y <- Surv(small$rfs_time, small$rfs_event)
cv_lasso <- cv.glmnet(x, y, family = "cox", alpha = 1, cox.ties = "breslow")
coef(cv_lasso, s = "lambda.min")         # the ones printed as . were dropped

# Stability of the penalty parameter: just rerun it
lam <- replicate(60, cv.glmnet(x, y, family = "cox", alpha = 1,
                               cox.ties = "breslow")$lambda.min)
range(lam); max(lam) / min(lam)

# Check on the held-out sample: read C-index and calibration slope together
b  <- as.numeric(coef(cv_lasso, s = "lambda.min"))
lp <- as.matrix(rest[, cand]) %*% b
concordance(Surv(rest$rfs_time, rest$rfs_event) ~ lp, reverse = TRUE)$concordance
coef(coxph(Surv(rfs_time, rfs_event) ~ lp, data = rest))

Verified with R 4.6.0 + survival 3.8.6 + glmnet 5.0. glmnet's Cox model handles ties with Breslow by default, unlike coxph's Efron; writing cox.ties out explicitly avoids drift across versions. glmnet standardises predictors internally, but the coefficients it returns are already back on the original scale.

import numpy as np
import pandas as pd
from lifelines import CoxPHFitter
from lifelines.utils import concordance_index

RD = "https://vincentarelbundock.github.io/Rdatasets/csv/"
d = pd.read_csv(RD + "survival/rotterdam.csv")

# The same derived columns as the R code above. rotterdam ships none of
# them, and without this every line below is a KeyError.
d["rfs_time"]  = d[["rtime", "dtime"]].min(axis=1)
d["rfs_event"] = ((d["recur"] == 1) | (d["death"] == 1)).astype(int)
d["size_mm"]   = d["size"].map({"<=20": 15, "20-50": 35, ">50": 60})
d["grade3"]    = (d["grade"] >= 3).astype(int)
d["log_pgr"]   = np.log1p(d["pgr"]); d["log_er"] = np.log1p(d["er"])
cand = ["age", "meno", "size_mm", "grade3", "nodes", "log_pgr", "log_er", "hormon"]
cols = cand + ["rfs_time", "rfs_event"]

rng = np.random.default_rng(1)
i = rng.choice(len(d), 200, replace=False)
small, rest = d.iloc[i], d.drop(d.index[i])

ml    = CoxPHFitter().fit(small[cols], "rfs_time", "rfs_event")
lasso = CoxPHFitter(penalizer=0.05, l1_ratio=1.0).fit(small[cols], "rfs_time", "rfs_event")
ridge = CoxPHFitter(penalizer=0.05, l1_ratio=0.0).fit(small[cols], "rfs_time", "rfs_event")

for name, m in [("ml", ml), ("lasso", lasso), ("ridge", ridge)]:
    lp = rest[cand].to_numpy() @ m.params_[cand].to_numpy()
    # mind the direction: concordance_index needs the risk score negated to be
    # in survival-time order
    print(name, concordance_index(rest["rfs_time"], -lp, rest["rfs_event"]))

scikit-learn has no penalised Cox model. lifelines' CoxPHFitter takes penalizer and l1_ratio, which correspond to glmnet's lambda and alpha, but the two packages define the lambda scale differently — do not move numbers across directly.

Common misuses

Misuse	Why it is wrong
Treating the variables lasso kept as “the important predictors”	They are the survivors of a trade-off at one penalty strength; another sample changes them
Reporting p-values or confidence intervals for penalised coefficients	Penalised estimates have no simple sampling distribution, and normal-approximation intervals are too narrow
Using penalised regression to fix an inadequate sample	It makes the predictions honest; it does not invent signal the data never held
Expecting shrinkage to raise the AUC	Uniform shrinkage does not change the ranking at all, so the C-index does not move a decimal place
Running `cv.glmnet` once and calling it settled	The folds are random, and rerunning gives a different penalty strength
Leaving the choice of penalty strength outside the resampling loop	That is also a decision made from the data, so the optimism comes straight back in
Assuming ridge is always safer than no penalty	The penalty cross-validation picks can be too heavy, pushing the calibration slope above 1
Comparing the size of penalised coefficients without standardising	The penalty treats all coefficients alike, so variables on large scales are punished disproportionately
Reading causation into a penalised model’s coefficients	The bias was introduced on purpose, for prediction; these are not effect estimates

Reproducing every number on this page

/opt/homebrew/bin/Rscript figures/scripts/B5-02-shrinkage.R

Read the figure

The answer comes from the same statistical output that produced this page's figures, not from a number typed in beside them.

The same eight-variable Cox model has an apparent C-index of 0.697 in a development sample of eighty people, and only 0.628 when applied to everyone left over. What does that gap say?

Show the answer and why

Correct answer: The gap is optimism of 0.069, and since the people left over come from the same cohort with the same measurements it cannot be a population difference

0.069 is apparent performance minus hold-out performance. The hold-out sample is deliberately drawn from the same cohort precisely to switch off the population-difference explanation, leaving overfitting as the only source of the gap. 0.002, the optimism in the largest row, shows the opposite of what the first option claims: same data, same model, only a larger sample, and overfitting all but disappears - so it is not a property of the data but a relationship between model complexity and sample size. Nor is the calibration slope of 0.527 an unrelated third thing: coefficients that are too extreme are exactly why the predictions are too extreme, and the discrimination gap and the falling slope are two faces of one mechanism.

The uniform shrinkage factor multiplies every coefficient by the same constant. Its hold-out C-index equals maximum likelihood exactly, while the calibration slope moves from 0.807 to 0.962. Why?

Show the answer and why

Correct answer: A C-index only sees ranking, and multiplying every coefficient by one positive number leaves the ranking untouched, so matching maximum likelihood at 0.6419 is arithmetic

0.6419 appearing twice is no accident: uniform shrinkage multiplies each coefficient by the same factor, so the ordering of any two linear predictors is unchanged, and the C-index only asks whether that ordering is right. Ridge at 0.6532 and lasso at 0.6501 do move, because they shrink different coefficients by different amounts - that is what changes the ordering, not shrinking more thoroughly. The division of labour is what to remember: shrinkage repairs calibration, moving the slope from 0.807 to 0.962, not discrimination. Expecting a penalty to raise the AUC usually ends in disappointment.

In the same hold-out sample, ridge has a calibration slope of 1.268. How should that number be read?

Show the answer and why

Correct answer: 1.268 is above 1, which means it has been shrunk too far: the predictions have become too conservative and the risk differences squeezed too close together, so heavier is not better

The target for a calibration slope is 1, and departures in either direction are problems. 0.962 is a few thousandths away from 1; 1.268 is nearly three tenths away, so calling the two distances comparable falls apart as soon as the numbers sit side by side. The directions differ too: below 1 the predictions are too extreme, above 1 they are too conservative and the high and low risks are pressed together. Ridge has not repaired more thoroughly, it has overshot; maximum likelihood at 0.807 and ridge at 1.268 are two opposite errors, not two positions along one road. And the penalty that cross-validation selects is chosen to minimise some prediction error, never to make the calibration slope equal 1.

Lasso pushes the coefficient for menopausal status to exactly zero, while the same variable has a maximum-likelihood coefficient of -0.183. Which statement is right?

Show the answer and why

Correct answer: -0.183 is the unpenalised estimate; the lasso zero only says that at this penalty, in this batch of people, the variable is not worth the penalty it costs

A lasso zero is the result of an optimisation, not the conclusion of a test. 0.022 pointing the opposite way from -0.183 is striking, but it is not a miscalculation: ridge distributes weight among correlated variables differently from maximum likelihood, and coefficient signs need not survive a penalty intact - which is one of the reasons penalised coefficients cannot be read as effect sizes. Uniform shrinkage at -0.153 is a different operation again: it multiplies the original estimate by a common factor and leaves every variable's relative weight alone, whereas lasso catches some variables on the kink in the penalty at zero. To say anything about whether this variable has an effect, go back to an unpenalised model and be ready to answer the whole confounding question.

The same data and the same code, rerun sixty times with only the cross-validation folds reshuffled, select penalties from 0.0166 to 0.0610, while the number of non-zero coefficients never changes at all. How should that be read?

Show the answer and why

Correct answer: The median of 0.0264 is the more trustworthy figure; the count holds still because the strong and weak variables here are far apart with nothing borderline in between

0.0264 is the median of sixty reruns, and taking it - or switching to the steadier lambda.1se - is the practical advice, rather than trusting the 0.0219 a single run happened to give. A number lying inside its own distribution is no evidence that it is trustworthy, especially when that range runs from 0.0166 to 0.0610. The unchanging count is not evidence of stability either: nodes, grade and tumour size are far stronger than the rest here, still uncrushed across the whole penalty range, while the weak ones were crushed long before it, so nothing sits on the edge waiting to flip. A swing of 3.6784 would change what is in the model in data whose variables were not so cleanly separated. And the larger point: choosing the penalty is itself a decision made from the data, so it has to go inside the resampling loop.

Van Houwelingen's uniform shrinkage factor needs only the development data. At the smallest sample size it gives 0.591, while the calibration slope actually measured in the hold-out sample is 0.527. What does that comparison mean?

Show the answer and why

Correct answer: 0.591 can be computed the moment the model is fitted, and it sits close to a calibration slope that is otherwise only measurable afterwards; its value is that it is cheap

0.591 and 0.527 sit close together, and so do their counterparts at all seven sample sizes, which is what the right-hand panel of this page's figure is checking. The row at 0.980 cannot be used to argue that the heuristic only works in large samples: the calibration slope in that row is close to 1 as well, so the two agree there too, and the row is really another way of writing down that overfitting is small when the sample is ample. Nor is 0.873 differing slightly from the slope in its row evidence of coincidence: the uniform shrinkage factor is a heuristic that assumes every coefficient should shrink by the same proportion, which is usually not true, so it delivers a good approximation rather than an identity. And the agreement is an observation on this dataset, not a guarantee.

Chapters that use this method

Prediction model study

Watch next

ENEhsan Karim· 30 minThe middle section covers overfitting and what it does to calibration, which is exactly the first two sections of this page.

ENNIHR Maudsley BRC· 62 minAn hour-long tour of the whole modelling process; it makes it easier to see where penalised regression sits in the pipeline.

Sources and licences

This page is original writing