# =============================================
# Boston housing data
# =============================================
#
# Summary: 
# These are records of housing values in 506 suburbs 
# of Boston
#
# Data description: 
# The data set consists of the following 13 variables
# (the first colums contains the number of the individual):
# Response variable:
#  medv: median value of owner-occupied homes in $1000s
# Covariates:
#  crim: per capita crime rate by town
#  zn: proportion of residential land zoned for lots over 25,000 sq.ft.
#  indus: proportion of non-retail business acres per town.
#  chas: Charles River dummy variable (= 1 if tract bounds river; 0 otherwise).
#  nox: nitrogen oxides concentration (parts per 10 million).
#  rm: average number of rooms per dwelling.
#  age: proportion of owner-occupied units built prior to 1940.
#  dis: weighted mean of distances to five Boston employment centres.
#  rad: index of accessibility to radial highways.
#  tax: full-value property-tax rate per $10,000.
#  ptratio: pupil-teacher ratio by town.
#  lstat: lower status of the population (percent)
#
#======================================================

library(ISLR2)
data(Boston)

# Divide into a training and a test set
set.seed(1001)
N <- nrow(Boston)
train <- sample(1:N,round(N/2))
Boston$chas <- as.factor(Boston$chas)
Boston$rad <- as.factor(Boston$rad)

# Fit a linear regression model to all the covariates
fit.lm <- lm(medv~.,data=Boston, subset=train)
summary(fit.lm)
pred.lm <- predict(fit.lm,newdata=Boston[-train,])
err.lm <- mean((Boston$medv[-train]-pred.lm)^2)
err.lm

# Forward subset selection
library(leaps)
nvmax=19
regfit.fwd = regsubsets(medv~.,data=Boston,subset=train,nvmax=nvmax,method="forward")
summary.regfit.fwd <- summary(regfit.fwd)
show(summary.regfit.fwd)
summary.regfit.fwd$aic = summary.regfit.fwd$bic-log(length(train))*(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,ylab="")
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]))
# Selecting the subset with the AIC
rad.gr <- as.character(Boston$rad)
rad.gr[(Boston$rad == "2")|(Boston$rad == "4")|(Boston$rad == "5")|(Boston$rad == "6")|(Boston$rad == "8")] <- "1"
rad.gr <- as.factor(rad.gr)
Boston$rad.gr <- rad.gr
fit.lm.sub <- lm(medv~crim+zn+chas+nox+rm+dis+rad.gr+tax+ptratio+lstat,data=Boston,subset=train)
summary(fit.lm.sub)
anova(fit.lm.sub,fit.lm,test="Chisq")
pred.lm.sub <- predict(fit.lm.sub,newdata=Boston[-train,])
err.lm.sub <- mean((Boston$medv[-train]-pred.lm.sub)^2)
err.lm.sub
Boston <- Boston[,-14] #remove grouped rad

# Fit a ridge regression model
library(glmnet)
# Standardise the numerical covariates
num.ind <- c(1:3,5:8,10:12)
Boston.stand <- Boston
Boston.stand[,num.ind] <- scale(Boston[,num.ind])
# Make design matrix (without intercept)
X <- model.matrix(medv~.,data=Boston.stand)[,-1]
y <- Boston.stand$medv
fit.ridge.cv <- cv.glmnet(X[train,],y[train],alpha=0)
summary(fit.ridge.cv)
fit.ridge <- glmnet(X[train,],y[train],alpha=0,lambda=fit.ridge.cv$lambda.min)
pred.ridge <- predict(fit.ridge,newx=X[-train,])
err.ridge <- mean((y[-train]-pred.ridge)^2)
err.ridge

# Fit a Lasso regression model
fit.lasso.cv <- cv.glmnet(X[train,],y[train],alpha=1)
summary(fit.lasso.cv)
fit.lasso <- glmnet(X[train,],y[train],alpha=1,lambda=fit.lasso.cv$lambda.min)
pred.lasso <- predict(fit.lasso,newx=X[-train,])
err.lasso <- mean((y[-train]-pred.lasso)^2)
err.lasso

