R: caret (elasticnet): ridge regression: understanding the returned parameters

Question

I wanted to play around with the ridge regression in caret (which apparently uses elasticnet), so I did two experiments:

• use the original data
• use the modified data where the values of x2 are multiplied by 0.5.

The value of ridgeFit$finalModel$beta.pure in the first case is

x1 x2 x3 0 0.000000 0.000000 0.0000000 1 0.000000 0.000000 0.4803075 2 0.000000 3.245819 2.3878478 3 1.464703 2.543341 3.4790604

Where does that come from? (The tested lambdas were only three: $0$, $10^{-4}$ and $10^{-1}$).

Apparently, the last line corresponds to the computed parameters (see beta_true values in the code below). It this always the case?

Moreover, if I compare the coefficients for the variable x2 (the values b1 and b2), it turns out that b2 = 2 b1. This seems to be wrong since the optimisation function of ridge regression should be $$ \sum_{(x, y)} (y - \sum_i \beta_i x_i)^2 + \lambda ||\beta||_2^2\text{,} $$ so making $\beta_2$ twice as big in the second case should keep the preditions $\hat{y} =\sum_i \beta_i x_i$ unchanged, but should increase the penalty term, so choosing a somewhat smaller $\beta_2$ should be preferred (the chosen lambda was not 0).

The same happens (b2 = 2 b1) when

• the number of examples is larger, e.g., 100 or 1000
• I specify the possible values of lambda in tuneGrid parameter, e.g., tuneGrid = data.frame(lambda = 11.1)

The code:

library(caret)

A = matrix(runif(30), ncol=3)
beta_true = c(1.5, 2.5, 3.5)
Y = A %*% beta_true
Y = Y + runif(length(Y)) * 0.1

data = as.data.frame(A)
data$y = Y
colnames(data) = c("x1", "x2", "x3", "y")

set.seed(123)
ridgeFit = train(y ~ ., data=data, method="ridge")
print(ridgeFit)
print(ridgeFit$finalModel$beta.pure)
b1 = ridgeFit$finalModel$beta.pure[4,2]

data$x2 = 0.5 * data$x2
set.seed(123)
ridgeFit = train(y ~ ., data=data, method="ridge")
print(ridgeFit)
print(ridgeFit$finalModel$beta.pure)
b2 = ridgeFit$finalModel$beta.pure[4,2]
print(sprintf("b2 - 2 b1 = %f", b2 - 2 * b1))

---

EDIT:

Using glmnet's method glment directly seems to reflect the changes in the data correctly. However, that does not solve the original question.

Answer 1

If you specify method="ridge", the model is fitted using elastic net, you can check it enter link description here. To answer your questions:

i) beta.pure, the first line of all zeros is for the null model, only intercept.

ii) you need to check again what is the lambda in the 2nd model. Caret selects the model with the least deviance and puts it as the finalModel. If i reran your example, I get a chosen lambda of 0 or 1e-04 which is really small. If you increase the lambda you will see that your beta2 with the scaled down data gets penalized more. And in your example, you need to have a huge lambda before you see it shrinking..

So we can check this:

library(caret)
set.seed(111)
A = matrix(runif(30), ncol=3)
beta_true = c(1.5, 2.5, 3.5)
Y = A %*% beta_true
Y = Y + runif(length(Y)) * 0.1

data = as.data.frame(A)
data$y = Y
colnames(data) = c("x1", "x2", "x3", "y")
data2 = data
data2$x2 = 0.5 * data2$x2

To check these models, a good practice might be to write a function to keep track of these:

fit = function(dat){
ridgeFit = train(y ~ ., data=dat, method="ridge")
beta = ridgeFit$finalModel$beta.pure
data.frame(
b2 = beta[nrow(beta),2],
lambda = ridgeFit$finalModel$lambda
)
}

sapply(list(data=data,data2=data2),function(i)fit(i))

       data     data2   
b2     2.537871 5.076027
lambda 1e-04    0 

We can use run the ridge regression from enet, which will show you the effect of lambda, basically only at 100 does it start going down:

library(elasticnet)
fit_enet = function(dat,lam){
    beta = enet(x=as.matrix(dat[,1:3]),y=dat[,4],lambda=lam)$beta.pure
    data.frame(b2 = beta[nrow(beta),2],lambda = lam)
    }    

sapply(c(0.01,0.1,1,10,100),function(i){
fit_enet(data2,i)
})

       [,1]     [,2]     [,3]    [,4]     [,5]    
b2     5.048046 4.834376 4.05609 3.596255 3.532956
lambda 0.01     0.1      1       10       100