#===================================================
# Advertising data
#===================================================

# Read the data
advertising <- read.csv("http://www.uio.no/studier/emner/matnat/math/STK2100/data/advertising.csv")

# Fit a model with TV and Radio
fit.lm2 <- lm(Sales~TV+Radio,data=advertising)

# Model with interaction effect
fit.lm4 <- lm(Sales~TV+Radio+TV:Radio,data=advertising)
summary(fit.lm4)
anova(fit.lm2,fit.lm4)

# =============================================
# Prostate data
# =============================================
#
# Summary:
# This data set comes from the study by Stamey
# et al. (1989), that examined the relationship
# between the level of prostate specific antigen
# (PSA) in 97 men who were about to receive
# a radical prostatectomy.
#
# Data description:
# The data set consists of the following 9 variables:
# Response variable:
#  lpsa: log of PSA
# Covariates:
#  lcavol: log cancer volume
#  lweight: log prostate weight
#  age: age in years
#  lbph: log of benignprostatic hyperplasia amount
#  svi: seminal vesicle invasion
#  lcp: log of capsular penetration
#  gleason: Gleason score
#  pgg45: percent of Gleason scores 4 and 5
#
#======================================================

# Read the data
prostate <- read.csv("http://www.uio.no/studier/emner/matnat/math/STK2100/data/prostate.data",sep="\t",header=TRUE)
prostate <- prostate[,-1]

# Fit full model
fit.lm = lm(lpsa~.-train,data=prostate,subset=(train==TRUE))
summary(fit.lm)

library(leaps)

### Forward selection
nvmax=8
regfit.fwd = regsubsets(lpsa~.-train,data=subset(prostate,train==TRUE),nvmax=nvmax,method="forward")
summary.regfit.fwd <- summary(regfit.fwd)
show(summary.regfit.fwd)

# R^2 when including each of the covariates
plot(regfit.fwd,scale="r2")

#R^2 and R_adj^2 against the no. of cov. in the model
matplot(cbind(summary.regfit.fwd$rsq,summary.regfit.fwd$adjr2),type="l")

#C_p, AIC and BIC against the no. of cov. in the model
n = nrow(prostate[prostate$train == TRUE,])
summary.regfit.fwd$aic = summary.regfit.fwd$bic-log(n)*(c(1:nvmax)+2)+2*(c(1:nvmax)+2)
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,)

### Backward selection
nvmax=8
regfit.bwd = regsubsets(lpsa~.-train,data=subset(prostate,train==TRUE),nvmax=nvmax,method="backward")
summary.regfit.bwd <- summary(regfit.bwd)
show(summary.regfit.bwd)

# R^2 when including each of the covariates
plot(regfit.bwd,scale="r2")

#R^2 and R_adj^2 against the no. of cov. in the model
matplot(cbind(summary.regfit.bwd$rsq,summary.regfit.bwd$adjr2),type="l")

#C_p, AIC and BIC against the no. of cov. in the model
summary.regfit.bwd$aic = summary.regfit.bwd$bic-log(n)*(c(1:nvmax)+2)+2*(c(1:nvmax)+2)
matplot(cbind(summary.regfit.bwd$cp,summary.regfit.bwd$aic,summary.regfit.bwd$bic),type="l",lty=1)
p.Cp = which.min(summary.regfit.bwd$cp)
p.aic = which.min(summary.regfit.bwd$aic)
p.bic = which.min(summary.regfit.bwd$bic)
points(c(p.Cp,p.aic,p.bic),
       c(summary.regfit.bwd$cp[p.Cp],summary.regfit.bwd$aic[p.aic],summary.regfit.bwd$bic[p.bic]))
legend("topright",c("C_p","AIC","BIC"),lty=1,col=1:3,)

#Full search
system.time({regfit.full = regsubsets(lpsa~.-train,data=subset(prostate,train==TRUE))})
summary.regfit.full = summary(regfit.full)
show(summary.regfit.full$rsq)

