2
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ glmnet normalises the input by default. have you checked the same is not done with caret? $\endgroup$ – seanv507 Mar 2 '20 at 17:42
  • $\begingroup$ @seanv507 The same IS done with caret: additional arguments can be passed to the train method; I tried and passed standardize = FALSE to the train calls above and nothing changed. $\endgroup$ – Antoine Mar 2 '20 at 18:06
1
$\begingroup$

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      
$\endgroup$
7
  • $\begingroup$ I am so sorry, but the question is not "Which $\lambda$ is the optimal one" but rather i) what is inside ridgeFit$finalModel$beta.pure, and ii) why b2 = 2 b1 (independently of $\lambda$). $\endgroup$ – Antoine Apr 9 '20 at 12:00
  • 1
    $\begingroup$ beta are the coefficients.. beta.pure, the first row is for a model with only intercept, it is used in calculating the null deviance $\endgroup$ – StupidWolf Apr 9 '20 at 12:02
  • $\begingroup$ Ok, so every next row adds one additional term to the model, I suppose. i) DONE (Thanks!) ii) remains open $\endgroup$ – Antoine Apr 9 '20 at 12:04
  • $\begingroup$ i don't see the b2 = 2 b1 in your example, basically, your assumptions only hold when the coefficients are penalized according to the lambda $\endgroup$ – StupidWolf Apr 9 '20 at 12:06
  • $\begingroup$ however if you use caret or cv.glmnet, it automatically goes to the lower lambda values, which minimize deviance, and throws you the coefficients which are close to what you simulated $\endgroup$ – StupidWolf Apr 9 '20 at 12:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.