Stat542: Variable Selection for Diabetes Data

Data

Brief Description of the Diabetes Data (Efron et al, 2004).

Ten baseline variables: age, sex, body mass index, average blood pressure, and six blood serum measurements were obtained for each of n = 442 diabetes patients, as well as the response of interest, a quantitative measure of disease progression one year after baseline.

data = read.table("https://liangfgithub.github.io/Data/diabetes.data", header=TRUE)
dim(data)
head(data)

Lasso

The purpose of this analysis is to demonstrate the following:

• Model size along the CV path plot is not monotonic;
• Lasso does not explore all possible sub-models, i.e., some variable sets are not listed on the path plot.

library(glmnet)
set.seed(1234)
X = data.matrix(data[, -11])
Y = data$Y

mylasso = glmnet(X, Y, alpha = 1)
plot(mylasso, xvar="lambda",label=TRUE)

Note 1

Check how the model sizes vary along the path.

mylasso$df  

When lambda (i.e., the regularization) decreases, the model size first increases from 0 to 10 and then decreases to 9 and then back to 10. So it is not true that the # of variables always increases when lambda decreases.

You can see how the coefficient of “S3” changes from non-zero to zero, and then back to non-zero, which corresponds to that white block in the right part of the matrix (see the image of the matrix).

lasso.coef = as.matrix(coef(mylasso))
# round(lasso.coef[, -c(1:60)], dig=2)
image(coef(mylasso))

lasso.varset is an 11-by-88 matrix denoting variable sets over 88 different values of lambda (denoted by s0, s1, …, s87).

lasso.coef = as.matrix(coef(mylasso))
lasso.varset = ifelse(lasso.coef==0, 0, 1)
lasso.varset

From the output above, we obtain the following summary of models on the path of Lasso:

• s0: intercept only (size = 0)
• s1: (BMI, S5) enters (size = 2)
• s8: BP enters (size = 3)
• s12: S3 enters (size = 4)
• s22: SEX enters (size = 5)
• s26: S6 enters (size = 6)
• s29: S1 enters (size = 7)
• s42: S4 enters (size = 8)
• s56: S2 enters (size = 9)
• s57: AGE enters (size = 10)
• s67: S3 leaves (size = 9)
• s73: S3 re-enters (size = 10)

S3 is an interesting variable; it seems to be related to other S variables. On Lasso path, S3 enters to size 4 model (4 non-intercept predictors), and stays for a long time and only leaves briefly leading to a size 9 model without S3.

Note 2

Next, you’ll see that the optimal model selected by AIC is not on the path plot.

library(leaps)
RSSleaps = regsubsets(X, Y, nbest=1, nvmax=10)
sumleaps = summary(RSSleaps, matrix=T)
sumleaps

Subset selection object
10 Variables  (and intercept)
    Forced in Forced out
AGE     FALSE      FALSE
SEX     FALSE      FALSE
BMI     FALSE      FALSE
BP      FALSE      FALSE
S1      FALSE      FALSE
S2      FALSE      FALSE
S3      FALSE      FALSE
S4      FALSE      FALSE
S5      FALSE      FALSE
S6      FALSE      FALSE
1 subsets of each size up to 10
Selection Algorithm: exhaustive
          AGE SEX BMI BP  S1  S2  S3  S4  S5  S6 
1  ( 1 )  " " " " "*" " " " " " " " " " " " " " "
2  ( 1 )  " " " " "*" " " " " " " " " " " "*" " "
3  ( 1 )  " " " " "*" "*" " " " " " " " " "*" " "
4  ( 1 )  " " " " "*" "*" "*" " " " " " " "*" " "
5  ( 1 )  " " "*" "*" "*" " " " " "*" " " "*" " "
6  ( 1 )  " " "*" "*" "*" "*" "*" " " " " "*" " "
7  ( 1 )  " " "*" "*" "*" "*" "*" " " "*" "*" " "
8  ( 1 )  " " "*" "*" "*" "*" "*" " " "*" "*" "*"
9  ( 1 )  " " "*" "*" "*" "*" "*" "*" "*" "*" "*"
10  ( 1 ) "*" "*" "*" "*" "*" "*" "*" "*" "*" "*"

If you look at optimal size (without intercept) 4, 6, 7, 8 models from regsubsets, none of them has S3. That means that if AIC/BIC ends up picking a model with size 4, 6, 7, or 8, then the AIC/BIC optimal model is not on the path of Lasso.

You can check which model returned by AIC/BIC. Note that msize0 denotes the model size including the intercept and msize = msize0 -1 denotes the model size without the intercept.

n = length(Y)
msize0 = apply(sumleaps$which, 1, sum);
msize = msize0 - 1
BIC = sumleaps$bic; BIC = BIC - min(BIC); BIC = BIC/max(BIC);
AIC = n*log(sumleaps$rss/n) + 2*msize0;
AIC = AIC - min(AIC); AIC = AIC/max(AIC);
cbind(msize, AIC, BIC)

   msize         AIC         BIC
1      1 1.000000000 1.000000000
2      2 0.322396838 0.230742097
3      3 0.199478409 0.123584530
4      4 0.132534927 0.083241551
5      5 0.029652507 0.000000000
6      6 0.000000000 0.004169454
7      7 0.005806613 0.050664745
8      8 0.013262202 0.099128337
9      9 0.027634213 0.155847729
10    10 0.043603930 0.214474226

plot(range(msize), c(0, 0.4), type="n", 
     xlab="Model Size (with Intercept)", ylab="Model Selection Criteria")
points(msize, AIC, col="blue", pch= 1)

points(msize, BIC, col="black", pch = 2)
legend("topright", lty=rep(1,3), col=c("blue", "black"), legend=c("AIC", "BIC"))