[R] obtaining p-values for lm.ridge() coefficients (package 'MASS')
Liviu Andronic
landronimirc at gmail.com
Wed Aug 24 13:44:57 CEST 2011
Dear Michael
Thank you for your pointers.
On Tue, Aug 23, 2011 at 4:05 PM, Michael Friendly <friendly at yorku.ca> wrote:
> First, you should be using rms::ols, as Design is old.
>
Good to know. I've always wondered why Design and rms, in many cases,
were providing similar functions and (upon cursory inspection)
identical documentation.
> Second, penalty in ols() is not the same as the ridge constant in lm.ridge,
> but rather a penalty on the log likelihood. The documentation
> for ols refers to ?lrm for the description.
>
I see. At the same time ?rms::ols contains: "The penalized maximum
likelihood estimate (penalized least squares or ridge estimate) of
beta is (X'X + P)^{-1} X'Y.", while ?rms::lrm defines P as "penalty
\times diag(pf) \times penalty.matrix \times diag(pf), where pf is the
vector of square roots of penalty factors computed from penalty". This
is suspiciously similar to the definition of the "classical" Ridge
Regression beta estimates (X'X + kI)^{-1} X'Y, where k is the lambda
constant and I the identity matrix.
Hence I was wondering if rms::ols could be forced to use 'P = kI', and
thus compute the "classical" Ridge estimates? I tried to hack rms::ols
to ensure that it passed 'kI' to lm.pfit(..., penalty.matrix = ...),
but I get strange results so I must be missing something. (See
attached.)
> Third, other ridge estimators (e.g., ElemStatLearn::simple.ridge also
> give slightly different results for the same data due to differences in the
> way scaling is handled.
>
And ?ElemStatLearn::prostate points to several other 'ridge'-related functions.
However, similar to MASS::lm.ridge, simple.ridge() roughly stops at
estimating the betas, hence my question: is the purpose of Ridge
Regression to stabilize the coefficients (when the design matrix is
collinear) and use the results to infer the relative importance of the
coefficients on the response, while keeping the unpenalized model for
prediction purposes? Or am I missing something?
Thank you
Liviu
More information about the R-help
mailing list