# Fit PCR (principal component regression) model 
num.ind <- c(1:3,5:8,10:12) #pick variables that are numeric
Z = prcomp(Boston[,num.ind],retx=TRUE,scale=TRUE)$x
# Fit regression model to the principal components in decreasing order 
M <- ncol(Z)
err.pca <- rep(0,M)
for(m in 1:M)
{
  Boston.pca = data.frame(Z=Z[,1:m],chas=Boston$chas,rad=Boston$rad,medv=Boston$medv)
  fit.pca = lm(medv~.,data=Boston.pca[train,])
  pred.pca = predict(fit.pca,newdata=Boston.pca[-train,])
  err.pca[m] = mean((Boston.pca$medv[-train]-pred.pca)^2)
}
plot(err.pca)

# Fit an additive model
library(gam)
fit.gam.1 <- gam(medv~s(crim)+s(zn)+s(indus)+chas+s(nox)+s(rm)+s(dis)+rad+s(tax)+s(ptratio)+s(lstat),data=Boston,subset=train)
summary(fit.gam.1)
anova(fit.lm,fit.gam.1,test="Chisq")
# Try a linear effect for ptratio
fit.gam.2 <- gam(medv~s(crim)+s(zn)+s(indus)+chas+s(nox)+s(rm)+s(dis)+rad+s(tax)+ptratio+s(lstat),data=Boston,subset=train)
summary(fit.gam.2)
anova(fit.gam.2,fit.gam.1,test="Chisq")
# Try a linear effect for tax
fit.gam.3 <- gam(medv~s(crim)+s(zn)+s(indus)+chas+s(nox)+s(rm)+s(dis)+rad+tax+ptratio+s(lstat),data=Boston,subset=train)
summary(fit.gam.3)
anova(fit.gam.3,fit.gam.2,test="Chisq")
# Try a linear effect for indus
fit.gam.4 <- gam(medv~s(crim)+s(zn)+indus+chas+s(nox)+s(rm)+s(dis)+rad+tax+ptratio+s(lstat),data=Boston,subset=train)
summary(fit.gam.4)
anova(fit.gam.4,fit.gam.3,test="Chisq")
# Try a linear effect for nox
fit.gam.5 <- gam(medv~s(crim)+s(zn)+indus+chas+nox+s(rm)+s(dis)+rad+tax+ptratio+s(lstat),data=Boston,subset=train)
summary(fit.gam.5)
anova(fit.gam.5,fit.gam.4,test="Chisq")
pred.gam <- predict(fit.gam.5,newdata=Boston[-train,])
err.gam <- mean((Boston$medv[-train]-pred.gam)^2)
err.gam
plot(fit.gam.5,se=TRUE,terms=c("s(crim)","s(zn)","s(rm)","s(dis)","s(lstat)"))

# Fit a regression tree
library(tree)
set.seed(1001)
N.tr <- length(train)
fit.t1 <- tree(medv~.,data=Boston[train,],split="deviance",control=tree.control(mincut=2,minsize=4,nobs=N.tr))
plot(fit.t1)
text(fit.t1, digits=2, pretty=3,cex=0.5)
pred.t1 = predict(fit.t1,Boston[-train,])
err.t1 <- mean((Boston$medv[-train]-pred.t1)^2)
err.t1
# Prune the tree with 10-fold cross validation
fit.p1.cv <- cv.tree(fit.t1,,FUN=prune.tree,K=10)
plot(fit.p1.cv)
best.size <- fit.p1.cv$size[which.min(fit.p1.cv$dev)]
best.size #No effect of pruning

# Fit a neural network with one hidden layer
library(nnet)
m <- 10:25
err.nnet.size <- rep(NA,length(m))
for(i in 1:length(m))
{
  print(i)
  fit.nnet.size <- nnet(medv~.,data=Boston.stand[train,],size=m[i],decay=0,linout=TRUE,MaxNWts=10000,maxit=500)
  pred <- predict(fit.nnet.size,newdata=Boston.stand[-train,])
  err.nnet.size[i] <- mean((Boston.stand$medv[-train]-pred)^2)
} 
err.nnet.size
decay <- seq(0.5,5,length=10)
m <- 20
err.nnet.decay <- rep(NA,length(decay))
for(i in 1:length(decay))
{
  print(i)
  fit.nnet.decay <- nnet(medv~.,data=Boston.stand[train,],size=m,decay=decay[i],linout=TRUE,MaxNWts=10000,maxit=500)
  pred <- predict(fit.nnet.decay,newdata=Boston.stand[-train,])
  err.nnet.decay[i] <- mean((Boston.stand$medv[-train]-pred)^2)
}  
err.nnet.decay