# R^2 when including each of the covariates
plot(regfit.full,scale="r2")

#R^2 and R_adj^2 against the no. of cov. in the model
matplot(cbind(summary.regfit.full$rsq,summary.regfit.full$adjr2),type="l")

#C_p, AIC and BIC against the no. of cov. in the model
summary.regfit.full$aic = summary.regfit.full$bic-log(n)*(c(1:nvmax)+2)+2*(c(1:nvmax)+2)
matplot(cbind(summary.regfit.full$cp,summary.regfit.full$aic,summary.regfit.full$bic),type="l",lty=1)
p.Cp = which.min(summary.regfit.full$cp)
p.aic = which.min(summary.regfit.full$aic)
p.bic = which.min(summary.regfit.full$bic)
points(c(p.Cp,p.aic,p.bic),
       c(summary.regfit.full$cp[p.Cp],summary.regfit.full$aic[p.aic],summary.regfit.full$bic[p.bic]))
legend("topright",c("C_p","AIC","BIC"),lty=1,col=1:3,)

# =============================================
# Hitters data
# =============================================
#
# Summary:
# These are Major League Baseball data from the
# 1986 and 1987 seasons. This dataset was taken
# from the StatLib library. The salary data were
# originally from Sports Illustrated, April 20, 1987.
# The 1986 and career statistics were obtained from
# The 1987 Baseball Encyclopedia Update.
#
# Data description:
# The data set consists of the following 20 variables:
# Response variable:
#  Salary: annual salary on opening day in thousands of dollars
# Covariates:
#  AtBat: Number of times at bat
#  Hits: Number of hits
#  HmRun: Number of home runs
#  Runs: Number of runs
#  RBI: Number of runs batted
#  Walks: Number of walks
#  Years: Number of years in the major leagues
#  CAtBat: Number of times at bat during his career
#  CHits: Number of hits during his career
#  CHmRun: Number of home runs during his career
#  CRuns: Number of runs during his career
#  CRBI: Number of runs batted in during his career
#  CWalks: Number of walks during his career
#  League: A factor with levels A and N indicating player’s league at the end of the first year
#  Division: A factor with levels E and W indicating player’s division at the end of the first year
#  PutOuts: Number of put outs
#  Assists: Number of assists
#  Errors: Number of errors
#  NewLeague: A factor with levels A and N indicating player’s league at the beginning of the second year
#
#======================================================

# The data are available in the R package ISLR
library(ISLR)
data("Hitters")
Hitters=na.omit(Hitters) #remove missing values
class(Hitters$League)
class(Hitters$Division)
class(Hitters$NewLeague)

# Fit full model
fit.lm = lm(Salary~.,data=Hitters)
summary(fit.lm)

library(leaps)

### Forward selection
nvmax=19
regfit.fwd = regsubsets(Salary~.,data=Hitters,nvmax=nvmax,method="forward")
summary.regfit.fwd <- summary(regfit.fwd)
show(summary.regfit.fwd)

# R^2 when including each of the covariates
plot(regfit.fwd,scale="r2")

#R^2 and R_adj^2 against the no. of cov. in the model
matplot(cbind(summary.regfit.fwd$rsq,summary.regfit.fwd$adjr2),type="l")

#C_p, AIC and BIC against the no. of cov. in the model
n = nrow(Hitters)
summary.regfit.fwd$aic = summary.regfit.fwd$bic-log(n)*(c(1:nvmax)+2)+2*(c(1:nvmax)+2)
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,)

# =============================================
# Yesterday data
# =============================================
#
# Summary:
# These are simulated data, generated from an
# additive error model Y=f(x)+epsilon, with
# epsilons independent from N(0,0.01^2). The data
# set is divided into a training set with N=30,
# called yesterday, and an independent test set
# with 30 observations, called tomorrow, with the
# same covariate values
#
# Data description:
# The data set consists of the following 2 variables:
# Response variable:
#  y.yesterday: responses in the training set
#  y.tomorrow: responses in the test set
# Covariates:
#  x (same values for y.yesterday and y.tomorrow)
#
#======================================================

