Dear Robert,
The differences have to do with diffent scaling defaults.
lrm by default standardizes the covariates to unit sd before applying
penalization. penalized by default does not do any standardization, but
if asked standardizes on unit second central moment. In your example:
x = c(-2, -2, -2, -2, -1, -1, -1, 2, 2, 2, 3, 3, 3, 3)
z = c(0, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 1, 1)
You got:
> g = lrm(z ~ x, penalty = 1)
> g$coefficients
Intercept x
-0.1620727 0.3241454
Now:
> coef(penalized(z ~ x, lambda2 = 1,standardize=T))
(Intercept) x
-0.1651565 0.3303130
is already more similar. To counter the difference between dividing by n
(penalized) or n-1 (lrm), do
> coef(penalized(z ~ x, lambda2 = 14/13,standardize=T))
(Intercept) x
-0.1620727 0.3241454
The glmnet case should be similar, but I don't know the details here. It
seems glmnet uses its own peculiar way of defining the penalty, but some
choice of scaling should be able to bring glmnet in line as well.
Kind regards,
Jelle
________________________________
From: Robert V (Bob) Sasseen [mailto:sass...@ai.sri.com]
Sent: 29 October 2009 19:32
To: Goeman, J.J. (MSTAT)
Subject: different L2 regularization behavior between lrm,
glmnet, and penalized? (original question)
Jelle,
Below is the original question. The formatting apparently got a
bit messed up when it was parsed by the list server. This may be easier
to read.
Bob
------------------------------------------------------------------------
---------
The following R code using different packages gives the same
results for a simple logistic regression without regularization, but
different results with regularization. This may just be a matter of
different scaling of the regularization parameters, but if anyone
familiar with these packages has insight into why the results differ,
I'd appreciate hearing about it. I'm new to R. Thanks. (Version info
below. Same results on Windows and Solaris 8, except that I haven't
gotten glmnet to compile on the latter.)
Robert V (Bob) Sasseen
sass...@ai.sri.com
> # Some x values (predictive var).
> x = c(-2, -2, -2, -2, -1, -1, -1, 2, 2, 2, 3, 3, 3, 3)
> # Some z values--the observed outcome.
> # Effect is that for
> # x = -2, p = 1/4;
> # x = -1, p = 1/3;
> # x = 2, p = 2/3;
> # x = 3, p = 3/4.
> z = c(0, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 1, 1)
> library(Design)
> g = lrm(z ~ x)
> g$coefficients
Intercept x
-0.2224842 0.4449685
> g = lrm(z ~ x, penalty = 1)
> g$coefficients
Intercept x
-0.1620727 0.3241454
> library(glmnet)
> g = glmnet(cbind(x), cbind(1-z, z), family = "binomial",
lambda = 0, standardize = FALSE)
> coef(g)
1
-0.2224843
x 0.4449687
> g = glmnet(cbind(x), cbind(1-z, z), family = "binomial",
lambda = 1, alpha = 0, standardize = FALSE)
> coef(g)
1
-0.1098361
x 0.2196721
> library(penalized)
> fit = penalized(z ~ x)
> coefficients(fit, "all")
(Intercept) x
-0.2224843 0.4449687
> fit = penalized(z ~ x, lambda2 = 1)
> coefficients(fit, "all")
(Intercept) x
-0.2060658 0.4121315
> sessionInfo()
R version 2.9.2 (2009-08-24)
i386-pc-mingw32
locale:
LC_COLLATE=English_United States.1252;LC_CTYPE=English_United
States.1252;LC_MONETARY=English_United
States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252
attached base packages:
[1] splines stats graphics grDevices utils datasets
methods
[8] base
other attached packages:
[1] penalized_0.9-26 glmnet_1.1-3 Matrix_0.999375-30
lattice_0.17-25
[5] Design_2.3-0 Hmisc_3.7-0 survival_2.35-4
loaded via a namespace (and not attached):
[1] cluster_1.12.0 grid_2.9.2
--
Robert V (Bob) Sasseen
Senior Software Engineer
Advanced Analytics
Artificial Intelligence Center
SRI International
3661 Valley Centre Drive, Suite 150
San Diego, CA 92130
robert.sass...@sri.com
(858) 350-2035
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
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.
