require(MASS)
prostate <- read.delim("../../data/prostate.txt")
prostate <- prostate[,-ncol(prostate)]
fit <- lm.ridge(lpsa~.,prostate,lambda=seq(0,100,len=501))
plot(fit)
fit$lambda[which.min(fit$GCV)]

## More interesting lambda values
lam <- c(0,exp(seq(log(0.01),log(1e8),len=201)))
fit <- lm.ridge(lpsa~.,prostate,lambda=lam)
## Calculating df
XX <- scale(model.matrix(lpsa~0+.,prostate),center=fit$xm,scale=fit$scales)
l <- eigen(crossprod(XX))$values
fit$df <- numeric(length(fit$lambda))
for (i in 1:length(fit$lambda)) fit$df[i] <- sum(l/(l+fit$lambda[i]))
## Cooler plot
matplot(c(fit$df,0),rbind(coef(fit)[,-1],0),type="l",lwd=2,lty=1,col="slateblue",xlab="Degrees of freedom",ylab="Coefficients",xlim=c(0,max(fit$df)+1))
x <- pretty(c(fit$df,0))
xx <- numeric(length(x))
xx[1] <- expression(infinity)
for (i in 2:(length(x)-1)) xx[i] <- round(fit$lambda[which.min(abs(fit$df-x[i]))])
xx[length(xx)] <- 0
axis(3,at=pretty(c(fit$df,0)),labels=xx)
mtext(expression(lambda),line=2.5)
text(max(x)+0.75,coef(fit)[1,-1],labels=colnames(XX))

## Model selection
n <- nrow(prostate)
X <- model.matrix(lpsa~.,prostate)
fit$RSS <- apply(prostate$lpsa-X%*%t(coef(fit)),2,crossprod)
AIC <- n*log(fit$RSS) + 2*fit$df
BIC <- n*log(fit$RSS) + log(n)*fit$df
abline(v=fit$df[which.min(AIC)],col="red")
abline(v=fit$df[which.min(BIC)],col="green")
abline(v=fit$df[which.min(fit$GCV)])

## Other plots
par(mfrow=c(1,2))
plot(fit$lambda,fit$df,log="x",type="l",xlab=expression(lambda),ylab="Degrees of freedom",lwd=3,xaxt="n")
axis(1,at=c(1e-2,1e1,1e4,1e7),labels=c(0.01,10,expression(10^4),expression(10^7)))
plot(fit$df,fit$RSS,type="l",xlab="Degrees of freedom",ylab="RSS",lwd=3)
plot(fit$df,AIC-min(AIC),type="l",ylim=c(0,20),lwd=3,col="red",ylab="Criterion - min",xlab="Degrees of freedom")
lines(fit$df,BIC-min(BIC),type="l",lwd=3,col="green")
g <- fit$GCV*n^2
lines(fit$df,g-min(g),type="l",lwd=3)

## Comparison with OLS
fit.ols <- lm(lpsa~.,prostate)
sig2 <- as.numeric(crossprod(fit.ols$residuals)/fit.ols$df.residual)
b.o <- coef(fit.ols)[-1]
se.o <- sqrt(diag(vcov(fit.ols)))[-1]

k <- which.min(fit$GCV)
XX <- cbind(1,scale(model.matrix(lpsa~0+.,prostate),center=fit$xm,scale=fit$scales))
b.r <- coef(fit)[k,-1]
l <- fit$lambda[k]
W <- solve(crossprod(XX)+diag(c(0,rep(l,ncol(XX)-1))))
V <- sig2*W%*%crossprod(XX)%*%W
se.r <- sqrt(diag(V))[-1]/fit$scales

z.o <- b.o/se.o
z.r <- b.r/se.r