I want to compare lassoshooting and glmnet for lasso.
No standardization option in lassoshooting; so i standardized data first, fit the model and re-standardized to original scale.
Results are different and it seems lassoshooting beta's is closer to the original beta's.
Do i have a mistake?
The code:
library(lassoshooting)
library(glmnet)
set.seed(327)
n = 500
p = 9
x = matrix(rnorm(n*p), ncol=p)
n = nrow(x)
b = c(.5, -.5, .25, -.25, .125, -.125, rep(0, 3))
y = x %*% b + rnorm(n, sd=.05)
xs = scale(x)
ys = scale(y)
lam = 0.1
glmnet_res = coef(glmnet(x, y), s=lam)[-1]
lassoshooting_res = lassoshooting(X=xs, y=ys, thr=1e-7, lambda=n*lam)$coefficients
# n in n*lam stems from difference between objective functions of two packages
# standard deviations for original scale
sds = apply(x,2,sd)
sdy = sd(y)
lasso_shooting_o = sdy*lassoshooting_res/sds
# compare
cbind(glmnet=glmnet_res, lassoshooting=lasso_shooting_o)
glmnet lassoshooting
[1,] 0.40123563 0.42270220
[2,] -0.38733635 -0.41195555
[3,] 0.14463257 0.16727953
[4,] -0.15914094 -0.17799495
[5,] 0.02942027 0.04958667
[6,] -0.01465437 -0.03777288
[7,] 0.00000000 0.00000000
[8,] 0.00000000 0.00000000
[9,] 0.00000000 0.00000000
# Is lassoshooting closer to true parameters ?
abs(lasso_shooting_o-b) <= abs(glmnet_res-b)
[1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
library(lassoshooting)
library(glmnet)
set.seed(327)
n = 500
p = 9
x = matrix(rnorm(n*p), ncol=p)
n = nrow(x)
b = c(5, -5, 25, -25, 125, -125, rep(0, 3))
y = x %*% b + rnorm(n, sd=.05)
# 1/n type standardization
xc = sweep(x, 2, colMeans(x))
sdc = sqrt(apply(xc, 2, crossprod)/nrow(x))
xs = sweep(xc, 2, sdc, "/")
ys = scale(y)*sqrt(n/(n-1))
lam = 0.1
# BOTH ARE STANDARDIZED: RESULTS ARE THE SAME
glmnet_std = coef(glmnet(xs, ys,standardize=F), s=lam)[-1]
lassoshooting_std = lassoshooting(X=xs, y=ys, thr=1e-7, lambda=n*lam)$coefficients
# n in n*lam stems from difference between objective functions of two packages
# compare
cbind(glmnet=glmnet_std, lassoshooting=lassoshooting_std)
#########################################################
glmnet lassoshooting
[1,] 0.00000000 0.00000000
[2,] 0.00000000 0.00000000
[3,] 0.04224107 0.04224107
[4,] -0.04178765 -0.04178765
[5,] 0.59188462 0.59188462
[6,] -0.59781943 -0.59781943
[7,] 0.00000000 0.00000000
[8,] 0.00000000 0.00000000
[9,] 0.00000000 0.00000000
#########################################################
# glmnet on ORIGINAL DATA with its own standardization
# lassoshooting on SCALED DATA (THEN RE-SCALED)
# VERY DIFFERENT RESULTS (selected variables are different too)
glmnet_o = coef(glmnet(x, y), s=lam)[-1]
# original scale
sdy = sd(y)/sqrt(n/(n-1))
lasso_shooting_o = sdy*lassoshooting_std/sdc # sdc is defined above
# compare
cbind(glmnet=glmnet_o, lassoshooting=lasso_shooting_o)
#########################################################
glmnet lassoshooting
[1,] 2.742806 0.000000
[2,] -2.412466 0.000000
[3,] 22.618378 7.379036
[4,] -23.014604 -7.205713
[5,] 122.877714 106.617676
[6,] -122.566183 -105.931439
[7,] 0.000000 0.000000
[8,] 0.000000 0.000000
[9,] 0.000000 0.000000
#########################################################
# OBVIOUSLY glmnet is correct.
note that lassoshooting gives the same results (to several digits anyway) as glmnet(xs, ys, standardize=FALSE), so really what we are interested in is why internal standardizing is different than external:
> coef(glmnet(xs, ys, intercept=FALSE, standardize=FALSE), s=lam)[-1]
[1] 0.54023720669 -0.51377928289 0.21423980260 -0.23094074895 0.06158780181
[6] -0.04769218136 0.00000000000 0.00000000000 0.00000000000
> lassoshooting(X=xs, y=ys, thr=1e-7, lambda=n*lam)$coefficients
[1] 0.54023682002 -0.51377917401 0.21423976696 -0.23094082042 0.06158781750
[6] -0.04769218772 0.00000000000 0.00000000000 0.00000000000
Internally, when glmnet standardizes, it uses a denominator of n, while scale() is using n - 1. We can adjust for that ourself:
> coef(glmnet(x, ys * sqrt((n-1)/n)), s=lam)[-1] * sdy / sqrt((n-1)/n)
[1] 0.42270249793 -0.41195563736 0.16727956071 -0.17799489896 0.04958665639
[6] -0.03777287171 0.00000000000 0.00000000000 0.00000000000
EDIT:
Also keep in mind that glmnet is adding an intercept by default, so
glmnet(x, y, intercept=FALSE) is very different from glmnet(x, y) - this is also a little hard to compare since it looks like lassoshooting puts it's lambda param on a different scale? anyway, try
lassoshooting(X=cbind(1,xs), y=y, thr=1e-7, lambda=n*sqrt(n)*lam, nopenalize=0) * sdc
which is close but not exactly the same, as in ?lassoshooting's example about intercepts.