[R] different L2 regularization behavior between lrm, glmnet, and penalized?

[Robert V (Bob) Sasseen]

The author of the penalized package, j.j.goeman at lumc.nl, kindly 
replied to my message. He also responded to another question I asked him.

------------------

The differences have to do with different 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.

------------------

I have another question about the penalized package. I'm not clear on 
the following behavior of penalized on the same data from my original 
question:

 > fit = penalized(z ~ x, lambda2 = 1)
 > penalty(fit)
        L1         L2
0.00000000 0.08492619

 > fit = penalized(z ~ x, lambda1 = 1, lambda2 = 1)
 > penalty(fit)
       L1        L2
0.3417345 0.1167825

What is the penalty function showing? I'd naively expect it to show the 
same lambda parameters that were given to the penalized() calls. I 
haven't been able to find this in any documentation.

------------------

The penalty returns the evaluated penalty term, lambda1*sum(abs(coef)) 
for L1; lambda2*sum(coef^2) for L2.
 
I'll make a note of documenting this.