I'm using lm.ridge from the library MASS and am getting (I think) nonsense when I extract the
coef attribute of the fitted model. This may be because I'm misunderstanding the statistical issues at play here, so I'm unsure whether this is a computing problem or a statistical problem. This is a widely used R package and one that is well vetted, so I'm guessing the problem is my understanding of what it's doing.
library(MASS) x1 <- rnorm(200) x2 <- rnorm(200) y <- x1+x2+rnorm(200) # Ordinary least squares. mod1 <- lm(y~x1+x2) coef(mod1) (Intercept) x1 x2 -0.06253124 1.01663166 0.99337139 mod1$coef (Intercept) x1 x2 -0.06253124 1.01663166 0.99337139
First, I think when $\lambda=0$, you get make the OLS estimator, yes? Well...
mod2 <- lm.ridge(y ~ x1+x2, lambda=0) mod2$coef x1 x2 0.9686132 1.0735558 coef(mod2) x1 x2 -0.06253124 1.01663166 0.99337139
So... why would the two different approaches to extracting the coefficients give different results? I did check whether the first call corresponded to leaving the intercept out of the model (because the intercept is not printed for some reason). No, that's not the explanation. Leaving out the intercept from the OLS model gives a whole other answer.
What I eventually want to use this for is to select the optimally penalized coefficient estimates, but the two totally different sets of estimates confuses me. E.g.
mod3 <- lm.ridge(y ~ x1+x2, lambda=c(0.0, 0.1, 0.2) ) coef(mod3) x1 x2 0.0 -0.06253124 1.016632 0.9933714 0.1 -0.06255998 1.016044 0.9928101 0.2 -0.06258869 1.015457 0.9922495 mod3$coef 0.0 0.1 0.2 x1 0.9686132 0.9680532 0.9674939 x2 1.0735558 1.0729493 1.0723434
Say I know that $\lambda=0.1$ is optimal (e.g. using the
select function). Which coefficient estimates do I report? I guess the ones obtained from the
coef() function, because those check out in the one case I can verify ($\lambda=0$)? If that's the case, should I let the
MASS people know that their
$coef component is messed up?