lassoshooting and glmnet different results

Question

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

EDIT: According to the comment below, the new 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)

# 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.

Neal Fultz · Accepted Answer · 2018-07-10T21:39:05.740

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.

(+1) Thanks for your reply. I attached a new code, according to your answer (I have edited the standardization). I run the code two times. First `glmnet` and `lassoshooting` on the data standardized manually (the same results), second `glmnet` (its own standardization) and `lassoshooting` (manual standardization and scale back to original). Results are very different. What can be the reason? — mert, Jul 10 '18 at 06:43
Obviously `glmnet` has correct beta's when it standardizes the data by itself. — mert, Jul 10 '18 at 06:51

lassoshooting and glmnet different results

EDIT: According to the comment below, the new code

1 Answers1