#Hitters data
library(ISLR)
data(Hitters)
head(Hitters)
dim(Hitters)
sum(is.na(Hitters$Salary))
Hitters = na.omit(Hitters)
dim(Hitters)
n = nrow(Hitters)
fit.lm = lm(Salary~.,data=Hitters)
summary(fit.lm)
library(leaps)
system.time({regfit.full = regsubsets(Salary~.,data=Hitters)})

#Forward selection
nvmax=19
regfit.fwd = regsubsets(Salary~.,data=Hitters,nvmax=nvmax,method="forward")
summary.regfit.fwd <- summary(regfit.fwd)
show(summary.regfit.fwd)
plot(regfit.fwd,scale="r2")
plot(regfit.fwd,scale="adjr2")
matplot(cbind(summary.regfit.fwd$rsq,summary.regfit.fwd$adjr2),type="l")
plot(summary.regfit.fwd$bic,type="l")
summary.regfit.fwd$aic = summary.regfit.fwd$bic-log(n)*(c(1:nvmax)+2)+2*(c(1:nvmax)+2)
plot(summary.regfit.fwd$aic,type="l")
plot(summary.regfit.fwd$cp,type="l")
matplot(cbind(summary.regfit.fwd$cp,summary.regfit.fwd$aic,summary.regfit.fwd$bic),type="l",lty=1)
p.Cp = which.min(summary.regfit.fwd$cp)
p.aic = which.min(summary.regfit.fwd$aic)
p.bic = which.min(summary.regfit.fwd$bic)
points(c(p.Cp,p.aic,p.bic),
       c(summary.regfit.fwd$cp[p.Cp],summary.regfit.fwd$aic[p.aic],summary.regfit.fwd$bic[p.bic]))
legend("topright",c("C_p","AIC","BIC"),lty=1,col=1:3,)

nvmax=50
regfit.fwd = regsubsets(Salary~.^2-NewLeague^2,data=Hitters,nvmax=nvmax,method="forward")
summary.regfit.fwd <- summary(regfit.fwd)
plot(regfit.fwd,scale="r2")
plot(regfit.fwd,scale="adjr2")
matplot(cbind(summary.regfit.fwd$rsq,summary.regfit.fwd$adjr2),type="l",ylab="R^2 and adj-R^2")
plot(summary.regfit.fwd$bic,type="l")
which.min(summary.regfit.fwd$bic)
summary.regfit.fwd$aic = summary.regfit.fwd$bic-log(n)*(c(1:nvmax)+2)+2*(c(1:nvmax)+2)
plot(summary.regfit.fwd$aic,type="l")
which.min(summary.regfit.fwd$aic)
plot(summary.regfit.fwd$cp,type="l")
matplot(cbind(summary.regfit.fwd$cp,summary.regfit.fwd$aic,summary.regfit.fwd$bic),type="l",lty=1)
p.Cp = which.min(summary.regfit.fwd$cp)
p.aic = which.min(summary.regfit.fwd$aic)
p.bic = which.min(summary.regfit.fwd$bic)
points(c(p.Cp,p.aic,p.bic),
       c(summary.regfit.fwd$cp[p.Cp],summary.regfit.fwd$aic[p.aic],summary.regfit.fwd$bic[p.bic]))

matplot(cbind(summary.regfit.fwd$aic,summary.regfit.fwd$cp-2*n),type="l")

#Plot av ulike seleksjonskriterier
reg.summary = summary(regfit.fwd)
names(reg.summary)
matplot(cbind(1:nvmax,1:nvmax),
        cbind(reg.summary$rsq,reg.summary$adjr2),
        type="l",lty=1,col=1:2,
        xlab="d",ylab="R^2/R^2 adjusted")
m1 = which.max(reg.summary$adjr2)
points(m1,reg.summary$adjr2[m1],col=2)
par(new=T)
matplot(cbind(1:nvmax,1:nvmax),
        cbind(reg.summary$cp,reg.summary$bic),
        type="l",lty=1,col=3:4,xaxt="n",yaxt="n",xlab="",ylab="",ylim=c(-350,100))
axis(4)
mtext("C_p/BIC",side=4,line=3)
legend("bottomright",lt=1,col=1:4,legend=c("R^2","R2^ adjusted","C_p","BIC"))
m2 = which.min(reg.summary$cp)
points(m2,reg.summary$cp[m2],col=3)
which.min(reg.summary$bic)
m3 = which.min(reg.summary$bic)
points(m3,reg.summary$bic[m3],col=4)
mod = summary.regfit.fwd$which[m3,]
show(mod[mod==TRUE])


