# Machine learning for econometricians
# Course for the ECB
# by Ewen Gallic
# December 2021
# Hands-on Session 2
# Trees


# 1.1 Data Used in the Notebook -------------------------------------------



library(tidyverse)

url <- "https://egallic.fr/Enseignement/ML/ECB/data/SeoulBikeData.csv"
bike <- read_csv(url, locale = locale(encoding = "latin1"))


table(bike$`Functioning Day`)

# Removing some observations
bike <- 
  bike %>% 
  filter(`Functioning Day` == "Yes") %>% 
  select(-`Functioning Day`)

# Renaming
bike <- 
  bike %>% 
  rename(
    "date" = `Date`,
    "rented_bike_count" = `Rented Bike Count`,
    "hour" = `Hour`,
    "temperature" = `Temperature(°C)`,
    "humidity" = `Humidity(%)`,
    "wind_speed" = `Wind speed (m/s)`,
    "visibility" = `Visibility (10m)`,
    "dew_point_temperature" = `Dew point temperature(°C)`,
    "solar_radiation" = `Solar Radiation (MJ/m2)`,
    "rainfall" = `Rainfall(mm)`,
    "snowfall" = `Snowfall (cm)`,
    "seasons" = `Seasons`,
    "holiday" = `Holiday`
  )


loc_time_english <- 
  ifelse(.Platform$OS.type == "unix", "en_US", "english_us")


library(lubridate)

# Creating new variables
bike <- 
  bike %>% 
  mutate(
    date = dmy(date),
    year = year(date),
    month = month(date, label = TRUE, locale = loc_time_english),
    month = factor(as.character(month),
                   levels = c("Jan", "Feb", "Mar", "Apr", 
                              "May", "Jun", "Jul", "Aug", 
                              "Sep", "Oct", "Nov", "Dec")),
    week_day = wday(date, label = TRUE, locale = loc_time_english),
    week_day = factor(as.character(week_day),
                      levels = c("Mon", "Tue", "Wed", 
                                 "Thu", "Fri", "Sat", "Sun")),
    seasons = factor(seasons, levels = c("Spring", "Summer",
                                         "Autumn", "Winter"))
  )

# Target variable
summary(bike$rented_bike_count)


ggplot(data = bike, mapping = aes(x = rented_bike_count)) +
  geom_histogram(colour = "black") +
  labs(x = "Rented bike count", y = "Count",
       title = "Distribution of rented bike count.") +
  theme(plot.title.position = "plot")


ggplot(data = bike, mapping = aes(x = rented_bike_count)) +
  geom_histogram(colour = "black") +
  labs(x = "Rented bike count", y = "Count", title = "Distribution of rented bike count by month.") +
  facet_wrap(~month) +
  theme(plot.title.position = "plot")


ggplot(data = bike, mapping = aes(x = rented_bike_count)) +
  geom_histogram(colour = "black") +
  labs(x = "Rented bike count", y = "Count", title = "Distribution of rented bike count by weekday.") +
  facet_wrap(~week_day) +
  theme(plot.title.position = "plot")


# Summary statistics
library(arsenal)
tableby_control <- tableby.control(
  numeric.stats=c("Nmiss", "meansd", "median", "q1q3")
)
tab <- tableby(week_day~rented_bike_count, data = bike, 
               control = tableby_control)

summary(tab, text = TRUE)


bike %>% 
  group_by(hour) %>% 
  summarise(rented_bike_count = mean(rented_bike_count)) %>% 
  ggplot(data = ., aes(x = hour, y = rented_bike_count)) +
  geom_bar(stat = "identity") +
  coord_polar() +
  labs(x = NULL, y = "Rented bike count",
       title = "Rented bike count per hour") +
  scale_x_continuous(breaks = seq(0,24, by = 2)) +
  theme(plot.title.position = "plot")



bike %>% 
  group_by(hour, seasons) %>% 
  summarise(rented_bike_count = mean(rented_bike_count)) %>% 
  ggplot(data = ., aes(x = hour, y = rented_bike_count)) +
  geom_bar(stat = "identity", mapping = aes(fill = seasons)) +
  coord_polar() +
  labs(x = NULL, y = "Rented bike count",
       title = "Rented bike count per hour and per season.") +
  scale_x_continuous(breaks = seq(0,24, by = 2)) +
  facet_wrap(~seasons) +
  scale_fill_manual("Seasons", 
                    values = c("Spring" = "#004D40", "Summer" = "#FFC107",
                               "Autumn" = "#D81B60", "Winter" = "#1E88E5")) +
  theme(plot.title.position = "plot")

# Weather variiables
library(cowplot)
p_temp <- 
  ggplot(data = bike, 
         mapping = aes(x = temperature, y = rented_bike_count)) +
  geom_point(alpha = .3) +
  geom_smooth() +
  labs(x = "Temperature (°C)", y = NULL,
       title = "Rented count vs temperature.") +
  theme(plot.title.position = "plot")

p_humidity <- 
  ggplot(data = bike, 
         mapping = aes(x = humidity, y = rented_bike_count)) +
  geom_point(alpha = .3) +
  geom_smooth() +
  labs(x = "Humidity (%)", y = NULL,
       title = "Rented count vs humidity")

p_wind <- 
  ggplot(data = bike, 
         mapping = aes(x = wind_speed, y = rented_bike_count)) +
  geom_point(alpha = .3) +
  geom_smooth() +
  labs(x = "Wind speed (m/s)", y = NULL,
       title = "Rented count vs wind speed")

p_visibility <- 
  ggplot(data = bike,
         mapping = aes(x = visibility, y = rented_bike_count)) +
  geom_point(alpha = .3) +
  geom_smooth() +
  labs(x = "Visibility (10m)", y = NULL,
       title = "Rented count vs visibility")

p_dew_point <- 
  ggplot(data = bike,
         mapping = aes(x = dew_point_temperature, y = rented_bike_count)) +
  geom_point(alpha = .3) +
  geom_smooth() +
  labs(x = "Dew point temperature (°C)", y = NULL,
       title = "Rented count vs dew point temperature")

p_solar_radiation <- 
  ggplot(data = bike, 
         mapping = aes(x = solar_radiation, y = rented_bike_count)) +
  geom_point(alpha = .3) +
  geom_smooth() +
  labs(x = "Solar Radiation (MJ/m2)", y = NULL,
       title = "Rented count vs solar radiation")


p_rainfall <- 
  ggplot(data = bike, 
         mapping = aes(x = rainfall, y = rented_bike_count)) +
  geom_point(alpha = .3) +
  geom_smooth() +
  labs(x = "Rainfall (mm)", y = NULL,
       title = "Rented count vs rainfall")

p_snowfall <- 
  ggplot(data = bike, 
         mapping = aes(x = snowfall, y = rented_bike_count)) +
  geom_point(alpha = .3) +
  geom_smooth() +
  labs(x = "Snowfall (cm)", y = NULL,
       title = "Rented count vs snowfall")


# Plots on a grid
cowplot::plot_grid(p_temp, p_humidity, p_wind,
                   p_visibility, p_dew_point, p_rainfall,
                   p_snowfall, ncol = 2)


# Correlation matrix
correl_mat <- 
  bike %>% 
  select_if(is.numeric) %>% 
  cor() %>% 
  round(2)
correl_mat


library(corrplot)
corrplot.mixed(correl_mat, order = 'AOE', main = "Correlation plot.")


#
tableby_control <- tableby.control(
  numeric.stats=c("Nmiss", "meansd", "range", "median", "q1q3")
)
tab <- tableby(~., data = bike, control = tableby_control)
summary(tab, text = NULL) %>% 
  kableExtra::kable(caption = "Summary statistics.", booktabs = T,
                    format = "latex", longtable = TRUE) %>%
  kableExtra::kable_classic(full_width = F, html_font = "Cambria") %>% 
  kableExtra::kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    font_size = 6)


# Defining a binary outcome variable
bike <- 
  bike %>% 
  mutate(y_binary = ifelse(rented_bike_count < 300, "Low", "High"))


table(bike$y_binary)
prop.table(table(bike$y_binary))


tableby_control <- tableby.control(
  numeric.stats=c("Nmiss", "meansd", "range", "median", "q1q3")
)
tab <- tableby(y_binary~., data = bike, control = tableby_control)
summary(tab, text = NULL) %>% 
  kableExtra::kable(
    caption = "Summary statistics depending on the binary response variable.", 
    booktabs = T,
    format = "latex", longtable = TRUE) %>%
  kableExtra::kable_classic(full_width = F, html_font = "Cambria") %>% 
  kableExtra::kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    font_size = 6)



# 1.2 Training and Validation Sets ----------------------------------------

# Sort by ascending values of date and then by ascending values of hour
bike <- 
  bike %>% 
  arrange(date, hour)


