Lasso Regression Example using glmnet package in R

More details please refer to the link below: (https://web.stanford.edu/~hastie/glmnet/glmnet_alpha.html#lin)

This post shows how to use glmnet package to fit lasso regression and how to visualize the output. The description of data is shown in here.

dt <- readRDS(url("https://s3.amazonaws.com/pbreheny-data-sets/whoari.rds"))
attach(dt)
fit <- glmnet(X, y)

1 Visualize the coefficients

plot(fit)

1.1 Label the path

plot(fit, label = TRUE)

The summary table below shows from left to right the number of nonzero coefficients (DF), the percent (of null) deviance explained (%dev) and the value of \lambda (Lambda).

We can get the actual coefficients at a specific \lambda whin the range of sequence:

coeffs <- coef(fit, s = 0.1) 
coeffs.dt <- data.frame(name = coeffs@Dimnames[[1]][coeffs@i + 1], coefficient = coeffs@x) 

# reorder the variables in term of coefficients
coeffs.dt[order(coeffs.dt$coefficient, decreasing = T),]

          name   coefficient
1  (Intercept)  3.010501e+00
7          hfa  4.634502e-01
12         aro  2.139111e-01
9          inc  1.298769e-01
6       convul  9.284738e-02
16         whz  9.018865e-02
14         afe  8.749033e-02
10         lcw  7.538421e-02
15        absu  5.758727e-02
13         qcr  2.952464e-02
11         str  2.298848e-02
3         temp  9.842032e-03
5           rr  7.309241e-03
8          hfe  4.870858e-03
2         wght -3.138053e-05
4          age -8.724823e-03

Also, it can allow people to make predictions at specific \lambda with new input data:

nx = matrix(rnorm(nrow(dt$X)*ncol(dt$X)), nrow = nrow(dt$X), ncol = ncol(dt$X))
pred <- predict(fit, newx = nx, s = c(0.1, 0.05)) 
head(pred, 20)

            s1         s2
 [1,] 3.465773  1.5135872
 [2,] 3.088525  1.3815177
 [3,] 3.587867  1.4431828
 [4,] 3.771435  1.5411323
 [5,] 3.595251  1.8061536
 [6,] 3.419286  1.3059683
 [7,] 3.195521  1.4421734
 [8,] 1.612829 -1.5908158
 [9,] 2.536198  0.1630636
[10,] 3.244446  0.5152836
[11,] 2.967109  1.5045971
[12,] 2.723644  1.0691778
[13,] 3.313086  0.9604486
[14,] 3.126201  0.9068035
[15,] 2.826930  0.1012069
[16,] 3.241242  0.2246355
[17,] 3.752990  2.2036950
[18,] 2.811372  1.0081468
[19,] 4.364803  3.1690336
[20,] 4.497694  3.1981861

cv.glmnet is the function to do cross-validation here.

X <- dt$X
y <- dt$y
cv.fit <- cv.glmnet(X, y)

Plotting the object gives the selected \lambda and corresponding Mean-Square Error.

plot(cv.fit)

We can view the selected \lambda’s and the corresponding coefficients, For example,

cv.fit$lambda.min

[1] 0.02081079

cv.fit$lambda.1se

[1] 0.1011946

lambda.min returns the value of \lambda that gives minimum mean cross-validated error. The other \lambda saved is lambda.lse, which gives the most regularized model such that error is within one standard error of the minimum. To use that, we only need to replace lambda.min with lambda.lse above.

# create a function to transform coefficient of glmnet and cvglmnet to data.frame
coeff2dt <- function(fitobject, s) {
  coeffs <- coef(fitobject, s) 
  coeffs.dt <- data.frame(name = coeffs@Dimnames[[1]][coeffs@i + 1], coefficient = coeffs@x) 

  # reorder the variables in term of coefficients
  return(coeffs.dt[order(coeffs.dt$coefficient, decreasing = T),])
}

coeff2dt(fitobject = cv.fit, s = "lambda.min") %>% head(20)

     name coefficient
22    str  0.67960398
37    whz  0.58664385
12    hfa  0.53643828
11 convul  0.39597719
28    aro  0.33919548
39 puskin  0.28612411
14    fde  0.23064257
41    oab  0.20502318
30    afe  0.16702394
18    inc  0.14657295
13    hfe  0.12027053
5    temp  0.08881899
20    lcw  0.08607789
38   smi2  0.06967305
29    qcr  0.05343341
9    omph  0.03806023
19    sr1  0.03314876
32    crs  0.03202430
23    gru  0.02295646
36    hyp  0.01462041

coeffs.table <- coeff2dt(fitobject = cv.fit, s = "lambda.min")
ggplot(data = coeffs.table) +
  geom_col(aes(x = name, y = coefficient, fill = {coefficient > 0})) +
  xlab(label = "") +
  ggtitle(expression(paste("Lasso Coefficients with ", lambda, " = 0.0275"))) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "none") 

2 Elastic net

As an example, we can set \alpha=0.2

fit2 <- glmnet(X, y, alpha = 0.2, weights = c(rep(1, 716), rep(2, 100)), nlambda = 20)

print(fit2, digits = 3)

Call:  glmnet(x = X, y = y, weights = c(rep(1, 716), rep(2, 100)), alpha = 0.2,      nlambda = 20) 

   Df  %Dev Lambda
1   0  0.00 3.2500
2   7 14.02 2.0000
3  12 26.98 1.2300
4  13 33.81 0.7590
5  18 38.03 0.4670
6  23 41.35 0.2880
7  29 43.24 0.1770
8  39 45.05 0.1090
9  47 46.17 0.0672
10 52 46.82 0.0414
11 57 47.15 0.0255
12 58 47.29 0.0157
13 60 47.35 0.0097
14 60 47.38 0.0060
15 64 47.39 0.0037
16 65 47.39 0.0023
17 66 47.40 0.0014
18 66 47.40 0.0009
19 66 47.40 0.0005

According to the default internal settings, the computations stop if either the fractional change in deviance down the path is less than 10^{-5} or the fraction of explained deviance reaches 0.999.

plot(fit2, xvar = "lambda", label = TRUE)

# plot against %deviance
plot(fit2, xvar = "dev", label = TRUE)

predict(fit2, newx = X[1:5, ], type = "response", s = 0.03)

           s1
[1,] 3.170055
[2,] 4.869413
[3,] 4.255960
[4,] 5.018313
[5,] 2.945014