## Read the data
yesterday <- read.table("http://azzalini.stat.unipd.it/Book-DM/yesterday.dat", header = TRUE)
attach(yesterday)

## Plot the data
y.range <- range(c(y.yesterday, y.tomorrow))
plot(x, y.yesterday, ylim = y.range, col=2, pch = 16, ylab = "y", cex = 2)
points(x, y.tomorrow, ylim = y.range, col=4, pch = 17, ylab = "y", cex = 2)

## Fit a linear model with polynomial in x of order p, for p=1,...,23
polynomial <- function(p) lm(y.yesterday ~ poly(x,p))
polynomials <- vector("list",length=23)
for(p in 1:23)
{
  polynomials[[p]] <- polynomial(p)
}

## Plot models for p=1,5,8,23
lines(x,polynomials[[1]]$fitted.values,col=2,lwd=2)
lines(x,polynomials[[5]]$fitted.values,col=3,lwd=2)
lines(x,polynomials[[8]]$fitted.values,col=7,lwd=2)
lines(x,polynomials[[23]]$fitted.values,col=8,lwd=2)

## Compute training and test error
error.train <- error.test <- rep(0,23)
for(p in 1:23)
{
  error.train[p] <- sum((y.yesterday-polynomials[[p]]$fitted.values)^2)
  error.test[p] <- sum((y.tomorrow-polynomials[[p]]$fitted.values)^2)
}
plot(1:23,error.train,type="l",lwd=2)
lines(1:23,error.test,col=2,lwd=2)
legend(16,0.015,c("Training error","Test error"),col=1:2,lty=rep(1,2))

## Compute bias and variance
sigma.eps <- 0.01
N <- 30
source("http://www.uio.no/studier/emner/matnat/math/STK2100/v25/undervisningsmateriell/r-kode/f_true.r") #evaluation of the true f in x
polynomial.true <- function(p) lm(f.true ~ poly(x,p))
polynomials.true <- vector("list",length=23)
for(p in 1:23)
{
  polynomials.true[[p]] <- polynomial.true(p)
}
bias2 <- variance <- rep(0,23)
for(p in 1:23)
{
  bias2[p] <- mean((fitted(polynomials.true[[p]]) - f.true)^2)
  variance[p] <- sigma.eps^2*p/N
}
err <- bias2 + variance
min.err <- which.min(err)
plot(bias2, type = "b", col = 4, ylab = "Error", xlab = "Model complexity", cex = 2, pch = 16)
lines(variance, type = "b", col = 2, cex = 2, pch = 17)
lines(err, type = "b", cex = 2, pch = 18)
abline(v=min.err,lty=3,lwd=2)
legend('topright', legend = c("Bias^2", "Variance", "Total"), col = c(4, 2, 1), pch = c(16, 17, 18))

## Compute information criteria: AIC, BIC

# Refit with all the data
y.all <- c(y.yesterday,y.tomorrow)
x.all <- rep(x,2)
polynomials.all <- vector("list",length=23)
for(p in 1:23)
{
  polynomials.all[[p]] <- lm(y.all~poly(x.all,p))
}
AICvec <- BICvec <- rep(0,23)
for(p in 1:23)
{
  AICvec[p] <- AIC(polynomials.all[[p]])
  BICvec[p] <- BIC(polynomials.all[[p]])
}
min.AIC <- which.min(AICvec)
min.BIC <- which.min(BICvec)
plot(1:23, AICvec, ylim=c(min(AICvec, BICvec), max(AICvec, BICvec)), col = 2, type = "b", pch = 16, cex = 2, xlab = "Model complexity", ylab = "AIC and BIC criterion")
lines(BICvec, col = 4, type = "b", pch = 18, cex = 2)
abline(v=min.AIC,lty=3,col=2)
abline(v=min.BIC,lty=3,col=4)
legend(20,-270,legend=c("AIC","BIC"),pch=c(16,18),col=c(2,4))