n_train <- round(.8*nrow(bike))

# Training set
df_train <- 
  bike %>% 
  slice(1:n_train)

# Validation set
df_valid <- 
  bike %>% 
  slice(-(1:n_train))


dim(df_train)
dim(df_valid)


# 1.3 Decision Trees ------------------------------------------------------

# Libraries used to create and plot decision trees
library(rpart)
library(rpart.plot)


# 1.3.1 Regression Trees --------------------------------------------------



part_tree <- 
  rpart(rented_bike_count ~., 
        data = df_train %>% select(-y_binary, -date),
        method = "anova")


part_tree

rpart.plot(part_tree)


mean(df_train$rented_bike_count)

# Left node
df_train %>% 
  filter(temperature < 12.15) %>% 
  summarise(mean = mean(rented_bike_count))

# Left of the left node
df_train %>% 
  filter(temperature < 12.15) %>% 
  filter(temperature < 5.65) %>% 
  count() %>% 
  mutate(prop = n / nrow(df_train))



df_train %>% 
  filter(temperature < 12.15) %>% 
  filter(temperature < 5.65) %>% 
  summarise(pred = mean(rented_bike_count))

# Right of the left node
df_train %>% 
  filter(temperature < 12.15) %>% 
  filter(temperature > 5.65) %>% 
  count() %>% 
  mutate(prop = n / nrow(df_train))


df_train %>% 
  filter(temperature < 12.15) %>% 
  filter(temperature > 5.65) %>% 
  summarise(pred = mean(rented_bike_count))


compute_rss <- function(observed, predicted){
  sum((observed-predicted)^2)
}


rss_init <- 
  compute_rss(df_train$rented_bike_count, 
              mean(df_train$rented_bike_count))
rss_init



variable_split <- "temperature"
threshold <- 12.15


tmp_data <- 
  df_train %>% 
  mutate(below_t = temperature < 12.15) %>% 
  select(rented_bike_count, below_t) %>% 
  group_by(below_t) %>% 
  mutate(pred = mean(rented_bike_count))

tmp_data


rss_after_split <- 
  compute_rss(tmp_data$rented_bike_count, 
              tmp_data$pred)
rss_after_split


# improvement due to the split, in terms of percentage deviation of the RSS
-(rss_after_split - rss_init)/rss_init


head(part_tree$splits)


#' @param data data frame
#' @param variable_to_predict name of the variable to predict
#' @param variable_split name of the variable used to create the partitions
#' @param threshold threshold used to create the two partitions
rss_split_numeric <- 
  function(data, variable_to_predict, 
           variable_split, threshold){
    tmp_data <- 
      data %>% 
      mutate(below_t = !!sym(variable_split) < threshold) %>% 
      select(!!sym(variable_to_predict), below_t) %>% 
      group_by(below_t) %>% 
      mutate(pred = mean(!!sym(variable_to_predict)))
    
    compute_rss(tmp_data[[variable_to_predict]], tmp_data[["pred"]])
  }


rss_split_numeric(
  data = df_train, variable_to_predict = "rented_bike_count",
  variable_split = "temperature", threshold = 12.15)


# Looping over many cutoff values

number_of_cuts <- 1000
thresholds <- 
  seq(min(df_train$temperature), max(df_train$temperature),
      length.out = number_of_cuts)



rss_tmp <- rep(NA, number_of_cuts)
for(i in 1:number_of_cuts){
  rss_tmp[i] <- 
    rss_split_numeric(
      data = df_train, variable_to_predict = "rented_bike_count",
      variable_split = "temperature", threshold = thresholds[i])
}


# Best cutoff
thresholds[which.min(rss_tmp)]


df_train %>% 
  mutate(below_t = temperature < thresholds[which.min(rss_tmp)]) %>% 
  select(rented_bike_count, below_t) %>% 
  group_by(below_t) %>% 
  summarise(pred = mean(rented_bike_count),
            n = n()) %>% 
  ungroup() %>% 
  mutate(prop = round(100*n/sum(n)))


# Now, turning to splits made with a categorical variable

#' @param classes vector of names of classes
possible_splits <- function(classes){
  l <- length(classes)
  if(l>2){
    resul <- resul <- rep(NA, 2^(l-1) - 1)
    i <- 1
    for(k in 1:(l-1)){
      tmp <- combn(c(classes), k)
      for(j in 1:ncol(tmp)){
        if(tmp[1, j] != classes[1]) break
        current <- str_c(tmp[,j], collapse = ",")
        resul[i] <- str_c(tmp[,j], collapse = ",")
        i <- i+1
      }
    }
  }else{
    # If only two classes
    resul <- classes[1]
  }
  resul
}


possible_splits(c("A", "B", "C", "D"))


rpart.plot(part_tree, main = "The same first decision tree.")


df_tmp <- 
  df_train %>% 
  filter(temperature >= 12.15 & hour < 15.5 & solar_radiation >= 0.205)

nrow(df_tmp)

# The different levels (or classes) of the month variable
classes <- as.character(unique(df_tmp$month))
classes

# All possible splits with this variable
list_possible_splits <- possible_splits(classes)
head(list_possible_splits)


length(list_possible_splits)

list_possible_splits[2]

str_split(list_possible_splits[2], ",")[[1]]

tmp_data <- 
  df_tmp %>% 
  mutate(var_in_categ = month %in% 
           str_split(list_possible_splits[2],",")[[1]])%>% 
  select(rented_bike_count, var_in_categ) %>% 
  group_by(var_in_categ) %>% 
  mutate(pred = mean(rented_bike_count))


table(tmp_data$var_in_categ)


compute_rss(tmp_data[["rented_bike_count"]], tmp_data[["pred"]])

# Let us wrap this up in a function

#' @param data data frame
#' @param variable_to_predict name of the variable to predict
#' @param variable_split name of the variable used to create the partitions
#' @param split_rule vector of classes used to create the two partitions
rss_split_categ <- 
  function(data, variable_to_predict, 
           variable_split, split_rule){
    tmp_data <- 
      data %>% 
      mutate(var_in_categ = !!sym(variable_split) %in% split_rule) %>% 
      select(!!sym(variable_to_predict), var_in_categ) %>% 
      group_by(var_in_categ) %>% 
      mutate(pred = mean(!!sym(variable_to_predict)))
    
    compute_rss(tmp_data[[variable_to_predict]], tmp_data[["pred"]])
  }


rss_split_categ(df_tmp, 
                variable_to_predict = "rented_bike_count",
                variable_split = "month", 
                split_rule = str_split(list_possible_splits[2], ",")[[1]])


# loop over all the possible splits, and store the RSS at each iteration.

rss_tmp <- rep(NA, length(list_possible_splits))
for(i in 1:length(list_possible_splits)){
  rss_tmp[i] <- 
    rss_split_categ(
      df_tmp, 
      variable_to_predict = "rented_bike_count",
      variable_split = "month", 
      split_rule = str_split(list_possible_splits[i], ",")[[1]])
}

# Best split
list_possible_splits[which.min(rss_tmp)]



df_tmp %>% 
  mutate(var_in_categ = month %in% 
           str_split(list_possible_splits[which.min(rss_tmp)],",")[[1]])%>% 
  select(rented_bike_count, var_in_categ) %>% 
  group_by(var_in_categ) %>% 
  summarise(pred = mean(rented_bike_count),
            n = n()) %>% 
  mutate(prop = round(100*n/nrow(df_train)))

rpart.plot(part_tree, main = "The same (again) first decision tree.")


# 1.3.2 Stopping the Recursive Splitting Process --------------------------

# Decision tree where splits are made if there is at least 1000 obs. 
# in the node and if the number of obs. in the resulting leaves are at least 500.

part_tree <- 
  rpart(rented_bike_count ~., 
        data = df_train %>% select(-y_binary, -date),
        method = "anova",
        minsplit = 1000,
        minbucket = 500,
        cp=0)
rpart.plot(part_tree)


# Growing the tree stopped earlier as we imposed restrictions on the improvement needed to make split.

part_tree <- 
  rpart(rented_bike_count ~., 
        data = df_train %>% select(-y_binary, -date),
        method = "anova",
        minsplit = 1000,
        minbucket = 500,
        cp=.01)
rpart.plot(part_tree)


# Decision tree when constraining its depth.
part_tree <- 
  rpart(rented_bike_count ~., 
        data = df_train %>% select(-y_binary, -date),
        method = "anova",
        maxdepth = 3,
        cp = 0)
rpart.plot(part_tree)

# Toying with dummy data

set.seed(123)
n <- 100
x <- runif(n=n, min=0, max=10)
eps <- rnorm(n, 0, 4)
f <- function(x) .1*x^3-1*x^2+x
y <- f(x)+eps


