11

I applied some data to find the best variables solution of regression model using ridge regression in R. I have used lm.ridge and glmnet (when alpha=0), but the results are very different especially when lambda=0. It suppose that both parameter estimators have the same values. So, what is the problem here? best regards

Andy
  • 18,070
  • 20
  • 77
  • 100

1 Answers1

15

glmnet standardizes the y variable and uses the mean squared errors instead of sum of squared errors. So you need to make the appropriate adjustments to match their outputs.

library(ElemStatLearn)
library(glmnet)
library(MASS)

dof2lambda <- function(d, dof) {
    obj <- function(lam, dof) (dof - sum(d ^ 2 / (d ^ 2 + lam))) ^ 2
    sapply(dof, function(x) optimize(obj, c(0, 1e4), x)$minimum)
}

lambda2dof <- function(d, lam) {
    obj <- function(dof, lam) (dof - sum(d ^ 2 / (d ^ 2 + lam))) ^ 2
    sapply(lam, function(x) optimize(obj, c(0, length(d)), x)$minimum)
}

dat   <- prostate
train <- subset(dat,  train, select = -train)
test  <- subset(dat, !train, select = -train)

train.x <- as.matrix(scale(subset(train, select = -lpsa)))
train.y <- as.matrix(scale(train$lpsa))

d   <- svd(train.x)$d
dof <- seq(1, 8, 0.1)
lam <- dof2lambda(d, dof)

ridge1 <- lm.ridge(train.y ~ train.x, lambda = lam)
ridge2 <- glmnet(train.x, train.y, alpha = 0, lambda = lam / nrow(train.x))

matplot(dof, t(ridge1$coef), type = 'l')
matplot(lambda2dof(d, ridge2$lambda * nrow(train.x)), t(ridge2$beta), type = 'l')
  • Hello. What R package do you advise to use for Ridge Regression? glmnet, bigRR, Mass, other? Any of them able to deal with repeated measures (random effects)? – skan Jun 04 '16 at 18:48
  • I've got some discrepancy between glmnet and MASS lm.ridge output that can't be explained by rescaling issues. However, lm.ridge gave me the result that coincides with hand calculation. If I have time in the future, I will post the complete example. Right now, I will go with MASS. Also, glmnet itself warns that results depend on the way you set up lambda parameters, that's one more strike against glmnet. – PA6OTA Feb 25 '18 at 01:00