\[
S(z, \lambda) = \begin{cases} z - \lambda & z > \lambda \\
z + \lambda & z < -\lambda \\ 0 & |z| \leq \lambda \end{cases}
%= \textrm{sign}(z)(|z| - \lambda)_+
\]

Iterating over this is called coordinate descent and gives the solution

Packages

There are two main R implementations for finding lasso

{glmnet}: lasso = glmnet(X, Y, alpha=1).

Setting alpha = 0 gives ridge regression (as does lm.ridge in the MASS package)

Setting alpha\(\in (0,1)\) gives a method called the “elastic net” which combines ridge regression and lasso (regularization \(\alpha \|\beta\|_1 + (1-\alpha)\|\beta\|^2_2\)). More on that later.

Default alpha = 1 (it does lasso)

{lars}: lars = lars(X, Y)

lars() also does other things called “Least angle” and “forward stagewise” in addition to “forward stepwise” regression

The path returned by lars() is more useful than that returned by glmnet().

But you should use {glmnet}.

Choosing the \(\lambda\)

You have to choose \(\lambda\) in lasso or in ridge regression

lasso selects variables (by setting coefficients to zero), but the value of \(\lambda\) determines how many/which.

All of these packages come with CV built in.

However, the way to do it differs from package to package

{glmnet} version (same procedure for lasso or ridge)

lasso <-cv.glmnet(X, Y) # estimate full model and CV no good reason to call glmnet() itself# 2. Look at the CV curve. If the dashed lines are at the boundaries, redo and adjust lambdalambda_min <- lasso$lambda.min # the value, not the location (or use lasso$lambda.1se)coeffs <-coefficients(lasso, s ="lambda.min") # s can be string or a numberpreds <-predict(lasso, newx = X, s ="lambda.1se") # must supply `newx`

\(\widehat{R}_{CV}\) is an estimator of \(R_n\), it has bias and variance

Because we did CV, we actually have 10 \(\widehat{R}\) values, 1 per split.

Calculate the mean (that’s what we’ve been using), but what about SE?

par(mfrow =c(1, 2), mar =c(5, 3, 3, 0))plot(lasso) # a plot method for the cv fitplot(lasso$glmnet.fit) # the glmnet.fit == glmnet(X,Y)abline(v =colSums(abs(coef(lasso$glmnet.fit)[-1, drop(lasso$index)])), lty =2)

Paths with chosen lambda

ridge <-cv.glmnet(X, Y, alpha =0, lambda.min.ratio =1e-10) # added to get a minimumpar(mfrow =c(1, 4))plot(ridge, main ="Ridge")plot(lasso, main ="Lasso")plot(ridge$glmnet.fit, main ="Ridge")abline(v =sum(abs(coef(ridge)))) # defaults to `lambda.1se`plot(lasso$glmnet.fit, main ="Lasso")abline(v =sum(abs(coef(lasso)))) # again, `lambda.1se` unless told otherwise

Degrees of freedom

Lasso is not a linear smoother. There is no matrix \(S\) such that \(\widehat{\y} = \mathbf{S}\y\) for the predicted values from lasso.

We can’t use cv_nice().

We don’t have \(\tr{\mathbf{S}} = \textrm{df}\) because there is no \(\mathbf{S}\).

However,

One can show that \(\textrm{df}_\lambda = \E[\#(\widehat{\beta}_\lambda \neq 0)] = \E[||\widehat{\beta}_\lambda||_0]\)

The proof is PhD-level material

Note that the \(\widehat{\textrm{df}}_\lambda\) is shown on the CV plot and that lasso.glmnet$glmnet.fit$df contains this value for all \(\lambda\).

Other flavours

The elastic net

generally used for correlated variables that combines a ridge/lasso penalty. Use glmnet(..., alpha = a) (0 < a < 1).

Grouped lasso

where variables are included or excluded in groups. Required for factors (1-hot encoding)

Relaxed lasso

Takes the estimated model from lasso and fits the full least squares solution on the selected covariates (less bias, more variance). Use glmnet(..., relax = TRUE).

Dantzig selector

a slightly modified version of the lasso

Lasso cinematic universe

SCAD

a non-convex version of lasso that adds a more severe variable selection penalty

\(\sqrt{\textrm{lasso}}\)

claims to be tuning parameter free (but isn’t). Uses \(\Vert\cdot\Vert_2\) instead of \(\Vert\cdot\Vert_1\) for the loss.

Generalized lasso

Adds various additional matrices to the penalty term (e.g. \(\Vert D\beta\Vert_1\)).

Arbitrary combinations

combine the above penalties in your favourite combinations

Warnings on regularized regression

This isn’t a method unless you say how to choose \(\lambda\).

The intercept is never penalized. Adds an extra degree-of-freedom.

Predictor scaling is very important.

Discrete predictors need groupings.

Centering the predictors is important

(These all work with other likelihoods.)

Software handles most of these automatically, but not always. (No Lasso with factor predictors.)

Next time…

What happens when we’re tired of all this linearity.