df_sim <- tibble(x=x, y=y)
ggplot(data = df_sim, mapping = aes(x = x, y = y)) +
  geom_line(data = tibble(x=seq(0, 10, by = .1), y=f(x)),
            size = 1.1, linetype = "dashed", colour = "#AA2F2F") +
  geom_point() +
  labs(title = "Generating Data Process and generated data.") +
  theme(plot.title.position = "plot")


# Decision tree built on the synthetic data with a maximum depth of 1.
part_tree <- 
  rpart(y ~ x, 
        data = df_sim,
        method = "anova",
        maxdepth = 1,
        cp = 0
  )
rpart.plot(part_tree)


ggplot(mapping = aes(x = x, y = y)) +
  geom_line(data = tibble(x=seq(0, 10, by = .1), y=f(x)),
            size = 1.1, linetype = "dashed", colour = "#AA2F2F") +
  geom_point(data = df_sim) +
  geom_segment(data = tibble(x = c(-Inf, 8.7), 
                             xend = c(8.7, Inf),
                             y = c(-4.9, 4.3), 
                             yend = c(-4.9, 4.3)),
               mapping = aes(x=x, y=y, xend=xend, yend=yend),
               colour = "#005A8B", size = 1.5) +
  labs(title = "Decision boundary of the grown tree.") +
  theme(plot.title.position = "plot")


# function that will return the decision boundaries
get_pred_part <- function(max_depth, min_split){
  part_tree <- 
    rpart(y ~ x, 
          data = df_sim,
          method = "anova",
          minsplit = min_split,
          maxdepth = max_depth,
          cp = 0
    )
  
  tibble(
    x = seq(0, 10, by = .1),
    y = predict(part_tree, newdata = tibble(x=seq(0, 10, by = .1))),
    max_depth = max_depth,
    min_split = min_split
  )
  
}

# looping over all the combinations that can be made up with the values for maxdepth and minsplit that we define, 
# and apply the get_pred_part() at each iteration.

predicted_vals <- 
  list(
    max_depth = c(1, 3, 5, 10),
    min_split = c(2, 10, 20)) %>% 
  cross() %>% 
  map_df(purrr::lift(get_pred_part))
predicted_vals


# Varying the parameters affect the decision boundary and may lead to overfitting.
library(ggtext)
ggplot(mapping = aes(x = x, y = y)) +
  geom_line(data = tibble(x=seq(0, 10, by = .1), y=f(x)),
            size = 1.1, linetype = "dashed", colour = "#AA2F2F") +
  geom_point(data = df_sim) +
  geom_step(data = predicted_vals %>% 
              mutate(max_depth = factor(
                max_depth, 
                levels = sort(unique(predicted_vals$max_depth)),
                labels = str_c("Max depth: ", 
                               sort(unique(predicted_vals$max_depth)))),
                min_split = factor(
                  min_split, 
                  levels = sort(unique(predicted_vals$min_split), 
                                decreasing = TRUE),
                  labels = str_c("Min node size: ", 
                                 sort(unique(predicted_vals$min_split), 
                                      decreasing = TRUE)))
              ),
            colour = "#005A8B", size = 1.5) +
  facet_grid(max_depth~min_split) +
  labs(title = str_c(
    "Change in **<span style='color:#005A8B'>",
    "decision boundaries</span>** induced by parameter variations")
  ) +
  theme(
    plot.title = element_markdown(lineheight = 1.1),
    plot.title.position = "plot",
    legend.text = element_markdown(size = 11)
  )


# 1.3.2.1 Pruning ---------------------------------------------------------

set.seed(123)
large_tree <- 
  rpart(rented_bike_count ~., 
        data = df_train %>% select(-y_binary, -date),
        method = "anova",
        control = rpart.control(xval = 10, cp = -1, minbucket = 20))

# Number of nodes
nodes_large_tree <- as.numeric(rownames(large_tree$frame))
max(rpart:::tree.depth(nodes_large_tree))


# Relative error depending on the complexity parameter (10-fold cross-validation results),
# for the Seoul bike data.
plotcp(large_tree)


head(large_tree$cptable)


# finding the minimum value of the cost-parameter
min_val <- 
  as_tibble(large_tree$cptable) %>%
  arrange(xerror) %>% 
  slice(1)

min_val

# all the values of the cost-parameter within a 1-standard error deviation from that value
candidates <- 
  as_tibble(large_tree$cptable) %>%
  filter( (xerror > min_val$xerror-min_val$xstd) & 
            (xerror < min_val$xerror+min_val$xstd)) %>% 
  arrange(desc(CP))
candidates


# Among those candidates, we can pick up the values that results in the smallest subtree:
cp_val <- 
  candidates %>% 
  slice(1) %>% 
  magrittr::extract2("CP")
cp_val


# tree can be pruned accordingly with this cost-complexity value
pruned_tree <- prune(large_tree, cp=cp_val)


nodes_pruned_tree <- as.numeric(rownames(pruned_tree$frame))
max(rpart:::tree.depth(nodes_pruned_tree))



small_tree <- 
  rpart(rented_bike_count ~., 
        data = df_train %>% select(-y_binary, -date),
        method = "anova",
        maxdepth = 2,
        cp = 0,
        xval = 0)

# Number of nodes
nodes_small_tree <- as.numeric(rownames(small_tree$frame))
max(rpart:::tree.depth(nodes_small_tree))


# Let us compute the MSE on the validation set:
pred_valid_small <- predict(small_tree, newdata = df_valid)
pred_valid_large <- predict(large_tree, newdata = df_valid)
pred_valid_pruned <- predict(pruned_tree, newdata = df_valid)

compute_mse <- function(observed, predicted){
  mean((observed-predicted)^2)
}

mse_small_tree <- 
  compute_mse(observed = df_valid$rented_bike_count,
              predicted = pred_valid_small)
mse_large_tree <- 
  compute_mse(observed = df_valid$rented_bike_count,
              predicted = pred_valid_large)
mse_pruned_tree <- 
  compute_mse(observed = df_valid$rented_bike_count,
              predicted = pred_valid_pruned)

mse <- scales::number(c(mse_small_tree, mse_large_tree,
                        mse_pruned_tree))
names(mse) <- c("Small", "Large", "Pruned")
mse



set.seed(234)
x_new <- runif(n=n, min=0, max=10)
eps <- rnorm(n, 0, 4)
y_new <- f(x)+eps
df_sim_new <- tibble(x=x_new, y=y_new)
ggplot(mapping = aes(x = x, y = y)) +
  geom_line(data = tibble(x=seq(0, 10, by = .1), y=f(x)),
            colour = "black", size = 1.1, linetype = "dashed") +
  geom_point(data = df_sim %>% mutate(sample = "Train") %>% 
               bind_rows(
                 df_sim_new %>% mutate(sample = "Test")
               ),
             mapping = aes(colour = sample)) +
  scale_colour_manual("Sample", 
                      values = c("Train" = "#D55E00", "Test" = "#009E73")) +
  labs(title = "Synthetic data.") +
  theme(plot.title.position = "plot")



part_tree_large <- 
  rpart(y ~ x, 
        data = df_sim,
        method = "anova",
        control = rpart.control(
          xval = 10, 
          cp = -1, 
          minbucket = 5)
  )
rpart.plot(part_tree_large, main = "Unpruned tree.")


# Relative error depending on the complexity parameter\nfor the synthetic data.
plotcp(part_tree_large)


# select the value that minimises the error:
min_val <- 
  as_tibble(part_tree_large$cptable) %>%
  arrange(xerror) %>% 
  slice(1)
min_val


# consider all values within a 1-standard error:
candidates <- 
  as_tibble(part_tree_large$cptable) %>%
  filter( (xerror > min_val$xerror-min_val$xstd) & 
            (xerror < min_val$xerror+min_val$xstd)) %>% 
  arrange(desc(CP))
candidates


# pick the values among the candidates that produces the shallowest tree:
cp_val <- 
  candidates %>% 
  slice(1) %>% 
  magrittr::extract2("CP")
cp_val



pruned_tree <- prune(part_tree_large, cp=cp_val)
rpart.plot(pruned_tree, main = "Pruned tree, synthetic data.")



nodes_pruned_tree <- as.numeric(rownames(pruned_tree$frame))
max(rpart:::tree.depth(nodes_pruned_tree))


# visualise the boundaries of each tree:

boundaries <- 
  tibble(
    x = seq(0, 10, by = .1),
    y = predict(part_tree_large,
                newdata = tibble(x=seq(0, 10, by = .1)))
  )
boundaries_pruned <- 
  tibble(
    x = seq(0, 10, by = .1),
    y = predict(pruned_tree,
                newdata = tibble(x=seq(0, 10, by = .1)))
  )

