6
$\begingroup$

Following is the ridge regression example in MASS package:

> head(longley)
        y     GNP Unemployed Armed.Forces Population Year Employed
1947 83.0 234.289      235.6        159.0    107.608 1947   60.323
1948 88.5 259.426      232.5        145.6    108.632 1948   61.122
1949 88.2 258.054      368.2        161.6    109.773 1949   60.171
1950 89.5 284.599      335.1        165.0    110.929 1950   61.187
1951 96.2 328.975      209.9        309.9    112.075 1951   63.221
1952 98.1 346.999      193.2        359.4    113.270 1952   63.639
> 
> mod = lm.ridge(y ~ ., longley,     lambda = seq(0,0.1,0.001))
> 
> plot(mod)

Following is the plot:

enter image description here

How do I interpret it? I understand these lines are for different independent variables but I want to know which of the independent variables are significant predictors of y in above dataset (i.e. I am interested in explanation and not prediction).

I tried to read up on internet but could not understand how to proceed.

Edit: Also what is the output of select(mod):

> select(mod)
modified HKB estimator is 0.006836982 
modified L-W estimator is 0.05267247 
smallest value of GCV  at 0.006 
$\endgroup$
2
$\begingroup$

The ridge regression will penalize your coefficients, such that those who are the least efficient in your estimation will "shrink" the fastest. Imagine you have a budget allocated and each coefficient can take some to play a role in the estimation. Naturally those who are more important will take more of the budget. As you increase the Lambda, you are "decreasing" the budget, aka penalizing more. For your plot, each line represents a coefficient whose value is going to zero as you are decreasing the budget or as you are penalizing more(increasing the lambda). To choose the best lambda though you should consult the MSE vs lambda plot. I would say though the faster a coefficient is shrinking the less important it is in prediction; e.g. I think the dotted dashed blue one should have more information than the solid black one. Try plotting a summary, a legend and an MSE vs lambda too. If you choose your best lambda and then look at your betas you can see which betas are more important by looking at their values at the optimum lambda.

$\endgroup$
  • $\begingroup$ I have added the legend. I am trying to get MSE vs lambda plot. I understand that I have to select the best lambda and beta for that lambda are the final result. $\endgroup$ – rnso May 30 '15 at 9:07
  • $\begingroup$ If i remember correctly you can use glmnet and use cv.glmnet for best lambda. Be sure to set your alpha correctly to get ridge and not lasso. Consult the glmnet documentation. $\endgroup$ – plumSemPy May 30 '15 at 9:20
  • $\begingroup$ I cannot find any MSE in the output. Are they same as GCV? GCV is lowest at lambda of about 0.005. $\endgroup$ – rnso May 30 '15 at 9:27
  • $\begingroup$ I do not know and I cannot seem to find a documentation on it. Give glmnet a try. It has way better documentation and examples. cran.r-project.org/web/packages/glmnet/glmnet.pdf $\endgroup$ – plumSemPy May 30 '15 at 9:34
  • $\begingroup$ I believe coefficient are best at minimum GCV, but how can I get their confidence intervals or their p values to indicate if they are significant predictors? $\endgroup$ – rnso May 30 '15 at 15:34
6
$\begingroup$

I know you say

I am interested in explanation and not prediction

but this is a road fraught with dangers. A model that does not predict its response well is not useful for explanation. Explanation and prediction are not independent desires, one should complement and enhance the other at all times, in both directions.

So given that, if you want to draw some explanation from a ridge regression, first find the most predictive model. Tune your regularization parameter $\lambda$ with cross validation or a hold out data set, there is lots of fantastic advice in this forum on how to do that. Once you have the most predictive $\lambda$, examine the coefficients of that model. If you want confidence intervals on these coefficients, bootstrap fit your fixed $\lambda$ model on your dataset, and empirically estimate the variance in the parameters.

Also what is the output of select(mod)?

The documentation of ridgelm is spotty, and doesn't say what select does, but you can get some clue from:

GCV: vector of GCV values

kHKB: HKB estimate of the ridge constant

kLW: L-W estimate of the ridge constant

So it looks like it's pulling out these three attributes (well the minimum of the first attribute) and wrapping them in a nice display. I can't recover the value it has for the minimum GCV though:

> min(mod$GCV) / nrow(longley)
[1] 0.007473027

so that remains a mystery.

$\endgroup$
  • $\begingroup$ why would you not bootstrap the crossvalidation step as well? $\endgroup$ – seanv507 Feb 24 '19 at 11:55
  • $\begingroup$ You probably should, I actually just re-read this last night and had the same thought! Not doing so estimates the variance of the test error given a fixed lambda, doing so estimates the variance of the test error of the entire procedure. You're probably right that the second thing is more honest. $\endgroup$ – Matthew Drury Feb 24 '19 at 18:03
2
$\begingroup$

Better check this?

http://www.mathworks.com/help/stats/ridge.html

There is a better picture an several more examples..... Though i dont know if you have Matlab...

EDIT - Command execution

Both the plot on the question and the image on the link show the ridge traces, showing the $k$ regularization coefficient on the x-axis, and the estimated coefficients on the y-axis.

As reference, the ridge regresion plot comes by doing regularization on linear regression $b=AX$ for improving conditioning, by placing a Ridge (Tikhonov Regularization) Matrix $K$: $$x_e=(A^TA+K^TK)^{-1}A^Tb$$

When having $K=kI$, we reduce the matrix selection to a scalar $k^2I$, which is done in both Matlab and R packages.

I tried to reproduce the results, for the same data, but apparently R do not make regularization on the constant term, and Matlab does it. Or R is lying you....

Ridge Example II

Regarding the question, it appears none of these algorithms deliver the p-values, the t-test or the MSE traces, so a coefficient assessment can't be done right away from there without an additional piece of code...

$\endgroup$
  • $\begingroup$ While this link may answer the question, it is better to include the essential parts of the answer here and provide the link for reference. Link-only answers can become invalid if the linked page changes. $\endgroup$ – gung - Reinstate Monica May 30 '15 at 8:45
  • $\begingroup$ No problem... lets wait for mso opinion. If he dont like it i will just delete it... $\endgroup$ – Brethlosze May 30 '15 at 8:47
  • $\begingroup$ The link is helpful. Thanks. Though i am using R. It will certainly be more helpful if you can add a few lines of important points. $\endgroup$ – rnso May 30 '15 at 9:02
2
$\begingroup$

One thing to watch out for is when a coefficient crosses zero, as Year does in the example.

If you made future predictions without regularization $(\lambda = 0)$, you have a negative trend over time; however, this is an artifact of multicollinearity in the data.

With not very much regularization $(\lambda > 0.005)$, the effect reverses, so we would expect $y$ to increase in the future, everything else being equal.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.