# What does the varImp function in the caret package actually compute for a glmnet (elastic net) object

I am fitting an elastic net model with glmnet via the caret package with 189 predictors and a binomial criteria (a,b)

lassocont <- trainControl(method='repeatedcv',
repeats=10,
returnResamp='final',
allowParallel = TRUE,
seeds = theseed10,
classProbs = TRUE,
summaryFunction = twoClassSummary
)

lasso <- train(x,
y,
method='glmnet',
metric = "Spec",
preProc = c("center", "scale"),
family="binomial",
tuneLength = 60,
#tuneGrid = lassotune,
trControl = lassocont)


The final model uses an alpha = 0.1 and a lambda = 0.1.

Consequently I print the confusion matrix

Confusion Matrix and Statistics

Reference
Prediction a     b
a  28     7
b   1    13

Accuracy : 0.8367
95% CI : (0.7034, 0.9268)
No Information Rate : 0.5918
P-Value [Acc > NIR] : 0.0002156

Kappa : 0.6456
Mcnemar's Test P-Value : 0.0770999

Sensitivity : 0.9655
Specificity : 0.6500
Pos Pred Value : 0.8000
Neg Pred Value : 0.9286
Prevalence : 0.5918
Detection Rate : 0.5714
Detection Prevalence : 0.7143
Balanced Accuracy : 0.8078

'Positive' Class : a


However, I would like to know which variables are most contributory to the model as well as which predictors do deviate from zero in the equation. Therefore I do request variable importance via

varImp(lasso,scale=F)


Now I can see which variables are most helpful to predict the positive class, which are zero and those that do not predict the positive class.

To cut a long story short:

-What do these variable importance measures actually mean?

-How are the calculated for glmnet objects?

-How can I interpret the number, let's say in an article?

-Are there any pitfalls and limitations associated with this measure?

-Does the glmnet varImp measure take correlation structures into account?

Your help is very much appreciated!

thx Clemens

• This is not related to your question, but I don't think you need to pass preProc = c("center", "scale") in your caret::train() function, because glmnet itself has a standardize option which defaults to TRUE. Description for standardize argument: Logical flag for x variable standardization, prior to fitting the model sequence. The coefficients are always returned on the original scale. Default is standardize=TRUE. – Boxuan Sep 22 '16 at 13:32
• True - Regularization without standardization makes no sense. – Clemens Oct 12 '16 at 13:37

For these models, they are regression the coefficients for the final Model. Big coefficients are associated with larger effects. Using scale = FALSE is good here so you can also get the signs too.

There are always pitfalls with these measures depending on how you want to measure importance. They don't measure lack of fit at all, so if your model is 51% accurate, they are not very reflective of the data. In the case of regression coefficients, main effects are misleading when interactions are present and so on.

As for correlation between predictors, Friedman et al. (2010, JSS) state:

Ridge regression is known to shrink the coefficients of correlated predictors towards each other, allowing them to borrow strength from each other. In the extreme case of $k$ identical predictors, they each get identical coefficients with $1/k^{th}$ the size that any single one would get if fit alone.[...]

Lasso, on the other hand, is somewhat indifferent to very correlated predictors, and will tend to pick one and ignore the rest.

We have a pretty good example of that in Section 6.4 of APM

Max

• I don't think this comment answers how VarImp calculates the question "How are the calculated for glmnet objects?" Could you please elaborate. – andor kesselman Mar 14 '16 at 17:27

caret actually looks at the final coefficients of the fit and then takes the absolute value to rank the coefficients. Then the ranked coefficients are stored as variable importance.

To view the source code, you can type

getModelInfo("glmnet")$glmnet$varImp


To summarize, these are the lines to calculate it:

function(object, lambda = NULL, ...) {

## skipping a few lines

beta <- predict(object, s = lambda, type = "coef")
if(is.list(beta)) {
out <- do.call("cbind", lapply(beta, function(x) x[,1]))
out <- as.data.frame(out)
} else out <- data.frame(Overall = beta[,1])
out <- abs(out[rownames(out) != "(Intercept)",,drop = FALSE])
out
}