ggplot(mapping = aes(x = x, y = y)) +
  geom_line(data = tibble(x=seq(0, 10, by = .1), y=f(x)),
            size = 1.1, linetype = "dashed", colour = "#AA2F2F") +
  geom_point(data = df_sim %>% mutate(sample = "Train") %>% 
               bind_rows(
                 df_sim_new %>% mutate(sample = "Test")
               ),
             mapping = aes(shape = sample)) +
  geom_step(data = boundaries %>% mutate(model = "Large Tree") %>% 
              bind_rows(
                boundaries_pruned %>% mutate(model = "Pruned Tree")
              ),
            mapping = aes(colour = model),
            # colour = "#005A8B",
            size = 1.5) +
  
  labs(title = str_c(
    "Decision boundaries for the **<span style='color:#E69F00'>",
    "large</span>** and for the **<span style='color:#56B4E9'>",
    "pruned</span>** trees")
  ) +
  scale_colour_manual("Boundary", 
                      values = c("Large Tree" = "#E69F00", 
                                 "Pruned Tree" = "#56B4E9")) +
  scale_shape_discrete("Dataset") +
  theme(
    plot.title = element_markdown(lineheight = 1.1),
    plot.title.position = "plot",
    legend.text = element_markdown(size = 11)
  )


# using these boundaries, let us compute the MSE on the unseen data (test set):
pred_valid_large <- predict(part_tree_large, newdata = df_sim_new)
pred_valid_pruned <- predict(pruned_tree, newdata = df_sim_new)

compute_mse <- function(observed, predicted){
  mean((observed-predicted)^2)
}


mse_large_tree <- 
  compute_mse(observed = df_sim_new$y, predicted = pred_valid_large)
mse_pruned_tree <- 
  compute_mse(observed = df_sim_new$y, predicted = pred_valid_pruned)

mse <- scales::number(c(mse_large_tree, mse_pruned_tree))
names(mse) <- c("Large", "Pruned")
print("MSE : \n")
mse


# 1.3.3 Classification Trees ----------------------------------------------


classif_tree_gini <- 
  rpart(y_binary ~., 
        data = df_train %>% select(-rented_bike_count, -date),
        method = "class",
        minsplit = 20,
        minbucket = 7,
        parms = list(split = "gini"))



rpart.plot(classif_tree_gini, extra = 106, main = "A first classification tree grown on Seoul bike data.")


prop.table(table(df_train$y_binary))



df_train %>% 
  filter(temperature >= 3.35) %>%
  group_by(y_binary) %>% 
  count() %>% 
  ungroup() %>% 
  mutate(prop = round(n/sum(n), 2))



df_train %>% 
  group_by(temperature < 3.35) %>% 
  count() %>% 
  ungroup() %>% 
  mutate(prop = n / sum(n))



df_train %>% 
  filter(temperature < 3.35) %>%
  group_by(y_binary) %>% 
  count() %>% 
  ungroup() %>% 
  mutate(prop = round(n/sum(n), 2))

rpart.plot(classif_tree_gini, extra = 102, main = "Showing the classification rate at the node.")


# Using Entropy instead of Gini
classif_tree_entropy <- 
  rpart(y_binary ~., 
        data = df_train %>% select(-rented_bike_count, -date),
        method = "class",
        parms = list(split = "information"))
rpart.plot(classif_tree_entropy, main = "Classification tree build using entropy instead of gini to measure impurity index.")

# Similarly to what we did with regression trees, let us have a closer look at how the splitting rule is performed

variable_split <- "temperature"
threshold <- 10

# Root node as current node
current_node <- df_train



prop_current_node <- 
  current_node %>% 
  group_by(y_binary) %>% 
  count() %>% 
  ungroup() %>% 
  mutate(p_k = n / sum(n))
prop_current_node


# Gini impurity index in that current node
gini_node <- 
  prop_current_node %>% 
  mutate(g_k = p_k*(1-p_k)) %>% 
  magrittr::extract2("g_k") %>% 
  sum()
gini_node


#' @param data_node tibble/data.frame with the data of the node
#' @param target name of the target variable (discrete variable)
compute_gini <- function(data_node, target){
  
  prop_current_node <- 
    data_node %>% 
    group_by(!!sym(target)) %>% 
    count() %>% 
    ungroup() %>% 
    mutate(p_k = n / sum(n))
  
  prop_current_node %>% 
    mutate(g_k = p_k*(1-p_k)) %>% 
    magrittr::extract2("g_k") %>% 
    sum()
}


compute_gini(current_node, "y_binary")


left_node <- 
  current_node %>% 
  filter(temperature < threshold)

right_node <- 
  current_node %>% 
  filter(temperature >= threshold)


# proportions of observations in each node
p_L <- nrow(left_node) / nrow(current_node)
p_R <- nrow(right_node) / nrow(current_node)
cat(str_c("Prop left: ", p_L, "\n", "Prop right: ", p_R))


gini_left <- compute_gini(left_node, "y_binary")
gini_right <- compute_gini(right_node, "y_binary")

c(gini_left, gini_right)


# weighted average of the Gini impurity indices:
gini_l_r <- p_L*gini_left + p_R * gini_right
gini_l_r


# Delta
gini_node - gini_l_r


# pct Delta
(gini_l_r - gini_node)/gini_node


gini_split <- function(data_node, variable_to_predict, 
                       variable_split, threshold,
                       minbucket){
  
  if(is.numeric(data_node[[variable_split]])){
    left_node <- 
      data_node %>% 
      filter(!!sym(variable_split) < threshold)
    
    right_node <- 
      data_node %>% 
      filter(!!sym(variable_split) >= threshold)
  }else{
    left_node <- 
      data_node %>% 
      filter(!!sym(variable_split) %in% threshold)
    
    right_node <- 
      data_node %>% 
      filter( !( !!sym(variable_split) %in% threshold ))
  }
  
  # If there is less than a given number of obs in the leaves
  # warning
  warning_min_bucket <- 
    any( table(left_node[[variable_to_predict]]) < minbucket) |
    any( table(right_node[[variable_to_predict]]) < minbucket)
  
  # Proportions of obs in the resulting leaves
  p_L <- nrow(left_node) / nrow(current_node)
  p_R <- nrow(right_node) / nrow(current_node)
  
  # Gini in leaves
  gini_left <- compute_gini(left_node, "y_binary")
  gini_right <- compute_gini(right_node, "y_binary")
  
  gini_l_r <- p_L*gini_left + p_R * gini_right
  
  list(gini_l_r = gini_l_r, warning_min_bucket = warning_min_bucket)
}


gini_split(
  data_node = df_train, variable_to_predict = "y_binary",
  variable_split = "temperature", threshold = 10, minbucket = 7)


ggini_l_r_temperature <- 
  gini_leaves_df %>% 
  filter(warning_min_bucket == "Permitted") %>% 
  arrange(desc(delta))
gini_l_r_temperature

# Results provided by rpart:
nrow(current_node) * (gini_node-gini_l_r_temperature$g_leaves[1])


# Let us consider multiple cutoffs:

number_of_cuts <- 1000
thresholds <- 
  seq(min(df_train$temperature), max(df_train$temperature),
      length.out = number_of_cuts)


# loop over these values of cutoff.
gini_leaves <- 
  map(thresholds, ~gini_split(
    data_node = df_train, variable_to_predict = "y_binary",
    variable_split = "temperature", threshold = ., minbucket = 7))


# results stored in a table:
gini_leaves_df <- 
  tibble(
    g_leaves = map_dbl(gini_leaves, "gini_l_r"),
    warning_min_bucket = map_lgl(gini_leaves, "warning_min_bucket"),
    threshold = thresholds) %>% 
  mutate(delta = gini_node - g_leaves,
         warning_min_bucket = factor(
           warning_min_bucket, levels = c(TRUE, FALSE),
           labels = c("Not permitted", "Permitted")))
gini_leaves_df


ggplot(data = gini_leaves_df, 
       mapping = aes(x = threshold, y = delta)) +
  geom_point(mapping = aes(colour = warning_min_bucket)) +
  scale_colour_manual(
    "Threshold", 
    values = c("Permitted" = "blue", "Not permitted" = "red")) +
  labs(x = "Temperature threshold", title = "Decrease in the impurity design.") +
  theme(plot.title.position = "plot")

gini_leaves_df %>% 
  filter(warning_min_bucket == "Permitted") %>% 
  arrange(desc(delta))

# best threshold for the split made with temperature is the one, among the permitted solutions, 
# that leads to the highest impurity decrease:
gini_leaves_df %>% 
  filter(warning_min_bucket == "Permitted") %>% 
  arrange(desc(delta))


# Let us just give a few examples of other competitor variables
head(classif_tree_gini$splits, 5)

(gini_initial <- compute_gini(df_train, "y_binary"))

