cv.glmnet vs glmnet

Question

I'm using glmnet to fit a ridge regression model on some data and evaluate the model's test MSE. The lambda value I select is derived from cross-validation. I'm using the College dataset from ISLR2, predicting applications to each college.

I have tried the following two approaches, and while in theory I should have the same result, I don't and I'm not sure why.

First way

• I use cv.glmnet() to perform cross validation on the data, extracting the lambda with the lowest validation MSE
• Then I fit a ridge regression model using glmnet() on the data with the previously computed lambda
• Predict response on the test data using the fitted model and compute the test MSE

# Perform the cross validation
cv.ridge <- cv.glmnet(model.matrix(Apps~.,train), train$Apps, alpha=0, nfold=100)
# Store the best lambda value
best.lambda <- cv.ridge$lambda.min

# Fit ridge regression model with that lambda
ridge.fit <- glmnet(model.matrix(Apps~.,train), train$Apps, alpha=0, lambda=best.lambda)

# Predict the test response 
pred.out <- predict(cv.ridge, newx = model.matrix(Apps~.,test), s=best.lambda)

# Compute test MSE
mean((pred.out- Y_test)^2)
>> 2206587

Second Way

• Once again I use cv.glmnet() to perform cross validation on the data
• I then predict the response on the test data directly with the cv.glmnet object, using the lambda value with the lowest validation MSE
• Compute the test MSE

# Perform the cross validation
cv.ridge <- cv.glmnet(model.matrix(Apps~.,train), train$Apps, alpha=0, nfold=100)
# Store the best lambda value
best.lambda <- cv.ridge$lambda.min

# Predict the test response
pred.out <- predict(cv.ridge, newx = model.matrix(Apps~.,test), s=best.lambda)

# Compute test MSE
mean((pred.out- Y_test)^2)
>>> 2204831

Why do these two approaches have different test MSE's? The only difference between the two ways is that in the second way I use the cv.glmnet object instead of the glmnet object in the predict() call.

I checked the coefficients of the models and they are not the same either.

# Coefficients from the glmnet() call on the specified lambda
coef(ridge.fit)

# Coefficients of the cv.glmnet() call given the same specified lambda
predict(cv.ridge, type='coefficients', s=best.lambda)

The coefficients are slightly different. Which I guess is why the test MSE's differ. But I'm not sure why this should be the case.

In both ways, since the lambda constraint specified is identical and the data used to fit the model is identical, shouldn't the resulting two ridge regression models be the same?

Follow up :

I have tried setting the s parameter in the predict call but that doesn't seem to work either.

cv.ridge MSE :

# Fit cross-validated ridge regression model
cv.ridge <- cv.glmnet(model.matrix(Apps~.,College), College$Apps, alpha=0, nfold=100)
# Make prediction using lambda that minimizes cross-val MSE
pred.out <- predict(cv.ridge, model.matrix(Apps~.,College), s="lambda.min")
# Compute MSE
mean((pred.out- College$Apps)^2)
>> 1358455

glmnet way :

# Explicitly fit a ridge regression model on the same data using the previously computed lambda that minimizes CV mse
ridge.fit <- glmnet(model.matrix(Apps~.,College), College$Apps, alpha=0, lambda=cv.ridge$lambda.min)
# Make prediction 
pred.out <- predict(ridge.fit, model.matrix(Apps~.,College))
# Compute MSE
mean((pred.out- College$Apps)^2)
>> 1359837

The issue persists.

Note that in the glmnet way I didn't set any s parameter since the provided glmnet object only contains the fitted model with the lambda which minimizes CV MSE. The same lambda value used when the cv.ridge object does the prediction as well.

Setting s='lambda.min' doesn't change the result either.

Answer 1

cv.glmnet returns an object of type cv.glmnet (using dots in class names is asking for trouble with method dispatching, but this does not seem to be the problem here), which has overwritten methods for predict and coef, namely predict.cv.glmnet and coef.cv.glmnet. When you read their docs with

