Re: [R] glmnet: converting coefficients back to original scale
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 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
__
R-help@r-project.org 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.
Re: [R] glmnet: converting coefficients back to original scale
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 __ R-help@r-project.org 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.
[R] glmnet: converting coefficients back to original scale
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
__
R-help@r-project.org 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.
