Dear R-help,
I'm having trouble understanding how glmnet converts its coefficient estimates back to the original scale. Here's an example with ridge regression: ######################################################################## library(glmnet) set.seed(1) n <- 20 # sample size d <- data.frame(x1 = rnorm(n, 1, 1), x2 = rnorm(n, 10, 2), y = rnorm(n, 1, 2)) # Sample means mx1 <- mean(d$x1) mx2 <- mean(d$x2) my <- mean(d$y) # Scaling factors ("1/n standard deviations") sx1 <- sd(d$x1)*sqrt((n - 1)/n) sx2 <- sd(d$x2)*sqrt((n - 1)/n) sy <- sd(d$y)*sqrt((n - 1)/n) # Scaled variables d$x1s <- (d$x1 - mx1)/sx1 d$x2s <- (d$x2 - mx2)/sx2 d$ys <- (d$y - my)/sy lam <- 0.5 # lambda value # Using scaled variables (same result for standardize=TRUE and standardize=FALSE) glmnet1 <- glmnet(as.matrix(d[, c("x1s", "x2s")]), d$ys, alpha=0, lambda = lam) # Using unscaled variables glmnet2 <- glmnet(as.matrix(d[, c("x1", "x2")]), d$y, alpha=0, lambda=lam) coef(glmnet2) ## s0 ## (Intercept) 2.5658491 ## x1 0.3471199 ## x2 -0.1703715 # I want to calculate the glmnet2 coef estimates from the glmnet1 coef estimates. # The following attempts are based on rearrangement of # (y - my)/sy = beta1*(x1 - mx1)/sx1 + beta2*(x2 - mx2)/sx2 my - coef(glmnet1)["x1s", "s0"]*mx1*sy/sx1 - coef(glmnet1)["x2s", "s0"]*mx2*sy/sx2 # 2.430971 # Not the same as coef(glmnet2)["(Intercept)", "s0"] coef(glmnet1)["x1s", "s0"]*sy/sx1 # 0.3096897 # Not the same as coef(glmnet2)["x1", "s0"] coef(glmnet1)["x2s", "s0"]*sy/sx2 # -0.1524043 # Not the same as coef(glmnet2)["x2", "s0"] ###################################################################### I can apply a similar method (with centring of y instead of standardisation) to successfully get the coefficient estimates on the original scale given by lm.ridge in the MASS package. I would appreciate any help anyone can give on where I'm going wrong with glmnet. Thanks, Mark ______________________________________________ [hidden email] mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. |
This is interesting, can you post your lm.ridge solution as well? I
suspect in glmnet, you need to use model.matrix with intercept, that could be the reason. -m ______________________________________________ [hidden email] mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. |
Thanks for your reply Mehmet. I've found that the problem was that I
didn't scale the lambda value. My original example did not follow the instruction not to give a single lambda value, but that in itself wasn't the problem. Example shown below. library(glmnet) library(MASS) set.seed(1) n <- 20 d <- data.frame(x1 = rnorm(n, 1, 1), x2 = rnorm(n, 10, 2), y = rnorm(n, 1, 2)) # Sample means mx1 <- mean(d$x1) mx2 <- mean(d$x2) my <- mean(d$y) # Scaling factors sx1 <- sd(d$x1)*sqrt((n-1)/n) sx2 <- sd(d$x2)*sqrt((n-1)/n) sy <- sd(d$y)*sqrt((n-1)/n) # Scaled variables d$x1s <- (d$x1 - mx1)/sx1 d$x2s <- (d$x2 - mx2)/sx2 d$ys <- (d$y - my)/sy # Centred y d$yc <- d$y - my lam <- 1 # lambda value for lm.ridge lmr1 <- lm.ridge(y ~ x1 + x2, data=d, lambda=lam) lmr2 <- lm.ridge(yc ~ x1s + x2s, data=d, lambda=lam) coef(lmr1) my - coef(lmr2)["x1s"]*mx1/sx1 - coef(lmr2)["x2s"]*mx2/sx2 # same as coef(lmr1)[1] coef(lmr2)["x1s"]/sx1 # same as coef(lmr1)["x1"] coef(lmr2)["x2s"]/sx2 # same as coef(lmr1)["x2"] glmnet1 <- glmnet(as.matrix(d[, c("x1", "x2")]), d[, "y"], alpha=0) glmnet2 <- glmnet(as.matrix(d[, c("x1s", "x2s")]), d[, "ys"], alpha=0) # Note: glmnet1$lambda is glmnet2$lambda*sy ind <- 80 # index of lambda values to look at coef(glmnet1)[, ind] my - coef(glmnet2)["x1s", ind]*mx1*sy/sx1 - coef(glmnet2)["x2s", ind]*mx2*sy/sx2 # same as coef(glmnet1)["(Intercept)", ind] coef(glmnet2)["x1s", ind]*sy/sx1 # same as coef(glmnet1)["x1", ind] coef(glmnet2)["x2s", ind]*sy/sx2 # same as coef(glmnet1)["x2", ind] On Sat, Apr 4, 2015 at 6:03 AM, Suzen, Mehmet <[hidden email]> wrote: > This is interesting, can you post your lm.ridge solution as well? I > suspect in glmnet, you need to use model.matrix with intercept, that > could be the reason. > > -m ______________________________________________ [hidden email] mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. |
Free forum by Nabble | Edit this page |