# lm.ridge returns different results that are from manual calculation

I am comparing the manual calculation of ridge regression again lm.ridge function, however, it seems like the answer of the 2 techniques do not match. They only match when $\lambda=0$.

set.seed(1)

x <- rnorm(1000,1,2)
x <- matrix(x,ncol=10,nrow=100)
y <- rnorm(100,2,5)

xs <- scale(x,T,T)
ys <- scale(y,T,T)

p <- dim(x)[2]
lam <- 2
# manual Calculation
bh <- solve(t(xs) %*% xs + lam * diag(p), t(xs) %*% ys)
# lm..ridge
fit <- lm.ridge(ys~xs-1, lambda=lam)
coef_fit <- as.matrix(coef(fit),nco1)
cbind(bh, coef_fit)


Does anyone know why the estimated coefficients only match when $\lambda=0$ but not for other values of $\lambda$ s?

Update: Apologies for not scaling the data before. I have scaled the data now but there are still discrepancies between the estimated coefficients.

> cbind(bh, coef_fit)
[,1]         [,2]
xs1  -0.144767582 -0.144799855
xs2  -0.114627840 -0.114652989
xs3  -0.019612430 -0.019612567
xs4   0.007292303  0.007293982
xs5   0.044335298  0.044354816
xs6  -0.034135483 -0.034137483
xs7   0.020260806  0.020265217
xs8   0.058511001  0.058520197
xs9  -0.124643955 -0.124671909
xs10  0.060076729  0.060097567

• Ridge regression is carried out on standardized versions of the variables. You do not appear to standardize yours.
– whuber
Commented Jul 4, 2017 at 13:47
• Should I move this question to stackoverflow ? But somehow I feel that this question is more relevant to cross validated because it requires understanding how lm.ridge estimates statistically the coefficients. Commented Jul 4, 2017 at 16:03
• Please note that lm.ridge most likely is the version in the MASS package. It standardizes the explanatory variables in a slightly different way.
– whuber
Commented Jul 4, 2017 at 16:08
• So where can I read up on this new way of standardising the explanatory variables because ?lm.ridge does not seem to say anything about this. Commented Jul 4, 2017 at 16:17
• It's right there in the code. See the lines Xscale <- drop(rep(1/n, n) %*% X^2)^0.5; X <- X/rep(Xscale, rep(n, p)) They divide each column by its root mean square. If you do that, you will get agreement between your method and this to about 13-14 decimal places. The remaining two or three significant figures of disagreement are floating point roundoff error.
– whuber
Commented Jul 4, 2017 at 18:30

To deal with scaling problem, I suggest you replace lambda with lambda*(n/(n-1)). This will resolve the discrepancies.

In your example, it would be lam*(100/99).

set.seed(1)
x <- rnorm(1000,1,2)
x <- matrix(x,ncol=10,nrow=100)
y <- rnorm(100,2,5)

xs <- scale(x,TRUE,TRUE)
ys <- scale(y,TRUE,TRUE)

p <- dim(x)[2]
lam <- 2
# manual Calculation
bh <- solve(t(xs) %*% xs + lam * diag(p), t(xs) %*% ys)
# lm..ridge
fit <- lm.ridge(ys~xs-1, lambda=lam*(100/99))
coef_fit <- as.matrix(coef(fit),nco1)

cbind(bh, coef_fit)


This produces exactly the same two columns as below.

             [,1]         [,2]
xs1  -0.144767582 -0.144767582
xs2  -0.114627840 -0.114627840
xs3  -0.019612430 -0.019612430
xs4   0.007292303  0.007292303
xs5   0.044335298  0.044335298
xs6  -0.034135483 -0.034135483
xs7   0.020260806  0.020260806
xs8   0.058511001  0.058511001
xs9  -0.124643955 -0.124643955
xs10  0.060076729  0.060076729


edit1 How does it work?

Figuring out how this lm.ridge() helps us to understand it.

It is assumed that x has been scaled and centered

The manual calculation of coefficients is $\hat \beta = ( X'X+\lambda I)^{-1} X'Y$.

lm.ridge() uses scaled X $X_s=\sqrt{\frac{n}{n-1}}X$

fit$coef returns the vector :$\begin{align} \hat \beta &= ( X_s'X_s+\lambda I)^{-1} X_s'Y\\ &= \left( X'X+\frac{n-1}{n}\lambda I\right)^{-1} X'Y \sqrt{\frac{n-1}{n}} \end{align}$. fit$scale returns the vector whose all elements are $\sqrt{\frac{n-1}{n}}$ and has the same length with $\hat \beta$

And coef(fit) returns the value of fit$coef/fit$scale .

So to get the same value, you should replace lambda with lambda*(n/(n-1)).