## Cross-validation

# LOOCV
N <- length(y.all)
Err.loocv <- rep(0,23)
for(p in 1:23)
{
  for(i in 1:N)
  {
    polynomials.loocv <- lm(y.all~poly(x.all,p),subset=(1:N)[-i])
    Err.loocv[p] <- Err.loocv[p]+(y.all[i]-predict(polynomials.loocv,data.frame(x.all=x.all[i])))^2/N
  }
}

# K-fold CV
Err.Kcv <- rep(0,23)
set.seed(20200219) # set seed, useful for data replication
K <- 10 # set the number of folds
index <- sample(rep(1:K, ceiling(N/ K))[1:N])# split the data in approximately equisize folds
for(p in 1:23)
{
  for(k in 1:K)
  {
    polynomials.Kcv <- lm(y.all~poly(x.all,p),subset=(1:N)[-which(index == k)])
    for(i in which(index==k))
    {
      Err.Kcv[p] <- Err.Kcv[p]+(y.all[i]-predict(polynomials.Kcv,data.frame(x.all=x.all[i])))^2/N
    }
  }
}

min.loocv <- which.min(Err.loocv)
min.Kcv <- which.min(Err.Kcv)
plot(1:23, Err.loocv, ylim=c(min(Err.loocv, Err.Kcv), max(Err.loocv, Err.Kcv)), col = 2, type = "b", pch = 16, cex = 2, xlab = "Model complexity", ylab = "Leave-one-out and 10-fold CV")
lines(1:23,Err.Kcv, col = 4, type = "b", pch = 18, cex = 2)
abline(v=min.loocv,lty=3,col=2)
abline(v=min.Kcv,lty=3,col=4)
legend(19,5.8e-04,legend=c("LOOCV","10-fold CV"),pch=c(16,18),col=c(2,4))

## Bootstrap
B <- 1000
error.train.all <- rep(0,23)
for(p in 1:23)
{
  error.train.all[p] <- mean(polynomials.all[[p]]$residuals^2)
}
Err.boot.1 <- rep(0, 11) # just for the first 10 degrees of the polynomial
n.test <- 0
for (b in 1:B)
{
  set.seed(b) #set seed, useful for data replication
  index <- sample(N, replace = TRUE)
  temp.train.data <- data.frame(y=y.all[index], x=x.all[index]) #the N selected observations from the training set
  temp.test.data <- data.frame(y=y.all[-index], x=x.all[-index]) #the observations that were not selected are in the test set
  n.test <- n.test+nrow(temp.test.data) #keep track of the number of observations in the test sets
  for(p in 1:11)
  {
    polynomials.boot <- lm(y~poly(x,p),data=temp.train.data)# fit the model on the training
    Err.boot.1[p] <- Err.boot.1[p]+sum((temp.test.data$y-predict(polynomials.boot,newdata=temp.test.data))^2)
  }
}
Err.boot.1 <- Err.boot.1/n.test
min.boot.1 <- which.min(Err.boot.1)
Err.boot.632 <- 0.368*error.train.all[1:11] + 0.632*Err.boot.1 #0.632 bootstrap
min.boot.632 <- which.min(Err.boot.632)
plot(1:11, Err.boot.1, ylim=c(min(Err.boot.1,Err.boot.632),max(Err.boot.1,Err.boot.632)),type = "b", xlab = "Model complexity", col = 2, pch = 16, ylab = "Error", cex = 2, lwd = 3)
lines(1:11, Err.boot.632, type = "b", pch = 18, col=4, lwd = 3, cex=2)
abline(v=min.boot.1,lty=3,col=2,lwd=2)
abline(v=min.boot.632,lty=3,col=4,lwd=2)
legend(6,0.0011,c("Err_boot^(1)","Err_boot^(0.632)"),pch=c(16,18),col=c(2,4))