# Fit gradient boosted trees
library(xgboost)
x.tr <- matrix(unlist(Boston[train,1:12]),nrow=length(train),ncol=12)
y.tr <- unlist(Boston$medv[train])
x.te <- matrix(unlist(Boston[-train,1:12]),nrow=N-length(train),ncol=12)
y.te <- unlist(Boston$medv[-train])
# Use tree depth 1 and 10-fold cross-validation to select M and the learning rate is set to 0.1
fit.boost.1.cv <- xgb.cv(data=xgb.DMatrix(data=x.tr,label=y.tr),nrounds = 1000,eta=0.1,nfold = 10,metrics=list("rmse"),max_depth=1,objective="reg:squarederror")
best.boost.1 <- which.min(fit.boost.1.cv$evaluation_log$test_rmse_mean)
fit.boost.1 <- xgboost(data=x.tr,label=y.tr,max.depth=1,nrounds=best.boost.1,eta=0.1,objective="reg:squarederror")
pred.boost.1 <- predict(fit.boost.1,newdata=x.te)
err.boost.1 <- mean((y.te - pred.boost.1)^2)
err.boost.1
# Use tree depth 2
fit.boost.2.cv <- xgb.cv(data=xgb.DMatrix(data=x.tr,label=y.tr),nrounds = 1000,eta=0.1,nfold = 10,metrics=list("rmse"),max_depth=2,objective="reg:squarederror")
best.boost.2 <- which.min(fit.boost.2.cv$evaluation_log$test_rmse_mean)
fit.boost.2 <- xgboost(data=x.tr,label=y.tr,max.depth=2,nrounds=best.boost.2,eta=0.1,objective="reg:squarederror")
pred.boost.2 <- predict(fit.boost.2,newdata=x.te)
err.boost.2 <- mean((y.te - pred.boost.2)^2)
err.boost.2
# Use tree depth 3
fit.boost.3.cv <- xgb.cv(data=xgb.DMatrix(data=x.tr,label=y.tr),nrounds = 1000,eta=0.1,nfold = 10,metrics=list("rmse"),max_depth=3,objective="reg:squarederror")
best.boost.3 <- which.min(fit.boost.3.cv$evaluation_log$test_rmse_mean)
fit.boost.3 <- xgboost(data=x.tr,label=y.tr,max.depth=3,nrounds=best.boost.3,eta=0.1,objective="reg:squarederror")
pred.boost.3 <- predict(fit.boost.3,newdata=x.te)
err.boost.3 <- mean((y.te - pred.boost.3)^2)
err.boost.3

# Fit random forests
library(randomForest)
x.tr <- matrix(unlist(Boston[train,1:12]),nrow=length(train),ncol=12)
y.tr <- unlist(Boston$medv[train])
x.te <- matrix(unlist(Boston[-train,1:12]),nrow=N-length(train),ncol=12)
y.te <- unlist(Boston$medv[-train])
# Fit a random forest with m=4 (the default value) and B=2500
m <- 4
B <- 2500
fit.rf.1 <- randomForest(x=x.tr,y=y.tr,ntree=B,mtry=m)
pred.rf.1 <- predict(fit.rf.1,newdata=x.te) 
err.rf.1 <- mean((y.te - pred.rf.1)^2)
err.rf.1
plot(fit.rf.1$mse,type="l",ylim=c(13,25),lwd=2,xlab="Number of trees",ylab="OOB error estimate")
# Random forest with m=6 
m <- 6
fit.rf.2 <- randomForest(x=x.tr,y=y.tr,ntree=B,mtry=m)
pred.rf.2 <- predict(fit.rf.2,newdata=x.te) 
err.rf.2 <- mean((y.te - pred.rf.2)^2)
err.rf.2
lines(fit.rf.2$mse,lwd=2,col=2)
# Random forest with m=6 
m <- 8
fit.rf.3 <- randomForest(x=x.tr,y=y.tr,ntree=B,mtry=m)
pred.rf.3 <- predict(fit.rf.3,newdata=x.te) 
err.rf.3 <- mean((y.te - pred.rf.3)^2)
err.rf.3
lines(fit.rf.3$mse,lwd=2,col=3)