Why is cv.glmnet returning absurd coefficients when intercept term is omitted?

Question

x is a numeric matrix and y is a numeric vector:

x = structure(c(53, 36, 51, 51, 54, 35, 56, 60, 60, 60, 35, 59, 62, 
36, 38, 61, 64, 60, 92, 92, 62, 42, 65, 89, 62, 61, 62, 62, 62, 
35, 35, 37, 3.32, 3.1, 3.18, 3.39, 3.2, 3.03, 4.78, 4.72, 4.6, 
4.53, 2.9, 4.4, 4.31, 4.27, 4.41, 4.39, 7.32, 7.32, 7.45, 7.27, 
3.91, 3.75, 6.48, 6.7, 4.3, 4.02, 4.02, 3.98, 4.39, 2.75, 2.59, 
2.73, 3.42, 3.26, 3.18, 3.08, 3.41, 3.03, 4.57, 4.72, 4.41, 4.53, 
2.95, 4.36, 4.42, 3.94, 3.49, 4.39, 6.7, 7.2, 7.45, 7.26, 4.08, 
3.45, 5.8, 6.6, 4.3, 4.1, 3.89, 4.02, 4.53, 2.64, 2.59, 2.59), .Dim = c(32L, 
3L), .Dimnames = list(NULL, c("PT", "ITP", "PP")))

y = c(29, 24, 26, 22, 27, 21, 33, 34, 32, 34, 20, 36, 34, 23, 24, 
32, 40, 46, 55, 52, 29, 22, 31, 45, 37, 37, 33, 27, 34, 19, 16, 
22)

Without intercept

Code:

fit.ridge = glmnet(x, y, alpha = 0, intercept = FALSE)
plot(fit.ridge, xvar = "lambda", label = TRUE)
cv.ridge = cv.glmnet(x, y, alpha = 0, intercept = FALSE)
plot(cv.ridge)
coef(cv.ridge)

#4 x 1 sparse Matrix of class "dgCMatrix"
#                       1
#(Intercept) .           
#PT          7.877576e-36
#ITP         7.371832e-35
#PP          7.871337e-35

With intercept

Code:

fit.ridge = glmnet(x, y, alpha = 0, intercept = TRUE)
plot(fit.ridge, xvar = "lambda", label = TRUE)
cv.ridge = cv.glmnet(x, y, alpha = 0, intercept = TRUE)
plot(cv.ridge)
coef(cv.ridge)

#4 x 1 sparse Matrix of class "dgCMatrix"
#                   1
#(Intercept) 5.821492
#PT          0.194511
#ITP         1.420347
#PP          1.884496

Why do I get these absurd coefficients?

Answer 1

I can explain what you're seeing but not necessarily why it is the way it is. glmnet is starting the no-intercept solution at a much higher initial regularization penalty $\lambda_{max}$ than the with-intercept solution, and then hitting an early-stop in the path before it can explore better solutions.

How $\lambda_{max}$ is chosen

For $0 \lt \alpha \leq 1$, $\lambda_{max}$ is chosen as the highest value of $\lambda$ that still produces one nonzero coefficient (other than the intercept, which is not regularised). Ridge pushes coefficients asymptotically towards zero, whereas LASSO can entirely zero out coefficients. As $\alpha \rightarrow 0$, the LASSO contribution diminishes and the $\lambda$ value that results in exactly one nonzero coefficient gets higher and higher.

For $\alpha = 0$ or L2 regularization, coefficients are never regularized all the way to zero, even as $\lambda \rightarrow \infty$. How glmnet chooses $\lambda_{max}$ here is hard to glean from the source code or paper, but it seems like it sets $\alpha$ to a very small positive number and finds $\lambda_{max}$ the conventional way. $\alpha=0$ has a $\lambda_{max}$ of 120,761.2, the same value as for $\alpha=0.001$.

How the full $\lambda$ vector is chosen

Usually the minimum lambda $\lambda_{min}$ is chosen as $0.001 * \lambda_{max}$ and the algorithm searches 100 evenly-spaced-on-the-log-scale points between $\lambda_{min}$ and $\lambda_{max}$. For some reason, the no-intercept model stops after only 31 values. I have no idea why this is. glmnet will early stop searching the $\lambda$ vector early if the most recently fit $\lambda$ doesn't significantly improve training deviance, but this is not the case for your example. The source code is a complete black box. Who knows what's happening.

Why no-intercept has a higher $\lambda_{max}$

Since the intercept is regularization-free, it's a no-cost way for the model to fit the data. An intercept leaves the coefficients with less work to do, as the no-coefficient model has much lower error than it would without an intercept. $\lambda_{max}$ for the with-intercept model is much lower, at 8496.6. It searches as far as 0.8497 and finds good solutions along the way. The no-intercept model searches from 120761.2 to 7409.8, barely grazing the top of the with-intercept model's $\lambda$ path.

How to get a better set of solutions

If you just transplant the with-intercept $\lambda$ path into the no-intercept model, you get much better solutions.

cv.ridge.wi = cv.glmnet(x, y, alpha = 0, intercept = TRUE)
cv.ridge.ni = cv.glmnet(x, y, alpha = 0, intercept = FALSE, lambda = cv.ridge.wi$lambda)
plot(cv.ridge.ni)

Why coef(cv.ridge) returns very small numbers

coef on a cv.glmnet by default targets the s = "lambda.1se" heuristic, described in the docs. Since the no-intercept model searches 31 $\lambda$ and their error is very flat, $\lambda_{1se}$ is $\lambda_{max}$. You can see this from plot(cv.ridge). The with-intercept model's $\lambda_{1se}$ is much lower and the coefficients are more developed.

Will edit in some images later.

Answer 2

I'll begin by stating this is not the answer to this problem, however I was having a similar problem and have identified the cause in my case which may help someone looking here in the future if they happened to make the same mistake as me.

I was fitting the LASSO to a reasonably large training set (p = 20, n = 100000) and knew there was a relationship between the majority of features and the response, with some level of collinearity. The optimal fit selected by lambda.1se had all zero coefficients which seemed odd given my knowledge of the relationships in play.

The LASSO model was fit using cv.glmnet with alpha = 1.

Plotting the MSE using plot(fit_object) showed lambda.min as the smallest lamdbda and lambda.1se as the largest, with all other values in between the two, the error bars at each point were large.

• the model using lambda.min had all coefficients as zero apart from 3 and the intercept
• the model using lambda.1se had all coefficients as zero apart from the intercept

There were some massive outliers in my response variable due to having missed a step in my initial cleansing. This resulted in very little difference between MSE across all values of lambda and a large overall level of MSE. The outliers were a number of magnitudes higher and few in value.