1.2.7 Choosing λ by cross-validation
A regularization path offers a fit at every λ. Cross-validation chooses among them by estimating how well each fit predicts observations it was not fitted to. Every family has a cross-validation procedure, the equivalent of R’s cv.glmnet: elnet-cv, logistic-cv, multinomial-cv, cox-cv, poisson-cv and mgaussian-cv. Each takes the arguments of the family’s path fitter and works as R’s does:
Fit the path to all the data. Its λ values are the candidates.
Split the observations at random into K folds, 10 by default. For each fold, fit the path to the other K − 1 folds, its training data.
Predict each observation from the path whose training data left it out, at every candidate λ, and measure the error of each prediction: by default the squared error for the Gaussian families and the deviance for the others.
Average the errors within each fold. At each λ, the mean of the K fold averages, weighted by the folds’ sizes, is the cross-validated error, and their spread gives its standard error.
Choose λ: lambda-min is the λ with the smallest cross-validated error, and lambda-1se the largest λ whose error is within one standard error of that smallest error.
Unless #:lambda fixes the sequence, each fold’s path chooses its own λ values, as in R, and is evaluated at the candidates by the interpolation of Predicting along a path.
Here are 60 observations of 8 predictors, of which only the first three carry signal, with y = 1 + 3x₁ − 2x₂ + x₃ plus noise:
> (random-seed 8) > (define (noise) (- (* 2 (random)) 1))
> (define X60 (for/list ([i (in-range 60)]) (for/list ([j (in-range 8)]) (noise)))) > (define beta '(3.0 -2.0 1.0 0.0 0.0 0.0 0.0 0.0))
> (define y60 (for/list ([row (in-list X60)]) (+ 1.0 (for/sum ([b (in-list beta)] [x (in-list row)]) (* b x)) (noise)))) > (define cv (elnet-cv X60 y60)) > cv
#<glmnet-cv:gaussian Mean-Squared Error
Lambda Index Measure SE Nonzero
min 0.05424 37 0.3984 0.05967 6
1se 0.1656 25 0.4468 0.06523 3>
The result is a glmnet-cv. It prints like R’s print.cv.glmnet, with the index counted from 0: the measure, then a row for each of the two choices, with its λ, its index among the candidates, the cross-validated error, its standard error and the number of nonzero coefficients (see Printing).
1.2.7.1 Reading lambda-min and lambda-1se
The smallest cross-validated error is only an estimate, and the standard error says how uncertain it is. lambda-min minimizes the estimate. lambda-1se takes the simplest model, the one with the most regularization, that the estimate cannot tell apart from the best, and it is what R uses by default. Here it keeps exactly the three predictors that matter, while lambda-min also keeps three of the noise predictors, with small coefficients:
> (glmnet-cv-lambda-min cv) 0.054239392702548954
> (glmnet-cv-lambda-1se cv) 0.16563926820045807
> (define (round3 x) (/ (round (* 1000 x)) 1000))
> (define (rounded v) (for/list ([b (in-vector v)]) (round3 b))) > (rounded (coef cv #:lambda 'lambda-min)) '(0.915 2.904 -1.705 0.869 -0.009 -0.111 0.0 0.0 -0.034)
> (rounded (coef cv #:lambda 'lambda-1se)) '(0.878 2.718 -1.471 0.696 0.0 0.0 0.0 0.0 0.0)
The whole curve is in the result: glmnet-cv-lambda holds the candidates, and glmnet-cv-cvm, glmnet-cv-cvsd and glmnet-cv-nzero the error, its standard error and the number of nonzero coefficients at each. Every fifth candidate, from the largest:
> (for ([lam (in-vector (glmnet-cv-lambda cv))] [err (in-vector (glmnet-cv-cvm cv))] [se (in-vector (glmnet-cv-cvsd cv))] [nz (in-vector (glmnet-cv-nzero cv))] [i (in-naturals)] #:when (zero? (remainder i 5))) (printf "~a: λ = ~a, error = ~a ± ~a, nonzero = ~a\n" i (round3 lam) (round3 err) (round3 se) nz))
0: λ = 1.695, error = 4.42 ± 0.464, nonzero = 0
5: λ = 1.065, error = 2.723 ± 0.303, nonzero = 1
10: λ = 0.669, error = 1.721 ± 0.205, nonzero = 2
15: λ = 0.42, error = 0.917 ± 0.122, nonzero = 3
20: λ = 0.264, error = 0.579 ± 0.082, nonzero = 3
25: λ = 0.166, error = 0.447 ± 0.065, nonzero = 3
30: λ = 0.104, error = 0.408 ± 0.06, nonzero = 4
35: λ = 0.065, error = 0.4 ± 0.06, nonzero = 5
40: λ = 0.041, error = 0.402 ± 0.06, nonzero = 6
45: λ = 0.026, error = 0.409 ± 0.062, nonzero = 6
50: λ = 0.016, error = 0.416 ± 0.065, nonzero = 8
55: λ = 0.01, error = 0.421 ± 0.067, nonzero = 8
60: λ = 0.006, error = 0.424 ± 0.068, nonzero = 8
The error falls quickly as the three real predictors enter, reaches its minimum, and then rises slowly as the noise predictors are fitted. plot-cv, from glmnet/plot, draws this curve as R’s plot draws a cv.glmnet result; see Cross-validation curves.
1.2.7.2 Predicting at the chosen λ
A glmnet-cv is a glmnet-model? whose path is the full-data path, so predict and coef work on it. As R’s predict.cv.glmnet does, they default to lambda-1se; #:lambda also takes 'lambda-min, 'lambda-1se or any number. deviance-ratio answers at lambda-1se:
> (define new-rows '((0.5 -0.5 0.0 0.0 0.0 0.0 0.0 0.0) (0.0 0.0 1.0 0.9 -0.9 0.9 -0.9 0.9))) > (predict cv new-rows) '(2.9728620818776497 1.5739136578358095)
> (predict cv new-rows #:lambda 'lambda-min) '(3.2201276679622577 1.84557015752822)
> (equal? (coef cv) (coef (glmnet-cv-path cv) #:lambda (glmnet-cv-lambda-1se cv))) #t
> (deviance-ratio cv) 0.9109831280347133
The true means of the two rows are 1 + 1.5 + 1 = 3.5 and 1 + 1 = 2. Both predictions fall short of them, those at lambda-1se, with the stronger penalty, by more.
1.2.7.3 Folds
#:nfolds sets the number of folds, at least 3. The folds are drawn by random-fold-ids from current-pseudo-random-generator, so a second call gives different folds and a slightly different curve. To repeat a result, seed the generator, or pass the folds with #:fold-ids: one fold id per observation, counting from 0. A result records its folds, so they can be reused:
> (glmnet-cv-lambda-min (elnet-cv X60 y60)) 0.03406391938001054
> (define (seeded seed) (parameterize ([current-pseudo-random-generator (make-pseudo-random-generator)]) (random-seed seed) (elnet-cv X60 y60 #:nfolds 5))) > (equal? (seeded 1) (seeded 1)) #t
> (define folds (glmnet-cv-fold-ids cv)) > (equal? (elnet-cv X60 y60 #:fold-ids folds) cv) #t
Fixed folds are also how to compare models fairly: cross-validate each value of #:alpha on the same folds, and the differences in error come from the models rather than from the split. Choosing α and λ does this.
1.2.7.4 Measures of error
#:type-measure chooses how a prediction’s error is measured, with R’s names and R’s default for each family:
Measure |
| Error of a prediction |
| Families |
'mse |
| squared error on the response scale |
| all but Cox; the default for the Gaussian families |
'deviance |
| the family's deviance |
| all; the default for the others |
'mae |
| absolute error on the response scale |
| all but Cox |
'class |
| misclassification rate |
| binomial, multinomial |
'auc |
| area under the ROC curve, per fold |
| binomial |
'C |
| Harrell's concordance index, per fold |
| Cox |
For the Gaussian families 'deviance is the squared error, as in R. For classifiers, 'mse and 'mae compare the predicted probabilities with 0/1 indicators of the classes. 'auc and 'C are better when larger, so lambda-min maximizes them.
> (define labels (for/list ([v (in-list y60)]) (if (> v 1.0) 1 0))) > (logistic-cv X60 labels #:type-measure 'class #:fold-ids folds)
#<glmnet-cv:binomial Misclassification Error
Lambda Index Measure SE Nonzero
min 0.07073 17 0.05 0.02546 3
1se 0.07763 16 0.06667 0.02722 3>
By default the errors are averaged within each fold first, which R calls grouped; #:grouped? #f computes the error and its standard error over the individual observations instead. As in R, the folds are never grouped when they average fewer than 3 observations, 'auc needs an average of 10 observations per fold and otherwise falls back to 'deviance, and each of these changes logs a warning.
For Cox models, the deviance of a fold is computed as R computes it. Grouped, it is the deviance of all the data less that of the fold’s training data, at the fold’s coefficients. Ungrouped, it is the deviance of the held-out fold alone, which is undefined when the fold has no event or when its first event is among its last two observations in time order; R stops then, and so does cox-cv, naming the fold. When the folds average fewer than 10 observations, the Cox deviance is grouped even with #:grouped? #f, with a warning, as in R.