# increase in the goodness of fit using temperature as the splitting variable
gini_split_temperature <- 
  gini_split(
    data_node = df_train, variable_to_predict = "y_binary",
    variable_split = "temperature", threshold = 3.35, minbucket = 7)

nrow(df_train) * (gini_initial - gini_split_temperature$gini_l_r)


# Using seasons as the splitting variable
gini_split_seasons <- 
  gini_split(
    data_node = df_train, variable_to_predict = "y_binary",
    variable_split = "seasons", threshold = c("Winter"), minbucket = 7)

nrow(df_train) * (gini_initial - gini_split_seasons$gini_l_r)


# Using month:
gini_split_month <- 
  gini_split(
    data_node = df_train, variable_to_predict = "y_binary",
    variable_split = "month", threshold = c("Jan","Feb","Dec"), minbucket = 7)

nrow(df_train) * (gini_initial - gini_split_month$gini_l_r)


# 1.3.4 Variable Importance -----------------------------------------------

# Root node as current
current_node <- df_train

# Winter gives the best split

tmp <- 
  current_node %>% 
  mutate(
    # Primary split
    leave_primary = case_when(
      temperature < 3.35 ~ "left prim.",
      temperature >= 3.35 ~ "right prim."),
    # Surrogate split
    leave_surrogate = case_when(
      seasons == "Winter" ~ "left surr.",
      seasons != "Winter" ~ "right surr."
    ))


confusion <- table(tmp$leave_primary, tmp$leave_surrogate, exclude = NULL)
confusion



correct_direction <- (confusion["left prim.", "left surr."] + confusion["right prim.", "right surr."])
correct_direction

# Proportion
correct_direction/sum(confusion)


# majority rule
maj_rule_correct <- sum(confusion["right prim.", ])
maj_rule_correct


# adjusted agreement
adjust <- (correct_direction - maj_rule_correct) / (sum(confusion) - maj_rule_correct)
adjust


classif_tree_gini$variable.importance


# Relative importance as reported in the summary
round(100 * classif_tree_gini$variable.importance / 
        sum(classif_tree_gini$variable.importance))


splits_df <- 
  as_tibble(classif_tree_gini$splits, rownames = "variable")
splits_df


head(classif_tree_gini$frame)


# Position of the primary variables in `frame`
ind_primary <- which(classif_tree_gini$frame$var != "<leaf>")

type_var_df <- 
  classif_tree_gini$frame %>% 
  slice(ind_primary) %>% 
  select(ncompete, nsurrogate)

type_var <- NULL
for(i in 1:nrow(type_var_df)){
  type_var <- c(type_var, "primary", 
                rep("competitor", type_var_df[i, "ncompete"]),
                rep("surrogate", type_var_df[i, "nsurrogate"])
  ) 
}
type_var



splits_df <- splits_df %>% mutate(type_var = type_var)
splits_df


# focus once again on the variable temperature
variable_of_interest <- "temperature"


# Each primary split for this variable is credited with the improvement value:
var_imp <- 
  splits_df %>% 
  filter(variable == variable_of_interest, type_var == "primary") %>% 
  magrittr::extract2("improve") %>% 
  sum()
var_imp


# And each surrogate split where temperature is used gets split$adj times the primary split’s value. 
# Hence adding a column that gives the primary split’s value for each surrogate variable.
splits_df <- 
  splits_df %>% 
  mutate(primary_split_value = 
           ifelse(type_var == "primary", improve, NA)) %>% 
  fill(primary_split_value) %>% 
  mutate(primary_split_value = 
           ifelse(type_var != "surrogate", NA, primary_split_value))
splits_df


# contribution of the variable temperature each time it is a surrogate variable to the importance metric:
var_imp_surrogate <- 
  splits_df %>% 
  filter(variable == variable_of_interest,
         type_var == "surrogate") %>%
  mutate(imp_contrib = primary_split_value*adj) %>% 
  magrittr::extract2("imp_contrib") %>% 
  sum()
var_imp_surrogate


# add the two parts to obtain the importance measure of the temperature variable:
var_imp+var_imp_surrogate

# To compute the importance of all variables at once:
splits_df %>% 
  mutate(importance = ifelse(type_var == "primary", improve, 0),
         importance = ifelse(
           type_var == "surrogate", 
           importance+primary_split_value*adj, importance)) %>% 
  group_by(variable) %>% 
  summarise(importance = sum(importance)) %>% 
  arrange(desc(importance))


# Same values reported by rpart:
classif_tree_gini$variable.importance


as_tibble(classif_tree_gini$variable.importance,
          rownames = "variable") %>% 
  mutate(rel_importace = 100*value / sum(value)) %>% 
  ggplot(data = ., 
         mapping = aes(x = rel_importace, 
                       y = fct_reorder(variable, rel_importace))) +
  geom_bar(stat = "identity") +
  labs(x = "relative importance", y = NULL,
       title= "Variable importance for the Classification Tree.") +
  theme(plot.title.position = "plot")


# 1.4 Ensemble Methods ----------------------------------------------------



# 1.4.1 Bagging -----------------------------------------------------------


set.seed(123)
n <- 100
x <- runif(n=n, min=0, max=10)
eps <- rnorm(n, 0, 4)
f <- function(x) .1*x^3-1*x^2+x
y <- f(x)+eps
df_sim <- tibble(x=x, y=y)


# 100 bootstrap resamples of the data
n <- nrow(df_sim)
nb_bagg <- 100
baggin_models <- vector(mode = "list", length = nb_bagg)
for(i in 1:nb_bagg){
  
  ind_x <- sample(1:n, size = n, replace = TRUE)
  df_tmp <- df_sim %>% slice(ind_x)
  
  cart <- 
    rpart(y ~ x, 
          data = df_tmp,
          method = "anova",
          control = rpart.control(cp = 0)
    )
  
  baggin_models[[i]] <- cart
  
}

# The first weak learner
baggin_models[[1]]


# function that will make predictions on each weak learner on some data
# The predictions for a given point obtained with the weak learners will then be averaged
predict_bagged <- function(baggin_models, newdata){
  map(baggin_models, ~predict(., newdata = newdata)) %>% 
    bind_cols(
      .name_repair = ~ vctrs::vec_as_names(
        ..., repair = "unique", quiet = TRUE)) %>% 
    rowMeans()
}


new_data <- tibble(x = seq(0, 10, by = .1))


pred_bagging <- predict_bagged(baggin_models, new_data)
head(pred_bagging)


df_plot <- 
  new_data %>% 
  mutate(pred = pred_bagging,
         expected = f(x))

ggplot(mapping = aes(x = x, y = y)) +
  geom_point(data = df_sim) +
  geom_line(data = df_plot %>% pivot_longer(-x, values_to = "y"), 
            mapping = aes(colour = name),
            size = 1.1) +
  labs(title = "Bagging with 100 bootstrap training sets.") +
  theme(plot.title.position = "plot")


# 1.4.1.1 Number of Trees -------------------------------------------------

# function that samples from the data and then grow an unpruned tree
train_weak_learner <- function(data, formula, method, control){
  n <- nrow(data)
  ind_x <- sample(1:n, size = n, replace = TRUE)
  df_tmp <- data[ind_x,]
  cart <- rpart(formula = formula, 
                data = df_tmp,
                method = method,
                control = control)
  cart
}


library(foreach)
library(doSNOW)


ncl <- parallel::detectCores()-1
(cl <- makeCluster(ncl))
registerDoSNOW(cl)


baggin_models_2 <- foreach(
  i = 1:500,
  .packages=c("rpart")
  ) %dopar% {
  train_weak_learner(
    data = df_sim,
    formula = y~x, 
    method="anova", 
    control = rpart.control(cp = 0))
}

stopCluster(cl)


class(baggin_models_2)
length(baggin_models_2)

# predictions on our sequence of values ranging from 0 to 1
# using either :
# 
# - 10 trees
# - 75 trees
# - 250 trees
#  -500 trees

df_plot_2 <- 
  map(c(10, 75, 250, 500), ~add_column(
  new_data, 
  nb_bagg = as.character(.),
  pred = predict_bagged(baggin_models_2[1:.], new_data))) %>% 
  bind_rows()


ggplot() +
  geom_line(
    data = tibble(x=seq(0, 10, by = .1), y=f(x)),
    mapping = aes(x = x, y = y),
    colour = "black", size = 1.1, linetype = "dashed") +
  geom_point(data = df_sim, 
             mapping = aes(x = x, y = y)) +
  geom_line(data = df_plot_2, size = 1.1,
            mapping = aes(x = x, y = pred, colour = nb_bagg)) +
  scale_colour_discrete("Number of weak learners") +
  theme(plot.title.position = "plot") +
  labs(title = "Predictions made from the aggregations.")