#Full subset selection
regfit.full = regsubsets(Salary~.,data=Hitters,nvmax=20)
summary(regfit.full)
reg.summary = summary(regfit.full)
matplot(cbind(1:19,1:19),
        cbind(reg.summary$rsq,reg.summary$adjr2),
        type="l",lty=1,col=1:2,
        xlab="d",ylab="R^2/R^2 adjusted")
par(new=T)
matplot(cbind(1:19,1:19),
        cbind(reg.summary$cp,reg.summary$bic),
        type="l",lty=1,col=3:4,xaxt="n",yaxt="n",xlab="",ylab="",ylim=c(-300,100))
axis(4)
mtext("C_p/BIC",side=4,line=3)
legend("bottomright",lt=1,col=1:4,legend=c("R^2","R2^ adjusted","C_p","BIC"))

which.max(reg.summary$adjr2)
which.min(reg.summary$cp)
which.min(reg.summary$bic)

plot(regfit.full,scale="r2")
plot(regfit.full,scale="adjr2")

#Backward selection
regfit.bwd = regsubsets(Salary~.,data=Hitters,nvmax=20,method="backward")

coef(regfit.full,7)
coef(regfit.fwd,7)
coef(regfit.bwd,7)

#Including second order interactions
nvmax=150
regfit.full = regsubsets(Salary~.^2-NewLeague^2,data=Hitters,nvmax=nvmax,method="forward")
summary(regfit.full)
plot(regfit.full,scale="r2")
plot(regfit.full,scale="adjr2")

#Plot av ulike seleksjonskriterier
reg.summary = summary(regfit.full)
names(reg.summary)
matplot(cbind(1:nvmax,1:nvmax),
        cbind(reg.summary$rsq,reg.summary$adjr2),
        type="l",lty=1,col=1:2,
        xlab="d",ylab="R^2/R^2 adjusted")
m1 = which.max(reg.summary$adjr2)
points(m1,reg.summary$adjr2[m1],col=2)
par(new=T)
matplot(cbind(1:nvmax,1:nvmax),
        cbind(reg.summary$cp,reg.summary$bic),
        type="l",lty=1,col=3:4,xaxt="n",yaxt="n",xlab="",ylab="",ylim=c(-350,100))
axis(4)
mtext("C_p/BIC",side=4,line=3)
legend("bottomright",lt=1,col=1:4,legend=c("R^2","R2^ adjusted","C_p","BIC"))
m2 = which.min(reg.summary$cp)
points(m2,reg.summary$cp[m2],col=3)
which.min(reg.summary$bic)
m3 = which.min(reg.summary$bic)
points(m3,reg.summary$bic[m3],col=4)
mod = reg.summary$which[m3,]
show(mod[mod==TRUE])

#Dividing into training and test set
set.seed(1)
train=sample(c(TRUE,FALSE), nrow(Hitters),rep=TRUE)
test=(!train)
regfit.best=regsubsets(Salary~.,data=Hitters[train,],nvmax=19)
test.mat=model.matrix(Salary~.,data=Hitters[test,])
val.errors=rep(NA,19)
for(i in 1:19){
  coefi=coef(regfit.best,id=i)
  pred=test.mat[,names(coefi)]%*%coefi
  val.errors[i]=mean((Hitters$Salary[test]-pred)^2)
}
val.errors
which.min(val.errors)
coef(regfit.best,10)

#Making function that do prediction for a regsubsets object
predict.regsubsets=function(object,newdata,id,...){
  form=as.formula(object$call[[2]])
  mat=model.matrix(form,newdata)
  coefi=coef(object,id=id)
  xvars=names(coefi)
  mat[,xvars]%*%coefi
}
coef(regfit.best,10)

#10-fold cross-validation
k=10
set.seed(1)
folds=sample(1:k,nrow(Hitters),replace=TRUE)
cv.errors=matrix(NA,k,19, dimnames=list(NULL, paste(1:19)))
for(j in 1:k){
  best.fit=regsubsets(Salary~.,data=Hitters[folds!=j,],nvmax=19)
  for(i in 1:19){
    pred=predict(best.fit,Hitters[folds==j,],id=i)
    cv.errors[j,i]=mean( (Hitters$Salary[folds==j]-pred)^2)
  }
}
mean.cv.errors=apply(cv.errors,2,mean)
mean.cv.errors
par(mfrow=c(1,1))
plot(mean.cv.errors,type='b')
reg.best=regsubsets(Salary~.,data=Hitters, nvmax=19)
coef(reg.best,11)