1
$\begingroup$

I've read a few Q&As about this, but am still not sure I understand, why the coefficients from glmnet and caret models based on the same sample and the same hyper-parameters are slightly different. Would greatly appreciate an explanation!

I am using caret to train a ridge regression:

library(ISLR)
Hitters = na.omit(Hitters)
x = model.matrix(Salary ~ ., Hitters)[, -1] #Dropping the intercept column.
y = Hitters$Salary

set.seed(0)
train = sample(1:nrow(x), 7*nrow(x)/10)

library(caret)
set.seed(0)
train_control = trainControl(method = 'cv', number = 10)
grid = 10 ^ seq(5, -2, length = 100)
tune.grid = expand.grid(lambda = grid, alpha = 0)
ridge.caret = train(x[train, ], y[train],
                    method = 'glmnet',
                    trControl = train_control,
                    tuneGrid = tune.grid)
ridge.caret$bestTune
# alpha is 0 and best lambda is 242.0128

Now, I use the lambda (and alpha) found above to train a ridge regression for the whole data set. At the end, I extract the coefficients:

ridge_full <- train(x, y,
                    method = 'glmnet',
                    trControl = trainControl(method = 'none'), 
                    tuneGrid = expand.grid(
                      lambda = ridge.caret$bestTune$lambda, alpha = 0)
                    )
coef(ridge_full$finalModel, s = ridge.caret$bestTune$lambda)

Finally, using exactly the same alpha and lambda, I try to fit the same ridge regression using glmnet package - and extract coefficients:

library(glmnet)
ridge_full2 = glmnet(x, y, alpha = 0, lambda = ridge.caret$bestTune$lambda)
coef(ridge_full2)
$\endgroup$
  • 1
    $\begingroup$ Please elaborate on what you mean by "slightly different." After all, when you rerun either of these procedures, you ought to get different results, because the cross-validation is based on random subsets. $\endgroup$ – whuber Feb 7 '18 at 0:57
  • $\begingroup$ Slightly different means: I thought, when I do NOT use the cross-validation (code blocks 2 and 3) the regression is run using specific alpha and lambda and the whole data set. So, the coefficients should be identical in caret and glmnet. But they are not. $\endgroup$ – user3245256 Feb 7 '18 at 13:32
  • $\begingroup$ Please give a specific example of what you mean by "slightly different," for otherwise we haven't adequate information about the problem--or even whether there is one. There are many reasons why results may differ, ranging from differences in your input to floating point roundoff error to different algorithms to different tolerance values for convergence and on and on--many of which are benign. $\endgroup$ – whuber Feb 7 '18 at 14:00
  • $\begingroup$ If you run the R code I provided in the question, you'll see the coefficients from glmnet and from caret. They are different (usually 1st or 2nd decimal place). $\endgroup$ – user3245256 Feb 7 '18 at 17:14
2
$\begingroup$

It seems like a bug in caret's implementation.

First some notes about glmnet package:

  • The documentation of glmnet() recommends against giving one single value for lambda. It is preferable to warm start with large values values of lambda.
  • predict.glmnet() lets you override the value(s) of lambda which was used to train the model (cf s argument). However, when some supplied values of s differ from the original lambda, the default behaviour is to use linear interpolation rather than re-fit a model with lambda=s. This can be controlled with the exact argument of predict().

Caret provides a loop() function for its glmnet wrapper (see https://github.com/topepo/caret/blob/master/models/files/glmnet.R). This loop function will only fit one model per value of alpha, with the max corresponding value of lambda. The other lambda values are considered as "submodels". Their evaluation is deferred to predict() rather than fit(), using the glmnet::predict(s=...) feature mentioned above. This is efficient and fine but predict should specify exact=TRUE to obtain the same results as glmnet.

With this setting, I don't think that ridge.caret$finalModel was effectively trained with the optimal lambda.

$\endgroup$
  • $\begingroup$ Pierre, thank you - very interesting. May I ask: are you saying my code should have been different? It'd be great to see what a better version of that code would look like. All I am trying to do is to train a model in caret using grid search and then, once the optimal lambda is found, run a regression with the total data set to get its regression coefficients. $\endgroup$ – user3245256 Feb 14 '18 at 13:22
  • $\begingroup$ @user3245256 Your code was correct. You could even ridge.caret$finalModel which, I think, is supposed to be the same as your ridge_full$finalModel. You just faced the comoplex interactions between caret and glmnet... In the present case you may have directly used cross-validation facility in glmnet: final_model=cv.glmnet(...) and then predict(final_model, ...). $\endgroup$ – Pierre Gramme Feb 14 '18 at 22:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.