Solved – how to calculate effective degrees of freedom in ridge regression in R

lmrridge regression

Here is example to workout:

require(lasso2)
data(Prostate)
fit = lm(lpsa~.,data=Prostate)
summary(fit)

lam = seq(0,10000,len=5000)
require(MASS)
ridgefits = lm.ridge(lpsa~.,data=Prostate,lam=lam)

# plot traces for each predictor
plot(range(lam), range(ridgefits$coef),type="n")
clrs <- c("red", "blue", "green", "pink", "green1", "tan", "yellow", "gray40")
for(i in 1:nrow(ridgefits$coef)){
   lines(lam,ridgefits$coef[i,], col = clrs[i])
   }
legend(6500, 0.7, rownames(ridgefits$coef), col = clrs,  lwd = 2)

Next I would like to generate graph something like this –
enter image description here

I think I see nice discussion and Matlab script in this Q/A. Is there anyway we can do same in R ? Once we get this how can we get AIC /BIC and use that in decision making in selection of $\lambda$ ?

Edits: based on following suggestions in answer:

    require(rms)
    ridgefits = ols(lpsa~lcavol+lweight+age+lbph+svi+lcp+gleason+pgg45,
       method="qr", data=Prostate,
    se.fit = TRUE, x=TRUE, y=TRUE)
    p <- pentrace(ridgefits, seq(.2,1,by=.05))
    effective.df(ridgefits,p)

     out <- p$results.all
   out
   penalty       df      aic      bic    aic.c
1     0.00 8.000000 87.15916 66.56148 85.52280
2     0.20 7.994347 87.16845 66.58532 85.53437
3     0.25 7.992884 87.17026 66.59089 85.53677
4     0.30 7.991400 87.17186 66.59631 85.53898
5     0.35 7.989898 87.17327 66.60159 85.54099
6     0.40 7.988376 87.17448 66.60671 85.54282
7     0.45 7.986836 87.17549 66.61170 85.54445
8     0.50 7.985277 87.17632 66.61654 85.54591
9     0.55 7.983700 87.17696 66.62124 85.54719
10    0.60 7.982104 87.17741 66.62580 85.54829
11    0.65 7.980491 87.17768 66.63022 85.54921
12    0.70 7.978860 87.17777 66.63451 85.54995
13    0.75 7.977212 87.17768 66.63867 85.55053
14    0.80 7.975546 87.17742 66.64269 85.55094
15    0.85 7.973864 87.17698 66.64658 85.55118
16    0.90 7.972165 87.17637 66.65035 85.55125
17    0.95 7.970449 87.17559 66.65399 85.55117
18    1.00 7.968717 87.17464 66.65750 85.55092

    plot(out$penalty, out$df)
    plot( out$df, out$penalty, type = "l", col = "red")

I am close know but not still could not plot the figure exactly like above.

enter image description here

Best Answer

In conjunction with the ols function in the R rms package you can use the effective.df function to get what you want. This also tells you how many d.f. are effectively going to different types of terms (nonlinear, interaction, nonlinear interaction). The pentrace function helps you choose a penalty, or more than one penalty. For example, you may decide to penalize complex terms more than simple linear terms.

Related Question