# Machine learning for econometricians # Course for the ECB # by Ewen Gallic # December 2021 # Hands-on Session 1 # Overfitting # 1 First Example: Default of Credit Card --------------------------------- library(tidyverse) library(readxl) url <- "https://egallic.fr/Enseignement/ML/ECB/data/default_credit_card_clients.xls" download.file(url = url, destfile = "default_credit_card_clients.xls") df <- read_excel("default_credit_card_clients.xls", skip = 1) # Rename the target variable: df <- df %>% rename(y = `default payment next month`) df # 1.1 Somme Summary Statistics on the Whole Dataset ----------------------- # Dimensions dim(df) # Unbalanced dataset df %>% group_by(y) %>% count() %>% arrange(desc(n)) # Barplot df %>% mutate(y=factor(y, levels = c(1,0), labels = c("Yes", "No"))) %>% ggplot(data = .) + geom_bar(stat = "count", aes(y = y)) + labs(x = "Frequency", y = NULL, title = "Default payment the next month") + theme(plot.title.position = "plot") # Table with variables, description and type c( "y", "Default payment the next month", "qualitative", "ID", "ID number", "numeric", "LIMIT_BAL", "Amount of the given credit", "numeric", "SEX", "Gender", "qualitative", "EDUCATION", "Education", "ordinal", "MARRIAGE", "Marital status", "qualitative", "AGE", "Age", "numeric", "PAY_0", "Repayment status September", "qualitative", "PAY_2", "Repayment status August", "qualitative", "PAY_3", "Repayment status July", "qualitative", "PAY_4", "Repayment status June", "qualitative", "PAY_5", "Repayment status May", "qualitative", "PAY_6", "Repayment status April", "qualitative", "BILL_AMT1", "Bill amount September", "numeric", "BILL_AMT2", "Bill amount August", "numeric", "BILL_AMT3", "Bill amount July", "numeric", "BILL_AMT4", "Bill amount June", "numeric", "BILL_AMT5", "Bill amount May", "numeric", "BILL_AMT6", "Bill amount April", "numeric", "PAY_AMT1", "Amount previous payment September", "numeric", "PAY_AMT2", "Amount previous payment August", "numeric", "PAY_AMT3", "Amount previous payment July", "numeric", "PAY_AMT4", "Amount previous payment June", "numeric", "PAY_AMT5", "Amount previous payment May", "numeric", "PAY_AMT6", "Amount previous payment April", "numeric") %>% matrix(ncol = 3, byrow = TRUE) %>% as_tibble() %>% magrittr::set_colnames(c("variable", "label", "type")) -> variable_names_df #' format_p_value #' Format the p-value so that it displays "<10^{i}", if i <= 3 #' rounds the value with 2 digits otherwise #' @param x p-value to format #' x <- 0.000133613 ; x <- 0.0011 format_p_value <- function(x){ if(x < 10^-3){ resul <- str_c("**< 10^{-3}**") }else if(x > 10^-3){ if(x < 0.05){ resul <- str_c("**", round(x, 3), "**") }else{ resul <- round(x, 3) } }else{ resul <- "**0.001**" } resul } #' Compute descriptive statistics on `df` grouped by `grouping_var` #' depending on the type of the variable of interest (qualitative, ordinal, numerical) #' Then performs test of equality of mean (or proportion) on the sub-samples. #' #' @param variable (string) name of the variable of interest #' @param grouping_var (string) name of the qualitative variable that defines groups #' @param variable_names table with three columns: 1- variable name, 2- label, 3-type #' variable <- "LIMIT_BAL" sample_diff <- function(variable, df, grouping_var, variable_names = variable_names_tlm){ label_variable <- variable_names %>% filter(variable == !!variable) %>% magrittr::extract2("label") type <- variable_names %>% filter(variable == !!variable) %>% magrittr::extract2("type") if(type %in% c("qualitative", "ordinal")){ whole_sample <- df %>% filter(!is.na(!!!syms(variable)), !is.na(!!!syms(grouping_var))) %>% mutate(label = label_variable) %>% select(label, !!!variable) %>% group_by(label, !!!syms(variable)) %>% summarise(n = n(), .groups = "drop") nb_tot <- sum(whole_sample$n) whole_sample <- whole_sample %>% mutate(pct = n / sum(n), pct = (round(100*pct, 2)) %>% str_c("(", ., "%)")) %>% unite(`Whole Sample`, n, pct, sep = " ") subsamples <- df %>% filter(!is.na(!!!syms(variable)), !is.na(!!!syms(grouping_var))) %>% mutate(label = label_variable) %>% select(label, !!!variable, !!!grouping_var) %>% group_by(label, !!!syms(variable), !!!(syms(grouping_var))) %>% summarise(n = n(), .groups = "drop") %>% group_by(!!!syms(grouping_var)) %>% mutate(pct = n / sum(n), pct = (round(100*pct, 2)) %>% str_c("(", ., "%)")) %>% unite(value, n, pct, sep = " ") res <- whole_sample %>% left_join(subsamples) %>% ungroup() %>% spread_(key_col = grouping_var, value_col = "value", fill = "0 (0.00%)") %>% mutate(label = ifelse(row_number() == 1, yes = str_c(label, " (n = ", format(nb_tot, big.mark = ","), ")"), no = label)) %>% unite(label, label, !!!variable, sep = " ") }else{ whole_sample <- df %>% filter(!is.na(!!!syms(variable)), !is.na(!!!syms(grouping_var))) %>% mutate(label = label_variable) %>% select(label, !!!variable) nb_tot <- nrow(whole_sample) whole_sample <- whole_sample %>% group_by(label) %>% summarise(mean = mean(!!!syms(variable), na.rm = TRUE) %>% formatC(format = "e", digits = 2), sd = sd(!!!syms(variable), na.rm = TRUE) %>% formatC(format = "e", digits = 2), .groups = "drop") %>% mutate(sd = str_c("(+/-", sd, ")")) %>% unite(`Whole Sample`, mean, sd, sep = " ") subsamples <- df %>% filter(!is.na(!!!syms(variable)), !is.na(!!!syms(grouping_var))) %>% mutate(label = label_variable) %>% select(label, !!!variable, !!!grouping_var) %>% group_by(label, !!!syms(grouping_var)) %>% summarise(mean = mean(!!!syms(variable)) %>% formatC(format = "e", digits = 2), sd = sd(!!!syms(variable)) %>% formatC(format = "e", digits = 2), .groups = "drop") %>% mutate(sd = str_c("(+/-", sd, ")")) %>% unite(value, mean, sd, sep = " ") res <- whole_sample %>% left_join(subsamples) %>% ungroup() %>% spread_(key_col = grouping_var, value_col = "value", fill = "0 (+/-0.00)") %>% mutate(label = ifelse(row_number() == 1, yes = str_c(label, " (n = ", format(nb_tot, big.mark = ","), ")"), no = label)) } if(type == "qualitative"){ chisq_test <- chisq.test(df %>% magrittr::extract2(grouping_var), df %>% magrittr::extract2(variable)) p_value <- chisq_test$p.value }else if(type == "ordinal"){ kruskal_test <- kruskal.test(str_c(grouping_var, "~", variable) %>% as.formula(), data = df) p_value <- kruskal_test$p.value }else{ # Numerical form <- str_c(variable, " ~ ", grouping_var) %>% as.formula() anov <- aov(form, data = df) anov_summary <- anov %>% summary() %>% .[[1]] p_value <- anov_summary$`Pr(>F)`[1] } # Adding stars for p-values p_value <- format_p_value(p_value) resul <- res %>% mutate(p_value = c(rep("", nrow(res)-1), p_value)) resul }# End of sample_diff() # Prepare table df_summary_stats <- variable_names_df$variable[!variable_names_df$variable %in% c("y", "ID")] %>% map_df(sample_diff, df = df, grouping_var = "y", variable_names = variable_names_df) # Print the summary table: library(kableExtra) df_summary_stats %>% kableExtra::kable(caption = "Driving habits before and after the claim") %>% kableExtra::kable_classic(full_width = F, html_font = "Cambria") %>% kableExtra::kable_styling( bootstrap_options = c("striped", "hover", "condensed", "responsive")) # 1.2 Fitting the Model --------------------------------------------------- # Create a dataset with only 50 observations here. set.seed(134) # Number of observations n <- 50 df_first_example <- sample_n(df %>% filter(y == 0), size = n/2, replace= FALSE) %>% bind_rows( sample_n(df %>% filter(y == 1), size = n/2, replace= FALSE) ) %>% select(-ID) # Fitting a logistic regression using all 23 predictors: reg <- glm(y ~ ., data = df_first_example, family = "binomial") #' Get the confusion table #' @param observed vector of observed values #' @param predicted vector of predicted values confusion_table <- function(observed, predicted){ confusion_matrix <- table(observed, predicted) confusion_matrix_prop <- prop.table(confusion_matrix) confusion_matrix_print <- confusion_matrix for(i in 1:ncol(confusion_matrix_print)){ for(j in 1:nrow(confusion_matrix_print)){ confusion_matrix_print[i,j] <- str_c(confusion_matrix_print[i,j], " (", confusion_matrix_prop[i,j], "%)") } } metrics <- compute_metric_confusion(confusion_matrix) list(confusion_matrix = confusion_matrix, confusion_matrix_print = confusion_matrix_print, metrics = metrics) } #' Computes errors metrics from the confusion matrix #' @param confusion_matrix Confusion matrix compute_metric_confusion <- function(confusion_matrix){ true_negative <- confusion_matrix["0", "0"] if(is.null(true_negative)) true_negative <- 0 false_positive <- confusion_matrix["0", "1"] if(is.null(false_positive)) false_positive <- 0 false_negative <- confusion_matrix["1", "0"] if(is.null(false_negative)) false_negative <- 0 true_positive <- confusion_matrix["1", "1"] if(is.null(true_positive)) true_positive <- 0 error_rate <- (false_positive+false_negative)/sum(confusion_matrix) # True positive rate TPR <- true_positive / (true_positive + false_negative) # False positive rate FPR <- false_positive / (false_positive + true_negative) # True negative rate TNR <- true_negative / (true_negative + false_positive) # False negative rate FNR <- false_negative / (false_negative + true_positive) # Overall error rate OER <- (false_negative + false_positive) / (false_negative + false_positive + true_positive + true_negative) list(`Erorr rate` = error_rate, `True positive rate` = TPR, `False positive rate` = FPR, `True negative rate` = TNR, `False negative rate` = FNR, `Overall error rate` = OER ) } # Predictive capabilities of our model: # Predictions: threshold <- .5 pred_prob <- predict(reg) pred <- ifelse(pred_prob>.5, 1, 0) confusion_matrix <- confusion_table(df_first_example$y, pred) confusion_matrix$confusion_matrix_print confusion_matrix$metrics # 1.3 Randomly Assigning the Classes -------------------------------------- set.seed(123) df_first_example <- df_first_example %>% mutate(y_new = sample(c(0, 1), size = n(), replace = TRUE)) df_first_example$y df_first_example$y_new # Half of the "observed" values have been modified here sum(df_first_example$y == df_first_example$y_new) # Fitting the model: reg_2 <- glm(y_new ~ ., data = df_first_example %>% select(-y), family = "binomial") # Predictive abilities: pred_prob_2 <- predict(reg_2) pred_2 <- ifelse(pred_prob_2>.5, 1, 0) conf_matrix_2 <- confusion_table(df_first_example$y_new, pred_2) conf_matrix_2$confusion_matrix_print conf_matrix_2$metrics # Repeat this random switching of the “true” class. # At each iteration, we will store the results in a list. random_fit <- function(){ df_first_example <- df_first_example %>% mutate(y_new = sample(c(0, 1), size = n(), replace = TRUE)) # Let us train the model again, this time on this randomly drawn response reg_2 <- glm(y_new ~ ., data = df_first_example %>% select(-y), family = "binomial") pred_prob_2 <- predict(reg_2) pred_2 <- ifelse(pred_prob_2>.5, 1, 0) conf_matrix_2 <- confusion_table(df_first_example$y_new, pred_2) conf_matrix_2 } results_random <- vector(mode = "list", length = 100) pb <- txtProgressBar(min = 1, max = 100, style = 3) for(i in 1:100){ results_random[[i]] <- random_fit() setTxtProgressBar(pb, i) } # overall error rates of each of the 100 repetitions: error_rates <- map(results_random, "metrics") %>% map_dbl("Overall error rate") summary(error_rates) ggplot(data = tibble(error_rate = error_rates)) + geom_boxplot(aes(x =error_rate)) + labs(x = "Overall error rate", title = "Overall error rate over 100 repetitions") # 1.4 Performances on Unseen Data ----------------------------------------- # Some unseen data: df_test <- sample_n(df %>% filter(y == 0), size = n/2, replace= FALSE) %>% bind_rows( sample_n(df %>% filter(y == 1), size = n/2, replace= FALSE) ) %>% select(-ID) # Predictions on these data, using our model: pred_prob_test <- predict(reg, newdata=df_test) pred_test <- ifelse(pred_prob_test>.5, 1, 0) # Confusion matrix on these unseen data: confusion_matrix_test <- confusion_table(df_test$y, pred_test) confusion_matrix_test$confusion_matrix_print confusion_matrix_test$metrics # The predictive ability of the model worsened. # Could it be bad luck? set.seed(123) # Sample of size 200 n <- 200 df_first_example_2 <- sample_n(df %>% filter(y == 0), size = n/2, replace= FALSE) %>% bind_rows( sample_n(df %>% filter(y == 1), size = n/2, replace= FALSE) ) %>% select(-ID) # Splitting into train/test to evaluate the capacity of the model # to generalise its results on unseen data: n_train <- round(.8*nrow(df_first_example_2)) ind_train <- sample(1:nrow(df_first_example_2), size = n_train, replace = FALSE) df_train <- df_first_example_2 %>% slice(ind_train) df_test <- df_first_example_2 %>% slice(-ind_train) dim(df_train) dim(df_test) # Fitting the model reg_train <- glm(y ~ ., data = df_train, family = "binomial") # Predictions on seen data: pred_train_prob <- predict(reg_train) pred_train <- ifelse(pred_train_prob>.5, 1, 0) # Confusion matrix: confusion_matrix_train <- confusion_table(df_train$y, pred_train) confusion_matrix_train$confusion_matrix_print confusion_matrix_train$metrics # Prediction on unseen data: pred_test_prob <- predict(reg_train, newdata = df_test) pred_test <- ifelse(pred_test_prob>.5, 1, 0) confusion_matrix_test <- confusion_table(df_test$y, pred_test) confusion_matrix_test$confusion_matrix_print confusion_matrix_test$metrics # We will repeat this procedure 100 times. # To make things easier, we can define a function that # -randomly constitutes a training and a testing datasets from df_first_example_2 # - fits the model # - then returns the confusion matrices as well as the goodness of fit metrics. train_test_estim <- function(){ n_train <- round(.8*nrow(df_first_example_2)) ind_train <- sample(1:nrow(df_first_example_2), size = n_train, replace = FALSE) df_train <- df_first_example_2 %>% slice(ind_train) df_test <- df_first_example_2 %>% slice(-ind_train) reg_train <- glm(y ~ ., data = df_train, family = "binomial") # Predictions on the train set pred_train_prob <- predict(reg_train) pred_train <- ifelse(pred_train_prob>.5, 1, 0) # On the test set pred_test_prob <- predict(reg_train, newdata = df_test) pred_test <- ifelse(pred_test_prob>.5, 1, 0) # Confusion matrices confusion_matrix_train <- confusion_table(df_train$y, pred_train) confusion_matrix_test <- confusion_table(df_test$y, pred_test) list(confusion_matrix_train = confusion_matrix_train, confusion_matrix_test = confusion_matrix_test) } results_random_2 <- vector(mode = "list", length = 100) pb <- txtProgressBar(min = 1, max = 100, style = 3) for(i in 1:100){ results_random_2[[i]] <- train_test_estim() setTxtProgressBar(pb, i) } # The overall error rates obtained on both sets over the 100 iterations: error_rates_train_2 <- map(results_random_2, "confusion_matrix_train") %>% map("metrics") %>% map_dbl("Overall error rate") error_rates_test_2 <- map(results_random_2, "confusion_matrix_test") %>% map("metrics") %>% map_dbl("Overall error rate") summary(error_rates_train_2) summary(error_rates_test_2) # Let us have a look at these overall error rates on boxplots: ggplot(data = tibble(error_rate = c(error_rates_train_2, error_rates_test_2), sample = rep(c("Train", "Test"), each = 100))) + geom_boxplot(aes(x =error_rate, y = sample)) + labs(x = "Overall Error rate", y = "Sample", title = "Model performance on 100 random draws of the data") # 1.5 K-fold Cross Validation --------------------------------------------- k <- 5 #' Splits a vector into folds #' @param x vector of observations to be splitted #' @param k number of desired folds split_into_folds <- function(x,k) split(x, cut(seq_along(x), k, labels = FALSE)) # Example of use: split_into_folds(seq(1, 20), k=3) # stratified sampling to make sure we end with approximatively the same distribution # of the target variable as that of the dataset. (classes <- unique(df_first_example_2$y)) folds <- vector(mode = "list", length = length(classes)) for(i in 1:length(classes)){ class <- classes[i] # row number of obs with target equal to !!class # shuffle the obs ind_class <- sample(which(df_first_example_2$y == class), replace=FALSE) # split into k folds folds[[i]] <- split_into_folds(ind_class, k=k) } folds <- map(1:k, ~map(folds, .)) %>% map(unlist) %>% map(sample) folds # We will consider each fold iteratively and make the estimation of the modle on the remaing folds. # At each iteration, we will store the overall error rates. error_rates <- NULL for(current_fold in 1:k){ # All the obs but those of the kth fold ind_obs_current <- unlist(folds[-current_fold]) df_train_current <- df_first_example_2 %>% slice(ind_obs_current) df_test_current <- df_first_example_2 %>% slice(-ind_obs_current) # Fitting the model reg_train <- glm(y ~ ., data = df_train_current, family = "binomial") # Predictions on the train set pred_train_prob <- predict(reg_train) pred_train <- ifelse(pred_train_prob>.5, 1, 0) # On the test set pred_test_prob <- predict(reg_train, newdata = df_test_current) pred_test <- ifelse(pred_test_prob>.5, 1, 0) # Confusion matrices confusion_matrix_train <- confusion_table(df_train_current$y, pred_train) confusion_matrix_test <- confusion_table(df_test_current$y, pred_test) error_rate_train <- confusion_matrix_train$metrics$`Overall error rate` error_rate_test <- confusion_matrix_test$metrics$`Overall error rate` error_rates[[current_fold]] <- list(error_rate_train = error_rate_train, error_rate_test = error_rate_test) } # Errors on Train set (seen data) map_dbl(error_rates, "error_rate_train") map_dbl(error_rates, "error_rate_train") %>% mean() # Errors on Test set (unseen data) map_dbl(error_rates, "error_rate_test") map_dbl(error_rates, "error_rate_test") %>% mean() # 1.6 Leave-one-out Cross Validation -------------------------------------- # Leave one out CV error_rates_loocv <- NULL for(i in 1:nrow(df_first_example_2)){ df_train_current <- df_first_example_2 %>% slice(-i) # Single obs: df_test_current <- df_first_example_2 %>% slice(i) # Fitting the model reg_train <- glm(y ~ ., data = df_train_current, family = "binomial") # Predictions on the train set pred_train_prob <- predict(reg_train) pred_train <- ifelse(pred_train_prob>.5, 1, 0) # On the test set pred_test_prob <- predict(reg_train, newdata = df_test_current) pred_test <- ifelse(pred_test_prob>.5, 1, 0) avg_error_train <- mean(!pred_train == df_train_current$y) avg_error_test <- mean(!pred_test == df_test_current$y) error_rates_loocv[[i]] <- list(avg_error_train = avg_error_train, avg_error_test = avg_error_test) } # Errors on train test map_dbl(error_rates_loocv, "avg_error_train") map_dbl(error_rates_loocv, "avg_error_train") %>% mean() # Errors on test test map_dbl(error_rates_loocv, "avg_error_test") map_dbl(error_rates_loocv, "avg_error_test") %>% mean() # 2 Second Example: Selling Price of Cars --------------------------------- cars_df <- read_csv("https://egallic.fr/Enseignement/ML/ECB/data/cars.csv") cars_df summary(cars_df$selling_price) sd(cars_df$selling_price) # Correlations cars_df %>% select(where(is_double)) %>% cor(use = "complete.obs") # Proportions of each category for discrete variables cars_df %>% select(selling_price, where(~is.character(.x) | is.factor(.x)), -name) %>% pivot_longer(cols = -selling_price) %>% group_by(name, value) %>% count() %>% group_by(name) %>% mutate(prop = round(100*n / sum(n), 2)) %>% arrange(name, desc(prop)) # Recoding some values cars_df <- cars_df %>% mutate(fuel = forcats::fct_lump_min(fuel, min = 100, other_level = "Other"), owner = forcats::fct_recode(owner, "First Owner" = "Test Drive Car")) # There are some missing values: cars_df %>% summarise(across(everything(), ~sum(is.na(.)))) %>% data.frame() %>% t() # quick and dirty method: replace NA by sample average cars_df <- cars_df %>% # replace NA values with average in data # (this method is really dirty and should be discussed) mutate( engine = replace_na(engine, mean(engine, na.rm=TRUE)), max_power = replace_na(max_power, mean(max_power, na.rm=TRUE)), seats = replace_na(seats, mean(seats, na.rm=TRUE)), mileage = ifelse(mileage < 0, yes = NA, no = mileage), mileage = replace_na(mileage, mean(mileage, na.rm=TRUE)) ) %>% # Transformation on data (will not be used after in this notebook, though) mutate(selling_price_log = log(selling_price), km_driven_log = log(km_driven), mileage_log = log(mileage+1), engine_log = log(engine), max_power_log = log(max_power+1)) # Let us only keep 500 observations to fasten the computations set.seed(456) cars_df <- cars_df %>% sample_n(1000) library(randomForest) # Defining training and validation samples set.seed(456) n_train <- round(.8*nrow(cars_df)) ind_train <- sample(1:nrow(cars_df), size = n_train, replace = FALSE) df_train <- cars_df %>% slice(ind_train) df_valid <- cars_df %>% slice(-ind_train) # Sanity check: dim(df_train) dim(df_valid) # 2.1 Predicting the Price with a Random Forest --------------------------- set.seed(456) mtry <- 4 model_1 <- randomForest::randomForest( selling_price ~ year + fuel + km_driven + seller_type + transmission + owner + mileage + engine + max_power + seats, data=df_train, importance=TRUE, ntree=50, mtry=mtry, nodesize = 10) predicted_train_1 <- predict(model_1, newdata = df_train) # MSE on train sample resid_train_1 <- df_train$selling_price - predicted_train_1 mse_train_1 <- mean(resid_train_1^2) mse_train_1 # With anotehr value for mtry # > a higher value can make the model more flexible and could thus reduce the MSE. set.seed(456) mtry <- 7 model_2 <- randomForest::randomForest( selling_price ~ year + fuel + km_driven + seller_type + transmission + owner + mileage + engine + max_power + seats, data=df_train, importance=TRUE, ntree=50, mtry=mtry, nodesize = 10) # MSE on train test predicted_train_2 <- predict(model_2, newdata = df_train) resid_train_2 <- df_train$selling_price - predicted_train_2 (mse_train_2 <- mean(resid_train_2^2)) # Comparing MSE: predicted_test_1 <- predict(model_1, newdata = df_valid) predicted_test_2 <- predict(model_2, newdata = df_valid) resid_test_1 <- df_valid$selling_price - predicted_test_1 resid_test_2 <- df_valid$selling_price - predicted_test_2 mse_test_1 <- mean(resid_test_1^2) mse_test_2 <- mean(resid_test_2^2) # Reminder on the values obtained on the train set: mse_train_1 mse_train_2 # With the test set: mse_test_1 # less flexible model mse_test_2 # more flexible model # Let us have a broader look at how the MSE varies as long as we make the model more or less flexible. # To do so, we will simply use a loop. # At each iteration, we will train the random forest on the same data, but we will make the number of # variables randomly sampled as candidates at each split vary. set.seed(123) mtry <- c(2,3,4,5,6,7,8) pb <- txtProgressBar(min = 1, max = length(mtry), style = 3) mse <- NULL for(j in 1:length(mtry)){ model <- randomForest::randomForest( selling_price ~ year + fuel + km_driven + seller_type + transmission + owner + mileage + engine + max_power + seats, data=df_train, importance=TRUE, ntree=50, mtry=mtry[j], nodesize = 10) predicted_train <- predict(model, newdata = df_train) predicted_test <- predict(model, newdata = df_valid) resid_train <- df_train$selling_price - predicted_train resid_test <- df_valid$selling_price - predicted_test mse_train <- mean(resid_train^2) mse_test <- mean(resid_test^2) mse <- mse %>% bind_rows(tibble(mse_train = mse_train, mse_test = mse_test, mtry = mtry[j])) setTxtProgressBar(pb, j) } mse mse %>% pivot_longer(cols = -mtry, names_to = "sample", values_to = "MSE") %>% mutate(sample = factor(sample, levels = c("mse_train", "mse_test"), labels = c("Train", "Test"))) %>% ggplot(data = ., mapping = aes(x = mtry, y = MSE)) + geom_line(mapping = aes(colour = sample)) + labs(x = "Number of variables randomly sampled as candidates at each split", y = "MSE", title = "Mean Squared Error on Train/Test data.") # 2.2 Cross-validation to Select the Hyperparameters ---------------------- nb_folds <- 5 #' Splits a vector into folds #' @param x vector of observations to be splitted #' @param k number of desired folds split_into_folds <- function(x,k) split(x, cut(seq_along(x), k, labels = FALSE)) set.seed(123) fold_ind <- split_into_folds(sample(1:nrow(df_train)), nb_folds) mse <- NULL mtry <- c(2,3,4,5,6,7,8) for(j in 1:length(mtry)){ for(k in 1:nb_folds){ df_train_current <- df_train %>% slice(unlist(fold_ind[-k])) df_test_current <- df_train %>% slice(fold_ind[[k]]) model <- randomForest::randomForest( selling_price ~ year + fuel + km_driven + seller_type + transmission + owner + mileage + engine + max_power + seats, data=df_train, importance=TRUE, ntree=50, mtry=mtry[j], nodesize = 10) predicted_train <- predict(model, newdata = df_train_current) predicted_test <- predict(model, newdata = df_test_current) resid_train <- df_train_current$selling_price - predicted_train resid_test <- df_test_current$selling_price - predicted_test mse_train <- mean(resid_train^2) mse_test <- mean(resid_test^2) mse <- mse %>% bind_rows(tibble(mse_train = mse_train, mse_test = mse_test, k = k, mtry = mtry[j])) } setTxtProgressBar(pb, j) } # The average MSE for each value of mtry: mse_summary <- mse %>% group_by(mtry) %>% summarise( mse_train = mean(mse_train), mse_test = mean(mse_test) ) mse_summary mse_summary %>% arrange(mse_train) mse_summary %>% arrange(mse_test) # 3 Third Example: Choice of Lambda in Lasso Regression ------------------- library(tidyverse) library(ISLR) # Removing NA values Hitters <- na.omit(Hitters) x <- model.matrix(Salary ~., Hitters)[,-1] head(x) y <- Hitters$Salary head(y) # Training / validation sets n_train <- round(.8*nrow(x)) ind_train <- sample(1:nrow(x), size = n_train, replace=FALSE) x_train <- x[ind_train,] x_valid <- x[-ind_train,] y_train <- y[ind_train] y_valid <- y[-ind_train] # Grid search: # Different values of lambda will be tested grid <- 10^seq(2, -2, length = 100) grid #' Cpmputes the Mean Squared Error #' @param observed vector of observed data #' @param predicted vector of predicted values, same length as `observed` compute_mse <- function(observed, predicted){ mean((observed - predicted)^2) } # 3.1 First Method: Without Cross-Validation ------------------------------ # Preparing list which will contain the results mse_lasso <- vector(mode = "list", length = length(grid)) library(glmnet) for(i in 1:length(grid)){ lambda <- grid[i] lasso_m_current <- glmnet(x = x_train, y = y_train, alpha=1, lambda = lambda, standardize = TRUE) predictions_current <- predict(lasso_m_current, newx = x_train) mse_lasso[[i]] <- tibble(lambda = lambda, mse = compute_mse(y_train, predictions_current)) } mse_lasso <- bind_rows(mse_lasso) %>% arrange(mse) mse_lasso best_lambda_1 <- mse_lasso$lambda[1] ggplot(data = mse_lasso, mapping = aes(x = lambda, y = mse)) + geom_line() + labs(title = "Mean Squared Error depending on the value of lambda") # 3.2 Second Method: With Cross-Validation -------------------------------- nb_folds <- 10 n <- nrow(x_train) folds <- sample(rep(1:nb_folds, length=n)) # List that will contain the results mse <- vector(mode = "list", length = nb_folds) for(k in 1:nb_folds){ ind_current <- folds != k # Train set x_train_current <- x_train[ind_current,] y_train_current <- y_train[ind_current] # Test set x_test_current <- x_train[-ind_current,] y_test_current <- y_train[-ind_current] # Looping over the values of lambda mse_fold <- vector(mode = "list", length = length(grid)) for(i in 1:length(grid)){ lambda <- grid[i] # Estimating on the training set lasso_m_current <- glmnet(x = x_train_current, y = y_train_current, alpha=1, lambda = lambda, standardize = TRUE) # Predictions on the test set pred_test <- predict(lasso_m_current, x_test_current) # MSE mse_test_current <- compute_mse(y_test_current, pred_test) # Storing the results mse_fold[[i]] <- tibble(nb_folds = k, lambda = lambda, mse_test = mse_test_current) } mse[[k]] <- bind_rows(mse_fold) } mse <- bind_rows(mse) mse # average and standard deviation of MSE on the k folds: mse_summary <- mse %>% group_by(lambda) %>% summarise( mse_sd_test = sd(mse_test), mse_test = mean(mse_test) ) best_lambda_cv <- mse_summary %>% arrange(mse_test) %>% slice(1) # computing how much is 1-SE away best_lambda_cv <- best_lambda_cv %>% mutate(mse_test_one_sd = mse_test+mse_sd_test) # Finding the candidates: # values of lambda for which the MSE is 1-SE away best_lambda_candidates <- mse_summary %>% filter(mse_test <= best_lambda_cv$mse_test_one_sd) %>% arrange(desc(lambda)) best_lambda_candidates # Highest lambda among the candidates: best_lambda_2 <- best_lambda_candidates$lambda[1] best_lambda_2 mse_summary %>% ggplot(data = mse_summary, mapping = aes(x = lambda, y = mse_test)) + geom_line() + labs(x = "Lambda", y = "MSE", title = "Average MSE computed on the left-aside fold") # 3.3 Comparing the Capacities of the Models on Unseen Data --------------- # Model where the value of lambda was selected based on # the MSE computed using all the observations at hand lasso_model_1 <- glmnet(x = x_train, y = y_train, alpha=1, lambda = best_lambda_1, standardize = TRUE) # Model where the value of lambda was selected using cross-validation lasso_model_2 <- glmnet(x = x_train, y = y_train, alpha=1, lambda = best_lambda_2, standardize = TRUE) # Predictions on unseen data for both models: pred_valid_1 <- predict(lasso_model_1, x_valid) mse_valid_1 <- compute_mse(y_valid, pred_valid_1) pred_valid_2 <- predict(lasso_model_2, x_valid) mse_valid_2 <- compute_mse(y_valid, pred_valid_2) mse_valid_1 mse_valid_2 # 3.4 Repeating the Comparison 100 times ---------------------------------- ---- eval=FALSE---------------------------------------------------------------------------- nb_repeat <- 100 mse_simulation <- vector(mode = "list", length = nb_repeat) # pb <- txtProgressBar(min = 1, max = nb_repeat, style = 3) for(i_repeat in 1:nb_repeat){ # Training and validation sets n_train <- round(.8*nrow(x)) ind_train <- sample(1:nrow(x), size = n_train, replace=FALSE) x_train <- x[ind_train,] x_valid <- x[-ind_train,] y_train <- y[ind_train] y_valid <- y[-ind_train] # 1. Lambda estimated without cross-validation mse_lasso <- vector(mode = "list", length = length(grid)) for(i in 1:length(grid)){ lambda <- grid[i] lasso_m_current <- glmnet(x = x_train, y = y_train, alpha=1, lambda = lambda, standardize = TRUE) predictions_current <- predict(lasso_m_current, newx = x_train) mse_lasso[[i]] <- tibble(lambda = lambda, mse = compute_mse(y_train, predictions_current)) } mse_lasso <- bind_rows(mse_lasso) %>% arrange(mse) best_lambda_1 <- mse_lasso$lambda[1] #2. Lambda estimated with cross-validation folds <- sample(rep(1:nb_folds, length=n)) mse <- vector(mode = "list", length = nb_folds) for(k in 1:nb_folds){ ind_current <- folds != k # Train set x_train_current <- x_train[ind_current,] y_train_current <- y_train[ind_current] # Test set x_test_current <- x_train[-ind_current,] y_test_current <- y_train[-ind_current] # Looping over the values of lambda mse_fold <- vector(mode = "list", length = length(grid)) for(i in 1:length(grid)){ lambda <- grid[i] # Estimating on the training set lasso_m_current <- glmnet(x = x_train_current, y = y_train_current, alpha=1, lambda = lambda, standardize = TRUE) # Predictions on the test set pred_test <- predict(lasso_m_current, x_test_current) # MSE mse_test_current <- compute_mse(y_test_current, pred_test) # Storing the results mse_fold[[i]] <- tibble(nb_folds = k, lambda = lambda, mse_test = mse_test_current) } mse[[k]] <- bind_rows(mse_fold) } mse <- bind_rows(mse) # Compute the average and standard deviation of MSE on the k folds mse_summary <- mse %>% group_by(lambda) %>% summarise( mse_sd_test = sd(mse_test), mse_test = mean(mse_test) ) # What is the value of lambda that gives the lowest MSE? best_lambda_cv <- mse_summary %>% arrange(mse_test) %>% slice(1) %>% mutate(mse_test_one_sd = mse_test+mse_sd_test) # Lambdas within 1 standard error best_lambda_candidates <- mse_summary %>% filter(mse_test <= best_lambda_cv$mse_test_one_sd) %>% arrange(desc(lambda)) # Among the candidates, let us pick the highest best_lambda_2 <- best_lambda_candidates$lambda[1] # 3. Performance on the validation set # Model where the value of lambda was selected based on # the MSE computed using all the observations at hand lasso_model_1 <- glmnet(x = x_train, y = y_train, alpha=1, lambda = best_lambda_1, standardize = TRUE) # Model where the value of lambda was selected using cross-validation lasso_model_2 <- glmnet(x = x_train, y = y_train, alpha=1, lambda = best_lambda_2, standardize = TRUE) # Predictions: pred_valid_1 <- predict(lasso_model_1, x_valid) mse_valid_1 <- compute_mse(y_valid, pred_valid_1) pred_valid_2 <- predict(lasso_model_2, x_valid) mse_valid_2 <- compute_mse(y_valid, pred_valid_2) mse_simulation[[i_repeat]] <- tibble(i = i_repeat, mse_valid_1 = mse_valid_1, mse_valid_2) # setTxtProgressBar(pb, i_repeat) } mse_simulation <- bind_rows(mse_simulation) mse_simulation # how many replications yielded lower MSE on the validation set when lambda # was selected by cross-validation: mse_simulation %>% mutate(lower_with_cv = mse_valid_2 < mse_valid_1) %>% group_by(lower_with_cv) %>% count() # descriptive statistics: mse_simulation %>% pivot_longer(cols = -i) %>% group_by(name) %>% summarise(mean = mean(value), sd = sd(value), q_1 = quantile(value, probs = .25), med = quantile(value, probs = .5), q_3 = quantile(value, probs = .75) ) # boxplots mse_simulation %>% pivot_longer(cols = -i, names_to = "sample", values_to = "MSE") %>% mutate(sample = factor(sample, levels = c("mse_valid_1", "mse_valid_2"), labels = c("Without CV", "With CV"))) %>% ggplot(data = ., mapping = aes(y = MSE)) + geom_boxplot(mapping = aes(fill = sample)) + labs(title = "MSE on the replications depending on how lambda was selected") + scale_fill_discrete("Choice of lambda")