# Unseen data
set.seed(234)
x_new <- runif(n=n, min=0, max=10)
eps <- rnorm(n, 0, 4)
y_new <- f(x)+eps
df_sim_new <- tibble(x=x_new, y=y_new)


# varying the number of weak learners
# before computing MSE
mse_test_bagging <- 
  map(1:500, ~add_column(
  df_sim_new, 
  nb_bagg = .,
  pred = predict_bagged(baggin_models_2[1:.], df_sim_new))) %>% 
  map(~mutate(., mse = (y-pred)^2)) %>% 
  map_df(~summarise(., nb_bagg=unique(nb_bagg), mse = mean(mse)))


ggplot(data = mse_test_bagging,
       mapping = aes(x = nb_bagg, y = mse)) +
  geom_line() +
  labs(x = "Number of trees (weak learners)", y = "MSE", title = "MSE vs number of trees.") +
  theme(plot.title.position = "plot")


# 1.4.1.2 Out-of-bag estimations ------------------------------------------


# Creating a weak learner
# and returning out-of-bag predictions

train_weak_learner <- function(data, formula, method, control){
  n <- nrow(data)
  ind_x <- sample(1:n, size = n, replace = TRUE)
  df_tmp <- data[ind_x,]
  cart <- rpart(formula = formula, 
                data = df_tmp,
                method = method,
                control = control)
  
  # oob predictions
  oob_pred <- predict(cart, newdata = data[-ind_x,])
  names(oob_pred) <- seq(1,n)[-ind_x]
  
  list(estim = cart, oob_pred = oob_pred)
}


# loop to build 500 trees
library(doSNOW)
ncl <- parallel::detectCores()-1
(cl <- makeCluster(ncl))
registerDoSNOW(cl)

baggin_models_3 <- foreach(
  i = 1:500,
  .packages=c("rpart")) %dopar% {
  
  train_weak_learner(
    data = df_sim,
    formula = y~x, 
    method="anova", 
    control = rpart.control(cp = 0))
}

stopCluster(cl)


# out-of-bag predictions
oob_pred <- 
  map(baggin_models_3, "oob_pred") %>% 
  map_df(~as_tibble(., rownames = "id")) %>% 
  bind_rows() %>% 
  mutate(id = as.numeric(id)) %>% 
  rename(pred = value) %>% 
  group_by(id) %>% 
  summarise(nb_seen = n(),
            pred = mean(pred))
oob_pred


# we can check how many times each observations has been used as an out-of-bag observation over the 500 trees:
summary(oob_pred$nb_seen)


# The out-of-bag MSE can then be computed:
df_sim %>% 
  mutate(id = row_number()) %>% 
  left_join(oob_pred) %>% 
  mutate(error_sq = (y-pred)^2) %>% 
  summarise(mse = mean(error_sq))


# Now, let us make the number of tree vary and look how it affects the out-of-bag MSE

#' Computes the OOB MSE using `nb_tree` trees
#' @param nb_tree number of trees used (weak learners)
compute_oob_mse_nb_trees <- function(nb_tree){
  oob_pred <- 
    map(baggin_models_3[1:nb_tree], "oob_pred") %>% 
    map_df(~as_tibble(., rownames = "id")) %>% 
    bind_rows() %>% 
    mutate(id = as.numeric(id)) %>% 
    rename(pred = value) %>% 
    group_by(id) %>% 
    summarise(nb_seen = n(),
              pred = mean(pred))
  
  df_sim %>% 
    mutate(id = row_number()) %>% 
    left_join(oob_pred, by = "id") %>% 
    mutate(error_sq = (y-pred)^2) %>% 
    summarise(mse = mean(error_sq, na.rm=TRUE)) %>% 
    mutate(nb_tree = nb_tree)
}


oob_mse_nb_trees <- map_df(1:300, compute_oob_mse_nb_trees)


ggplot(data = oob_mse_nb_trees,
       mapping = aes(x = nb_tree, y = mse)) +
  geom_line() +
  labs(x = "Number of trees grown (weak learners", y = "MSE",
       title = "Out-of-bag MSE - Bagging.") +
  theme(plot.title.position = "plot")



# 1.4.1.3 Variable Importance Measures ------------------------------------

# grow many unpruned trees.

library(doSNOW)
ncl <- parallel::detectCores()-1
(cl <- makeCluster(ncl))
registerDoSNOW(cl)

bagging_models_bike <- foreach(
  i = 1:200,
  .packages=c("rpart", "dplyr")
  ) %dopar% {
  
  train_weak_learner(
    data = df_train %>% select(-y_binary, -date),
    formula = rented_bike_count~., 
    method="anova",
    control = rpart.control(cp = 0))
}

stopCluster(cl)


# extract the variable importance measures from all these trees, and average the values over the trees:

variable_importance <- 
  map(bagging_models_bike, "estim") %>% 
  map("variable.importance") %>% 
  bind_rows() %>% 
  colMeans()


# relative measure
variable_importance_df <- 
  as_tibble(100*variable_importance/sum(variable_importance),
          rownames = "variable")
variable_importance_df


ggplot(data = variable_importance_df,
       mapping = aes(x = value, y = fct_reorder(variable, value))) +
  geom_bar(stat="identity")+
  labs(x = NULL, y = NULL, title = "Variable relative importance.") +
  theme(plot.title.position = "plot")



# 1.4.1.4 Pre-built Function ----------------------------------------------



library(ipred)

mod_bag <- bagging(
  formula = rented_bike_count ~ .,
  data = df_train %>% select(-y_binary, -date),
  nbagg = 100,  
  coob = TRUE,
  control = rpart.control(cp = 0)
)

# The trees are stored in the element mtrees
class(mod_bag$mtrees)
length(mod_bag$mtrees)

# the indices of the observations used to grow the tree (row number of the initial data)
names(mod_bag$mtrees[[1]])
# the estimated weak learner
head(mod_bag$mtrees[[1]]$bindx)

# out-of-bag error
mod_bag$err


# Variable importance
library(caret)
variable_importance <- varImp(mod_bag)
variable_importance


ggplot(data = as_tibble(variable_importance, rownames = "variable"),
       mapping = aes(x = Overall, y = fct_reorder(variable, Overall))) +
  geom_bar(stat="identity")+
  labs(x = NULL, y = NULL,
       title = "Variable relative importance.") +
  theme(plot.title.position = "plot")



# 1.4.2 Random Forests ----------------------------------------------------

# Trees may be correlated when use bagging.

tree_1_pruned <- prune(mod_bag$mtrees[[1]]$btree, cp=.1)
tree_2_pruned <- prune(mod_bag$mtrees[[2]]$btree, cp=.1)
tree_3_pruned <- prune(mod_bag$mtrees[[3]]$btree, cp=.1)
tree_4_pruned <- prune(mod_bag$mtrees[[4]]$btree, cp=.1)

op <- par()
par(mfrow = c(2,2))
rpart.plot(tree_1_pruned, main = "Tree 1")
rpart.plot(tree_2_pruned, main = "Tree 2")
rpart.plot(tree_3_pruned, main = "Tree 3")
rpart.plot(tree_4_pruned, main = "Tree 4")
par(op)


# 1.4.2.1 A First Example with {randomForest} -----------------------------

library(randomForest)


# Growing 200 trees
# 5 variables samples to perform each split

mod_rf <- 
  randomForest(
  formula = rented_bike_count ~ .,
  data = df_train %>% select(-y_binary, -date),
  ntree = 200,
  mtry = round(ncol(select(df_train, -y_binary, -date))/3),
  nodesize = 5,
  maxnodes = NULL
)


# out-of-bag MSE computed using 1 to 100 trees are stored in the mse elements of the model:
df_plot_mse_oob_rf <- 
  tibble(mse = mod_rf$mse) %>% 
  mutate(nb_tree = row_number())
df_plot_mse_oob_rf


# MSE obtained on a test sample
predictions_valid_tree_wise <- 
  predict(mod_rf, newdata = df_valid, predict.all = TRUE)


dim(predictions_valid_tree_wise$individual)

# All we need to do is to compute the MSE using up to nb_tree:

#' Computes the MSE using `nb_tree` trees
compute_mse_trees <- function(nb_tree){
  pred <- predictions_valid_tree_wise$individual[,1:nb_tree]
  if(nb_tree > 1) pred <- rowMeans(pred)
  tibble(mse = mean((df_valid$rented_bike_count - pred)^2), 
         nb_tree = nb_tree)
}


df_plot_mse_validation_rf <- map_df(1:200, compute_mse_trees)
df_plot_mse_validation_rf

# let us compare these MSE to that obtained with bagging.

mod_bag <- bagging(
  formula = rented_bike_count ~ .,
  data = df_train %>% select(-y_binary, -date),
  nbagg = 200,  
  coob = TRUE,
  control = rpart.control(cp = 0)
)


