d = read.table("../data/prostate.data",header=T)
fit.lm = lm(lpsa~.-train,data=d,subset=(train==TRUE))
summary(fit.lm)
#Covariates needs to be scale in order to reproduce Table 3.2
for(j in c(1,2,3,4,5,6,7,8))
 d[,j] = scale(d[,j])
fit.lm2 = lm(lpsa~.-train,data=subset(d,train==TRUE))
summary(fit.lm2)

fit.lm3 = lm(lpsa~.-train-gleason,data=subset(d,train==TRUE))
summary(fit.lm3)



library(leaps)
#Full search
system.time({regfit.full = regsubsets(lpsa~.-train,data=subset(d,train==TRUE))})
reg.summary = summary(regfit.full)
show(reg.summary$rsq)
par(mfrow = c(1,1))
plot(reg.summary$rss,xlab="Number of Variables",ylab="RSS",type="l")
plot(reg.summary$adjr2,xlab="Number of Variables",
         ylab="Adjusted RSq",type="l")
q = which.max(reg.summary$adjr2)
show(reg.summary$which[q,])

#Forward search
regfit.frw = regsubsets(lpsa~.-train,data=subset(d,train==TRUE),method="forward")
reg.summary = summary(regfit.frw)
q = which.max(reg.summary$adjr2)
show(reg.summary$which[q,])

#Backwards search
regfit.back = regsubsets(lpsa~.-train,data=subset(d,train==TRUE),method="backward")
reg.summary = summary(regfit.back)
q = which.max(reg.summary$adjr2)
show(reg.summary$which[q,])