[R] lm.ridge

[Berwin A Turlach]

G'day Daniel,

>>>>> "DR" == daniel  <daniel at sintesys.com.ar> writes:

    DR> First: I think coefficients from lm(Employed~.,data=longley)
    DR> should be equal coefficients from
    DR> lm.ridge(Employed~.,data=longley, lambda=0) why it does not
    DR> happen?
Which version of R and which version of MASS are you using?

> lm(Employed~.,data=longley)

Call:
lm(formula = Employed ~ ., data = longley)

Coefficients:
 (Intercept)  GNP.deflator           GNP    Unemployed  Armed.Forces  
  -3.482e+03     1.506e-02    -3.582e-02    -2.020e-02    -1.033e-02  
  Population          Year  
  -5.110e-02     1.829e+00  

> lm.ridge(Employed~.,data=longley, lambda=0)
               GNP.deflator           GNP    Unemployed  Armed.Forces 
-3.482259e+03  1.506187e-02 -3.581918e-02 -2.020230e-02 -1.033227e-02 
   Population          Year 
-5.110411e-02  1.829151e+00 


These coefficients look pretty identical to me, except that they are
printed to different numbers of significant digits.

In fact, the following shows that they are identical (upto numerical
precision): 

> fm1 <- lm(Employed~.,data=longley)
> fm2 <- lm.ridge(Employed~.,data=longley, lambda=0)
> coef2 <- print(fm2)
               GNP.deflator           GNP    Unemployed  Armed.Forces 
-3.482259e+03  1.506187e-02 -3.581918e-02 -2.020230e-02 -1.033227e-02 
   Population          Year 
-5.110411e-02  1.829151e+00 
> max(abs(coef(fm1)-coef2))
[1] 7.275958e-12


    DR> Second: if I have for example Ridge<-lm.ridge(Employed~.,
    DR> data=longley, lambda = seq(0,0.1,0.001)), I suppose intercept
    DR> coefficient is defined implicit,
Yes.

    DR> why it does not appear in Ridge$coef?
If you look at the code of lm.ridge, you will see that, if an
intercept is included in the model, all non-constant regressors are
centered (i.e. made orthogonal to the intercept term) and scaled to
have the same variance.  Further more, the intercept term is typically
*not* penalised.  The components in Ridge$coef are the coefficients on
this transformed scale.  No need of including the intercept here,
since it is the same for all values of lambda.  If you print the
model, then the ridge coefficients on the original scale are
calculated, see:

> getAnywhere("print.ridgelm")
A single object matching 'print.ridgelm' was found
It was found in the following places
  registered S3 method for print from namespace MASS
  namespace:MASS
with value

function (x, ...) 
{
    scaledcoef <- t(as.matrix(x$coef/x$scales))
    if (x$Inter) {
        inter <- x$ym - scaledcoef %*% x$xm
        scaledcoef <- cbind(Intercept = inter, scaledcoef)
    }
    print(drop(scaledcoef), ...)
}
<environment: namespace:MASS>

    DR> Third: I suppose that if I define
    DR> 1) y<-longley$Employed
    DR> 2) X<-as.matrix(cbind(1,Longley[,1:6])
    DR> 3) I = identity matrix the

    DR> following should be true: Coef=(X'X+kI)^(-1) X'y
No, as noted above, the intercept term is usually not penalised.

    DR> and if a take k=Ridge$kHKV, Coef should be approx equal to
    DR> Ridge$Coef[near value of kHKV]
No, as noted above the estimates in the "coef" component of an object
returned by lm.ridge are the coefficients on a different scale.

    DR> and it does not seem to happen, why?
Because the intercept is not penalised by lm.ridge and the
non-constant columns of the design matrix are rescaled; hence the
returned coefficients are on another scale.

    DR> Any help, suggestion or orientation?
HTH.

Cheers,

        Berwin