Add error for ridge regression with glmnet

[juliasilge]

Based on simulations, it looks like pure ridge regression (mixture = 0) is not correct using the glmnet engine. Let's add an error if a user tries to use mixture = 0 with glmnet.

[topepo]

The defaults for glmnet::glmnet() give poor (if not wrong results) when a pure ridge regression model is used.

To show this, let's use a simulated data set with known coefficients:

library(AppliedPredictiveModeling)
library(tidymodels)

# lp = 0 - 4 * x1 + 4 * x2
# cor[x1, x2] about 0.62
set.seed(1)
dat <- easyBoundaryFunc(1000, interaction = 0) %>% select(-prob)
x <- as.matrix(dat[, 1:2])
y <- dat$class

An unpenalized model gives these coefficients:

glm_mod <- glm(class ~ X1 + X2, data = dat, family = binomial)
coefs <- coef(glm_mod)
coefs

## (Intercept)          X1          X2 
##  0.04210609  3.95833315 -3.99492631

We'll show the coefficients for the glmnet model:

lr_coefs <- function(pen, alpha, ...) {
  mod <- glmnet::glmnet(x, y, family = "binomial", alpha = alpha, ...)
  betas <- coef(mod, s = pen)
  betas <- as.matrix(betas)
  betas <- as.data.frame(betas)
  term <- rownames(betas)
  betas <- stack(betas)
  betas$lambda <- rep(pen, each = length(term))
  betas$alpha <- alpha
  betas$label <- paste("alpha =", signif(alpha))
  betas$term <- rep(term, length(pen))
  betas$ind <- NULL
  betas <- betas[betas$term != "(Intercept)",]
  betas
}

penalties <- 10^seq(-4, 1, length.out = 50)

ridge_coefs <- lr_coefs(penalties, alpha = 0)
lasso_coefs <- lr_coefs(penalties, alpha = 1)

glmn_res <- bind_rows(lasso_coefs, ridge_coefs) 

glmn_res %>% 
  ggplot(aes(x = lambda, y = values, col = term)) + 
  geom_line() + 
  scale_x_log10() +
  geom_hline(yintercept = coefs[-1], lty = 3) +
  facet_wrap(~ label) +
  labs(y = "coefficients")

The ridge regression coefficients not appear to converge the unpenalized model estimates.

[topepo]

One possible fix is to use the lambda.min.ratio argument. ?glmnet has:

Smallest value for lambda, as a fraction of lambda.max, the (data derived) entry value (i.e. the smallest value for which all coefficients are zero). The default depends on the sample size nobs relative to the number of variables nvars. If nobs > nvars, the default is 0.0001, close to zero. If nobs < nvars, the default is 0.01. A very small value of lambda.min.ratio will lead to a saturated fit in the nobs < nvars case. This is undefined for "binomial" and "multinomial" models, and glmnet will exit gracefully when the percentage deviance explained is almost 1.

It turns out that using lambda.min.ratio = 0 will fix this issue:

We could use the translate() function to add this argument if lambda.min.ratio was specified. The potential problem might occur when nobs < nvars. translate() cannot determine this since the data are not exposed to that function.

Alternatively, we could document this in the model functions that use this engine and hope that people read the comments:

[bnaras]

@topepo to match the $\lambda$ sequence, specify it in the call to glmnet thus in lr_coefs:

lr_coefs <- function(pen, alpha, ...) {
  mod <- glmnet::glmnet(x, y, family = binomial, alpha = alpha, lambda = pen, ...)
  betas <- coef(mod, s = pen)
  betas <- as.matrix(betas)
  betas <- as.data.frame(betas)
  term <- rownames(betas)
  betas <- stack(betas)
  betas$lambda <- rep(pen, each = length(term))
  betas$alpha <- alpha
  betas$label <- paste("alpha =", signif(alpha))
  betas$term <- rep(term, length(pen))
  betas$ind <- NULL
  betas <- betas[betas$term != "(Intercept)",]
  betas
}

Furthermore, since you also want the coefs for the unpenalized model in your sequence, add zero to the $\lambda$ sequence:

penalties <- c(0, 10^seq(-4, 1, length.out = 50))

This will, of course, emit a warning as the documentation specified but will yield this figure below.

[github-actions]

This issue has been automatically locked. If you believe you have found a related problem, please file a new issue (with a reprex: https://reprex.tidyverse.org) and link to this issue.