#' Computes out-ob-bag predictions for the ith tree
compute_oob_bagging <- function(i, model, data){
  n <- nrow(data)
  current_tree <- model$mtrees[[i]]$btree
  ind_x_current <- model$mtrees[[i]]$bindx
  oob_pred <- predict(current_tree, newdata = data[-ind_x_current,])
  as_tibble(oob_pred) %>% 
    mutate(ind_x = seq(1,n)[-ind_x_current], tree = i)
}

oob_bagging <- 
  map_df(1:length(mod_bag$mtrees), 
         ~compute_oob_bagging(.,mod_bag, df_train))
oob_bagging


#' Computes out-of-bag MSE using out-of-bag predictions
#' from `nb_tree` trees
compute_mse_oob_ntree <- function(nb_tree){
  
  oob_pred <- 
    oob_bagging %>% 
    filter(tree <= nb_tree) %>% 
    group_by(ind_x) %>% 
    summarise(oob_pred = mean(value))
  
  df_train %>% 
    mutate(ind_x = row_number()) %>% 
    left_join(
      oob_pred,
       by = "ind_x"
    ) %>% 
    mutate(error_sq = (rented_bike_count-oob_pred)^2) %>% 
    summarise(mse = mean(error_sq, na.rm=TRUE)) %>% 
    mutate(nb_tree = nb_tree)
  
}

df_plot_mse_oob_bagging <- 
  map_df(1:length(mod_bag$mtrees), compute_mse_oob_ntree)
df_plot_mse_oob_bagging


#' Computes predictions for the ith tree
#' (helper function)
compute_pred_tree <- function(i, model, newdata){
  current_tree <- model$mtrees[[i]]$btree
  pred <- predict(current_tree, newdata = newdata)
  as_tibble(pred, rownames = "ind_x") %>%
    mutate(tree = i)
}

pred_validation_bagging <-
  map_df(1:length(mod_bag$mtrees), 
         ~compute_pred_tree(., mod_bag, df_valid))


#' Compute the MSE on a test set (`df_validation`)
#' @param nb_tree number of trees used to compute the MSE
compute_mse_bagging_ntree <- function(nb_tree){

  pred <-
    pred_validation_bagging %>%
    mutate(ind_x = as.numeric(ind_x)) %>%
    filter(tree <= nb_tree) %>%
    group_by(ind_x) %>%
    summarise(oob_pred = mean(value))

  df_train %>%
    mutate(ind_x = row_number()) %>%
    left_join(
      pred,
      by = "ind_x"
    ) %>%
    mutate(error_sq = (rented_bike_count-oob_pred)^2) %>%
    summarise(mse = mean(error_sq, na.rm=TRUE)) %>%
    mutate(nb_tree = nb_tree)

}
df_plot_mse_validation_bagging <- 
  map_df(1:length(mod_bag$mtrees), compute_mse_bagging_ntree)
df_plot_mse_validation_bagging



df_plot <- 
  df_plot_mse_oob_bagging %>% 
  mutate(type = "OOB", model = "Bagging") %>% 
  bind_rows(
    df_plot_mse_validation_bagging %>% 
      mutate(type = "Test", model = "Bagging")
  ) %>% 
  bind_rows(
    df_plot_mse_oob_rf %>% 
      mutate(type = "OOB", model = "Random Forests")
  ) %>% 
  bind_rows(
    df_plot_mse_validation_rf %>% 
      mutate(type = "Test", model = "Random Forests")
  )

ggplot(data = df_plot,
       mapping = aes(x = nb_tree, y = mse)) +
  geom_line(mapping = aes(colour = model)) +
  labs(x = "Number of trees", y = "Error (MSE)",
       title = "Goodness of fit depending on the number of trees.") +
  theme(plot.title.position = "plot") +
  facet_wrap(~type, scales = "free_y", ncol = 2)



# 1.4.2.2 Varying the hyperparameters -------------------------------------



# 1.4.2.2.1 Varying mtry --------------------------------------------------


library(doSNOW)
ncl <- parallel::detectCores()-1
(cl <- makeCluster(ncl))
registerDoSNOW(cl)

# With a progress bar
pb <- txtProgressBar(min=1, max = 5, style = 3)
progress <- function(n) setTxtProgressBar(pb, n)
opts <- list(progress = progress)

nb_predictors <- ncol(df_train %>% select(-y_binary, -date))

mod_rf_mtry <- foreach(
  mtry=1:5,
  .packages=c("randomForest", "dplyr"),
  .options.snow = opts) %dopar% { # With probressbar
    randomForest(
      formula = rented_bike_count ~ .,
      data = df_train %>% select(-y_binary, -date),
      ntree = 200,
      mtry = mtry,
      nodesize = 5,
      maxnodes = NULL
    )
  }
close(pb)
stopCluster(cl)


# From the results, we can extract the out-of-bag MSE:
df_plot_mtry <- 
  map(mod_rf_mtry, "mse") %>% 
  map_df(~tibble(mse=.) %>% 
           mutate(nb_tree = row_number()), .id = "mtry")
df_plot_mtry


ggplot(data = df_plot_mtry,
       mapping = aes(x = nb_tree, y = mse)) +
  geom_line(mapping = aes(colour = mtry)) +
  labs(x = "Number of trees", y = "Error (MSE)",
       title = "Varying the number of variables samples to perform a split.")


# 1.4.2.2.2 Varying nodesize ----------------------------------------------

library(doSNOW)
ncl <- parallel::detectCores()-1
(cl <- makeCluster(ncl))
registerDoSNOW(cl)

pb <- txtProgressBar(min=1, max = 4, style = 3)
progress <- function(n) setTxtProgressBar(pb, n)
opts <- list(progress = progress)

nb_predictors <- ncol(df_train %>% select(-y_binary, -date))


mod_rf_nodesize <- foreach(
  nodesize = c(5, 10, 15, 20),
  .packages=c("randomForest", "dplyr"),
  .options.snow = opts # progressbar
  ) %dopar% {
    randomForest(
      formula = rented_bike_count ~ .,
      data = df_train %>% select(-y_binary, -date),
      ntree = 200,
      mtry = round(nb_predictors/3),
      nodesize = nodesize,
      maxnodes = NULL
    )
  }
# close(pb)
stopCluster(cl)


# out-of-bag errors can once
df_plot_nodesize <- 
  map(mod_rf_nodesize, "mse") %>% 
  map_df(~tibble(mse=.) %>% mutate(nb_tree = row_number()), .id = "ind_nodesize") %>% 
  mutate(ind_nodesize = as.numeric(ind_nodesize)) %>% 
  left_join(
    tibble(nodesize = c(5, 10, 15, 20)) %>% mutate(ind_nodesize = row_number()),
    by = "ind_nodesize"
  )


df_plot_nodesize

ggplot(data = df_plot_nodesize,
         mapping = aes(x = nb_tree, y = mse)) +
  geom_line(mapping = aes(colour = factor(nodesize))) +
  labs(x = "Number of trees", y = "Error (MSE)",
       title = "Varying the minimum size of terminal nodes.") +
  scale_colour_discrete("nodesize") +
  theme(plot.title.position = "plot")

# 1.4.2.2.3 Varying maxnodes ----------------------------------------------

library(doSNOW)
ncl <- parallel::detectCores()-1
(cl <- makeCluster(ncl))
registerDoSNOW(cl)

pb <- txtProgressBar(min=1, max = 6, style = 3)
progress <- function(n) setTxtProgressBar(pb, n)
opts <- list(progress = progress)

nb_predictors <- ncol(df_train %>% select(-y_binary, -date))

mod_rf_maxnodes <- foreach(
  max_nodes = c(5, 10, 20, 30, 40, 50),
  .packages=c("randomForest", "dplyr"),
  .options.snow = opts # with progressbar
  ) %dopar% {
    randomForest(
      formula = rented_bike_count ~ .,
      data = df_train %>% select(-y_binary, -date),
      ntree = 500,
      mtry = round(nb_predictors/3),
      nodesize = 5,
      maxnodes = max_nodes
    )
  }
close(pb)
stopCluster(cl)


# out-of-bag MSE
df_plot_maxnodes <- 
  map(mod_rf_maxnodes, "mse") %>% 
  map_df(~tibble(mse=.) %>%
           mutate(nb_tree = row_number()), .id = "ind_maxnodes") %>% 
  mutate(ind_maxnodes = as.numeric(ind_maxnodes)) %>% 
  left_join(
    tibble(maxnodes = c(5, 10, 20, 30, 40, 50)) %>% 
      mutate(ind_maxnodes = row_number()),
    by = "ind_maxnodes"
  )
df_plot_maxnodes


