Re: [R] lm.ridge

From: Berwin A Turlach <berwin_at_maths.uwa.edu.au>
Date: Thu 25 Aug 2005 - 20:18:43 EST

G'day Daniel,

>>>>> "DR" == daniel <daniel@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



R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Thu Aug 25 20:23:24 2005

This archive was generated by hypermail 2.1.8 : Fri 03 Mar 2006 - 03:39:55 EST