#####################################################################################################
##
##                                Classification and Regression Tree
##                                    Bagging and Random Forests
##                                              Boosting
##
## E. Flachaire, 13/12/2021
##
#####################################################################################################

rm(list=ls())

# Classification Tree -----------------------------------------------------

library(PASWR)      # Data contained in that package
data(titanic3)

library(rpart)      # CART package
library(rpart.plot) # fancy plots

tree <- rpart(survived~sex+age+pclass+sibsp+parch,
              data=titanic3, method = "class")
prp(tree, extra=1, faclen=5,
    box.col = c("indianred1", "aquamarine")[tree$frame$yval])


# Variable importance
tree$variable.importance
barplot(tree$variable.importance, horiz=TRUE, col="yellow3")


# Regression Tree ---------------------------------------------------------

library(ISLR)
# remove observations that are missing Salary values
df <- Hitters[complete.cases(Hitters$Salary),]

library(rpart)      # load CART library
library(rpart.plot) # fancy plots

# Estimate the tree
tree <- rpart(log(Salary)~Years+Hits+AtBat, data=df, cp=0)
# plot the tree
prp(tree, extra=1, faclen=5)


# Tree pruning
tree <- rpart(log(Salary)~Years+Hits+AtBat, data=df) # CV
prp(tree, extra=1, faclen=5, main = "Pruned tree")

# Bagging -----------------------------------------------------------------

library(rpart)      # load CART library
library(rpart.plot) # fancy plots
n <- nrow(titanic3)
par(mfrow=c(3,3))
for(i in 1:9){
  idx <- sample(1:n, size=n, replace=TRUE)
  cart <- rpart(as.factor(survived)~sex+age+pclass+sibsp+parch,
                data=titanic3[idx,], cp=0) # unpruned
  prp(cart, type=1, extra=1)
}
par(mfrow=c(1,1))

# Random Forest -----------------------------------------------------------

library(rpart)        # load CART library
library(rpart.plot)   # fancy plots
library(randomForest) # Random Forest

set.seed(1)
n <- 200
x <- runif(n)
y <- sin(12*(x+.2))/(x+.2) + rnorm(n)/2

fit.tr <- rpart(y~x) # CART
fit.ba1 <- randomForest(y~x) # no depth control
fit.ba2 <- randomForest(y~x, maxnodes=20) #depth control

u <- seq(min(x), max(x), length.out=1000)
plot(x,y,col="gray", main="n=200")
lines(u, predict(fit.ba1, data.frame(x=u)), col="green")
lines(u, predict(fit.ba2, data.frame(x=u)), col="red")
lines(u, predict(fit.tr, data.frame(x=u)), col="blue")


set.seed(1)
n <- 1000
x <- runif(n)
y <- sin(12*(x+.2))/(x+.2) + rnorm(n)/2

fit.tr <- rpart(y~x) # CART
fit.ba1 <- randomForest(y~x) # no depth control
fit.ba2 <- randomForest(y~x, maxnodes=20) #depth control

u <- seq(min(x), max(x), length.out=1000)
plot(x,y,col="gray", main="n=1000")
lines(u, predict(fit.ba1, data.frame(x=u)), col="green")
lines(u, predict(fit.ba2, data.frame(x=u)), col="red")
lines(u, predict(fit.tr, data.frame(x=u)), col="blue")


# Tuning parameters
randomForest(y~x)$mse[500] # OOB-MSE, no depth control
maxnode <- c(10,50,100, 500, 1000, 2000)
for(i in 1:length(maxnode)){
  # OOB-MSE depth control
  aa <- randomForest(y~x, maxnodes=maxnode[i])$mse[500]
  print(c(maxnode[i], aa))
}

# Variable importance
rf <- randomForest(as.factor(survived)~sex+age+pclass+sibsp+parch,
                   data = titanic3, na.action=na.omit)
importance(rf)


# Boosting ----------------------------------------------------------------


set.seed(1)
n <- 200
x <- runif(n)
y <- sin(12*(x+.2))/(x+.2) + rnorm(n)/2

library(rpart)
library(gbm)
nb <- 500
# By default: interaction.depth=1 and shrinkage=0.1
fit.bo <- gbm(y~x, distribution="gaussian", n.tree=nb)
fit.tre <- rpart(y~x)

u <- seq(min(x), max(x), length.out=1000)
plot(x,y, col="gray", main="n=200", xlab=NA, ylab=NA)
lines(u, predict(fit.tre, data.frame(x=u)), col="blue")
lines(u, predict(fit.bo, data.frame(x=u), n.tree=nb), col="purple")


set.seed(1)
n <- 1000
x <- runif(n)
y <- sin(12*(x+.2))/(x+.2) + rnorm(n)/2


nb <- 500
# By default: interaction.depth=1 and shrinkage=0.1
fit.bo <- gbm(y~x, distribution="gaussian", n.tree=nb)
fit.tre <- rpart(y~x)

u <- seq(min(x), max(x), length.out=1000)
plot(x,y, col="gray", main="n=1000", xlab=NA, ylab=NA)
lines(u, predict(fit.tre, data.frame(x=u)), col="blue")
lines(u, predict(fit.bo, data.frame(x=u), n.tree=nb), col="purple")