df_plot_maxnodes

ggplot(data = df_plot_maxnodes,
       mapping = aes(x = nb_tree, y = mse)) +
  geom_line(mapping = aes(colour = factor(maxnodes))) +
  labs(x = "Number of trees", y = "Error (MSE)",
       title = "Varying the maximum number of terminal nodes trees.") +
  scale_colour_discrete("maxnodes") +
  theme(plot.title.position = "plot")


# 1.4.2.3 Variable importance ---------------------------------------------


mod_rf <- 
  randomForest(
    formula = rented_bike_count ~ .,
    data = df_train %>% select(-y_binary, -date),
    ntree = 100,
    mtry = round(ncol(select(df_train, -y_binary, -date))/3),
    nodesize = 5,
    maxnodes = NULL,
    keep.inbag = TRUE,
    importance = TRUE
  )


compute_mse <- function(observed, predicted){
  mean((observed-predicted)^2)
}

# Let us focus on the first weak learner:
tree_number <- 1


# out-of-bag observations from the first tree
ind_x_oob <- 
  mod_rf$inbag[, tree_number][mod_rf$inbag[, tree_number] == 0] %>% 
  names() %>% 
  as.numeric()

df_oob <- df_train[ind_x_oob, ]


# first predictor:
predictor_names <- attr(mod_rf$terms, "term.labels")
j <- 1
predictor_names[j]


predictions_init <- 
  predict(mod_rf, newdata = df_oob, predict.all = TRUE)

# Then we can extract only the predicted values from the first tree, and compute the initial MSE, i.e., 
# the MSE computed without noising up the first predictor.

predictions_init_t <- predictions_init$individual[, tree_number]
mse_init_t <- compute_mse(df_oob$rented_bike_count, predictions_init_t)
mse_init_t


# noise-up the first predictor.
df_oob_noised_j <- df_oob %>% 
  mutate(!!sym(predictor_names[j]) := 
           sample(!!sym(predictor_names[j]), replace=FALSE))


cbind(init = df_oob[[predictor_names[j]]], 
      noised = df_oob_noised_j[[predictor_names[j]]]) %>% 
  head()


# out-of-bag MSE with this noised-up data can then be computed:
predictions_xj <-
  predict(mod_rf, newdata = df_oob_noised_j, predict.all = TRUE)
predictions_noised_j_t <- 
  predictions_xj$individual[, tree_number]
mse_noised_j_t <- 
  compute_mse(df_oob_noised_j$rented_bike_count, predictions_noised_j_t)
mse_noised_j_t


cbind(mse_before = mse_init_t, mse_after = mse_noised_j_t)


# function that recreate the out-of-bag data for a given tree
get_oob_data <- function(tree_number, model){
  # Indices of out-of-bag observations
  ind_x_oob <- 
    mod_rf$inbag[, tree_number][model$inbag[, tree_number] == 0] %>% 
    names() %>% 
    as.numeric()
  
  # oob data
  df_train[ind_x_oob, ]
}


# wrap-up the previous code to expand the computation to all trees and all variables


#' Computes tree-wise permutation importance
#' @param tree_number index of the t-th tree
#' @param model object of class randomForest
compute_pi_tree <- function(tree_number, model){
  
  # Out-of-bag data for tree `tree_number`
  df_oob <- get_oob_data(tree_number, model)
  
  # tree-wise oob MSE before noising-up predictors
  predictions_init <- predict(model, newdata = df_oob, predict.all = TRUE)
  predictions_init_t <- predictions_init$individual[, tree_number]
  mse_init_t <- compute_mse(df_oob$rented_bike_count, predictions_init_t)
  
  predictor_names <- attr(mod_rf$terms, "term.labels")
  
  # Looping over indices of predictors
  mse_noised_t <- rep(NA, length(predictor_names))
  for(j in seq_along(predictor_names)){
    # Noising-up data
    df_oob_noised_j <- df_oob %>% 
      mutate(!!sym(predictor_names[j]) := sample(
        !!sym(predictor_names[j]),
        replace=FALSE))
    
    # tree-wise oob MSE after noising-up jth predictor
    predictions_xj <- 
      predict(model, newdata = df_oob_noised_j, predict.all = TRUE)
    predictions_noised_j_t <- 
      predictions_xj$individual[, tree_number]
    mse_noised_j_t <- 
      compute_mse(df_oob_noised_j$rented_bike_count, predictions_noised_j_t)
    mse_noised_t[[j]] <- mse_noised_j_t
  }
  
  # tree-wise permuation importance
  tree_wise_pi <- mse_noised_t-mse_init_t
  names(tree_wise_pi) <- predictor_names
  tree_wise_pi
}


# tree-wise permuation importance of the variables computed over the first tree
permutation_importance_t <- compute_pi_tree(tree_number = 1, model = mod_rf)
permutation_importance_t


# tree-wise permuation importance values for all trees
tree_wise_pi <- 
  map(1:mod_rf$ntree, ~compute_pi_tree(., model = mod_rf))


# results in a table
tree_wise_pi_df <- 
  tree_wise_pi %>% 
  map_df(~as_tibble(., rownames = "variable"), .id="tree")
tree_wise_pi_df


# forest-wise permutation importance of each variable
permutation_importance <- 
  tree_wise_pi_df %>% 
  group_by(variable) %>% 
  summarise(permutation_importance = mean(value)) %>% 
  arrange(desc(permutation_importance))


ggplot(data = permutation_importance,
       mapping = aes(x = permutation_importance,
                     y = fct_reorder(variable, permutation_importance))) +
  geom_bar(stat = "identity") +
  labs(x = "Permutation importance", y = NULL)


# Values reported by {randomForest}
ggplot(data = as_tibble(mod_rf$importance, rownames = "variable"),
       mapping = aes(x = IncNodePurity,
                     y = fct_reorder(variable, IncNodePurity))) +
  geom_bar(stat = "identity") +
  labs(x = "Variable importance (Increase in node Purity)", y = NULL) +


# Permutation importance computed using importance()
imp <- importance(mod_rf, type=1, scale = F)

ggplot(data = as_tibble(imp, rownames = "variable"),
       mapping = aes(x = `%IncMSE`,
                     y = fct_reorder(variable, `%IncMSE`))) +
  geom_bar(stat = "identity") +
  labs(x = "Permutation importance (% Increase in MSE)", y = NULL)


# 1.4.2.4 Fine Tuning -----------------------------------------------------



number_predictors <- 
  df_train %>% select(-y_binary, -date) %>% ncol()
number_predictors


hyperparameters_grid <- 
  expand_grid(
    ntree = c(number_predictors*10, number_predictors*20),
    mtry = c(3,5),
    nodesize = c(3, 5, 10),
    maxnodes = c(-1, 10, 50, 100, 200)
  )
hyperparameters_grid



library(doSNOW)
ncl <- parallel::detectCores()-1
(cl <- makeCluster(ncl))
registerDoSNOW(cl)

# Progress bar
pb <- txtProgressBar(min=1, max = nrow(hyperparameters_grid), style = 3)
progress <- function(n) setTxtProgressBar(pb, n)
opts <- list(progress = progress)

mse_search <- foreach(
  i=1:nrow(
    hyperparameters_grid),
  .packages=c("randomForest", "dplyr"),
  .options.snow = opts # with progressbar
  ) %dopar% {
    
    current_hyper <- hyperparameters_grid[i,]
    
    mod_rf_tmp <- 
      randomForest(
        formula = rented_bike_count ~ .,
        data = df_train %>% select(-y_binary, -date),
        ntree = current_hyper$ntree,
        mtry = current_hyper$mtry,
        nodesize = current_hyper$nodesize,
        maxnodes = current_hyper$maxnodes
      )
    
    current_hyper %>% mutate(mse = mean(mod_rf_tmp$mse))
  }
close(pb)
stopCluster(cl)



mse_search <- 
  bind_rows(mse_search) %>% 
  arrange(mse)
mse_search


best_hyper_param <- mse_search %>% slice(1)

# "best" random forest can then be grown again:

mod_rf_best <- 
  randomForest(
    formula = rented_bike_count ~ .,
    data = df_train %>% select(-y_binary, -date),
    ntree = best_hyper_param$ntree,
    mtry = best_hyper_param$mtry,
    nodesize = best_hyper_param$nodesize,
    maxnodes = best_hyper_param$maxnodes
  )



pred_rf_best <- predict(mod_rf_best, newdata = df_valid)
mse_rf_best <- 
  compute_mse(observed = df_valid$rented_bike_count,
              predicted = pred_rf_best)

pred_bagging <- predict(mod_bag, newdata = df_valid)
mse_bagging <- 
  compute_mse(observed = df_valid$rented_bike_count,
              predicted = pred_bagging)

mse_bagging
mse_rf_best