I am trying to implement ridge regression in R, but the results are wrong.

$\hat\beta^{Ridge} = argmin\sum_{i=1}^N(y-\beta_0 - \sum_{j=1}^{p} x_{i,j}\beta_j)^2 + \lambda\sum_{j=1}^{p}\beta_j^2$


$RSS(\lambda) = (y - X\beta)^T(y-X\beta) + \lambda\beta^T\beta = |y|^2-2|X^Ty\beta|+|X\beta|^2 + \lambda|\beta|^2$

$dRSS(\lambda)/d\beta = -2X^Ty + 2X^TX\beta+2\lambda\beta$, so $2\beta(X^TX+\lambda I) = 2X^Ty => \hat\beta^{Ridge} = (X^TX+\lambda I)^{-1}X^Ty$

Now I want to use this in a code to get $\hat\beta$'s for Ridge regression

 y <- longley$GNP.deflator
 x <- longley[,2:ncol(longley)]
 X <- model.matrix(~as.matrix(x))
 BetaRidge <- function(lam){
         solve( crossprod(X) +
                diag(ncol(X))*lam) %*% crossprod(X,y)

     data=longley,lambda = 1)

The output is considerably different! What is wrong with the logic?

  • 2
    $\begingroup$ there is at least one error in your code and calculations: model.matrix will include an intercept, in ridge regression however you don't penalize the intercept. second thing which comes to my mind is that in penalized regression methods one often estimates the coefficients using standardized regressors (implying $\beta_0 = 0$). $\endgroup$
    – BloXX
    Oct 22 '18 at 10:02
  • $\begingroup$ I had tried without intercept too, but in that way if I remember it right the difference between betas was 10 times bigger. $\endgroup$
    – user13
    Oct 22 '18 at 11:17
  • $\begingroup$ I will try to standardize then use the model. $\endgroup$
    – user13
    Oct 22 '18 at 11:19

As said in the comments:

1.) you don't penalize the intercept in ridge regression and

2.) penalized least square methods like ridge regression or lasso are most meaningful applied to standardized regressors and a centered dependent variable (to get rid of the intercept). So, internally, lm.ridge will perform it's optimization task using $X_{i,j}^*$ and $y_i^*$, where $$y_i^* = y_i - \overline{y}$$ and $$x_{i,j}^* = (x_{i,j}-\overline{x}_j)/\widehat{\sigma}_j,$$ where (in the case of lm.ridge) $$\widehat{\sigma}_j^2 = \frac{1}{n}\sum_{i=1}^n (x_{i,j}-\overline{x}_j)^2.$$

Indeed: using standardized variables, your R-Code leads to identical results (besides some minor rounding errors, coming from different inversion methods):

 BetaRidge <- function(lam){solve( crossprod(X) +diag(ncol(X))*lam) %*% crossprod(X,y)}
 y <- longley$GNP.deflator
 X <- as.matrix(longley[,2:ncol(longley)])
 X<-diag(rep(sqrt(16/15),16))%*%scale(X) #standardization such that \sum_{i}(x_{i,j}^*)^2 = n

all(round(coef(lm.ridge(y~X,lambda = 1))[-1],12)== round(BetaRidge(1),12))


of course, after estimating the coefficients you need (could!) to transform them back on their original scale.


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.