> ?coef.cv.glmnet
> ?predict.cv.glmnet

you will notice that they have an option s for specifying the lambda used. The coefficients seen in your question refer to s="lambda.1se". Set it to s="lambda.min", and the result is identical with your model trained with lambda.min:

> cv.ridge <- cv.glmnet(model.matrix(Apps~.,College), College$Apps, alpha=0, nfold=100)
> coef(cv.ridge, s="lambda.min")
19 x 1 sparse Matrix of class "dgCMatrix"
                        1
(Intercept) -1.468326e+03
(Intercept)  .           
PrivateYes  -5.278781e+02
Accept       1.004588e+00
...
> pred.out <- predict(cv.ridge, model.matrix(Apps~.,College), s="lambda.min")
> mean((pred.out- College$Apps)^2)
[1] 1358455

Remark: What do your data sets train and test contain? Apparently something different from College because your MSE and coefficients are different.

Answer 2

The answer is, in part, about numerical accuracy. The thresh argument of glmnet controls the convergence threshold of the gradient descent algorithm used to fit the model. According to the documentation, the optimization stops when

the maximum change in the objective after any coefficient update is less than thresh times the null deviance.

For gaussian models the deviance is mean squared error. For this particular example, the outcome is number of applications received by US colleges in 1995, so the null deviance is rather high: it's on the same order as the mean squared errors reported in the question, ie., on the order of 1,000,000.

Let's vary the threshold to demonstrate how the difference between the two fitted models gets smaller as we require a more accurate solution. (Full R code listing is attached at the end.)

# null deviance on the original scale
mean((y - mean(y))^2)
#> [1] 14959182

# `abs_mse_diff` is the absolute difference in MSE between the cv.glmnet and glmnet models

abs_mse_diff(x, y, thresh = 1e-7)
#> [1] 1101.85
abs_mse_diff(x, y, thresh = 1e-8)
#> [1] 517.3258
abs_mse_diff(x, y, thresh = 1e-9)
#> [1] 113.6585
abs_mse_diff(x, y, thresh = 1e-10)
#> [1] 14.41409

The final cv.glmnet object is fitted to the full dataset, so there must be another source of randomness to obtain two numerically different solutions. However, they are "equivalent" within the precision specified by thresh. If you transform the response y to have unit scale, the difference in MSE between the cv.glmnet and the glmnet solutions is much smaller.

---

library("glmnet")
library("tidyverse")

mse <- function(actual, predicted) {
  mean((actual - predicted)^2)
}

abs_mse_diff <- function(x, y, ...) {
  cv.ridge <- cv.glmnet(
    x, y,
    alpha = 0,
    ...
  )
  ridge.fit <- glmnet(
    x, y,
    alpha = 0,
    lambda = cv.ridge$lambda.min,
    ...
  )

  abs(
    mse(y, predict(cv.ridge, x, s = "lambda.min"))
    - mse(y, predict(ridge.fit, x))
  )
}

data(College, package = "ISLR2")

College <- College %>%
  select(where(is.numeric)) %>%
  select(Apps, everything())

College <- as.matrix(College)

x <- College[, -1]
y <- College[, 1]

# null deviance on the original scale
mean((y - mean(y))^2)

abs_mse_diff(x, y, thresh = 1e-7)
abs_mse_diff(x, y, thresh = 1e-8)
abs_mse_diff(x, y, thresh = 1e-9)
abs_mse_diff(x, y, thresh = 1e-10)

y <- scale(y, center = FALSE, scale = TRUE)

# null deviance after scaling the response
mean((y - mean(y))^2)

abs_mse_diff(x, y, thresh = 1e-7)
abs_mse_diff(x, y, thresh = 1e-8)
abs_mse_diff(x, y, thresh = 1e-9)
abs_mse_diff(x, y, thresh = 1e-10)