Chapter 2: Modeling Process

Note: Some results may differ from the hard copy book due to the changing of sampling procedures introduced in R 3.6.0. See http://bit.ly/35D1SW7 for more details. Access and run the source code for this notebook here.

Hidden chapter requirements used in the book to set the plotting theme and load packages used in hidden code chunks.

knitr::opts_chunk$set(
  echo = TRUE,
  fig.align = "center",
  message = FALSE,
  warning = FALSE,
  collapse = TRUE,
  cache = FALSE
)

# underlying code dependencies
library(tidyverse)

# Set the graphical theme
theme_set(theme_light())

Figure 2.1:

knitr::include_graphics("images/modeling_process.png")

Prerequisites

This chapter leverages the following packages.

Data used:

# Ames housing data
ames <- AmesHousing::make_ames()
ames.h2o <- as.h2o(ames)

# Job attrition data
churn <- rsample::attrition %>% 
  mutate_if(is.ordered, .funs = factor, ordered = FALSE)
churn.h2o <- as.h2o(churn)

Data splitting

Figure 2.2:

knitr::include_graphics("images/data_split.png")

Simple random sampling

# Using base R
set.seed(123)  # for reproducibility
index_1 <- sample(1:nrow(ames), round(nrow(ames) * 0.7))
train_1 <- ames[index_1, ]
test_1  <- ames[-index_1, ]

# Using caret package
set.seed(123)  # for reproducibility
index_2 <- createDataPartition(ames$Sale_Price, p = 0.7, 
                               list = FALSE)
train_2 <- ames[index_2, ]
test_2  <- ames[-index_2, ]

# Using rsample package
set.seed(123)  # for reproducibility
split_1  <- initial_split(ames, prop = 0.7)
train_3  <- training(split_1)
test_3   <- testing(split_1)

# Using h2o package
split_2 <- h2o.splitFrame(ames.h2o, ratios = 0.7, 
                          seed = 123)
train_4 <- split_2[[1]]
test_4  <- split_2[[2]]

Figure 2.3:

p1 <- ggplot(train_1, aes(x = Sale_Price)) + 
    geom_density(trim = TRUE) + 
    geom_density(data = test_1, trim = TRUE, col = "red") +
  ggtitle("Base R")

p2 <- ggplot(train_2, aes(x = Sale_Price)) + 
    geom_density(trim = TRUE) + 
    geom_density(data = test_2, trim = TRUE, col = "red") +
    theme(axis.title.y = element_blank(),
          axis.ticks.y = element_blank(),
          axis.text.y = element_blank()) +
    ggtitle("caret") 

p3 <- ggplot(train_3, aes(x = Sale_Price)) + 
    geom_density(trim = TRUE) + 
    geom_density(data = test_3, trim = TRUE, col = "red") +
    theme(axis.title.y = element_blank(),
          axis.ticks.y = element_blank(),
          axis.text.y = element_blank()) +
    ggtitle("rsample")

p4 <- ggplot(as.data.frame(train_4), aes(x = Sale_Price)) + 
    geom_density(trim = TRUE) + 
    geom_density(data = as.data.frame(test_4), trim = TRUE, col = "red") +
    theme(axis.title.y = element_blank(),
          axis.ticks.y = element_blank(),
          axis.text.y = element_blank()) +
    ggtitle("h2o")

# Side-by-side plots
gridExtra::grid.arrange(p1, p2, p3, p4, nrow = 1)

# clean up
rm(p1, p2, p3, p4)

Stratified sampling

# orginal response distribution
table(churn$Attrition) %>% prop.table()

       No       Yes 
0.8387755 0.1612245 

# stratified sampling with the rsample package
set.seed(123)
split_strat  <- initial_split(churn, prop = 0.7, 
                              strata = "Attrition")
train_strat  <- training(split_strat)
test_strat   <- testing(split_strat)

# consistent response ratio between train & test
table(train_strat$Attrition) %>% prop.table()

      No      Yes 
0.838835 0.161165 

table(test_strat$Attrition) %>% prop.table()

       No       Yes 
0.8386364 0.1613636 

Many formula interfaces

This code chunk is for illustrative purposes only; it is not designed to execute.

# Sale price as function of neighborhood and year sold
model_fn(Sale_Price ~ Neighborhood + Year_Sold, 
         data = ames)

# Variables + interactions
model_fn(Sale_Price ~ Neighborhood + Year_Sold + 
           Neighborhood:Year_Sold, data = ames)

# Shorthand for all predictors
model_fn(Sale_Price ~ ., data = ames)

# Inline functions / transformations
model_fn(log10(Sale_Price) ~ ns(Longitude, df = 3) + 
           ns(Latitude, df = 3), data = ames)

This code chunk is for illustrative purposes only; it is not designed to execute.

# Use separate inputs for X and Y
features <- c("Year_Sold", "Longitude", "Latitude")
model_fn(x = ames[, features], y = ames$Sale_Price)

This code chunk is for illustrative purposes only; it is not designed to execute.

model_fn(
  x = c("Year_Sold", "Longitude", "Latitude"),
  y = "Sale_Price",
  data = ames.h2o
)

Many engines

lm_lm <- lm(Sale_Price ~ ., data = ames)
lm_glm <- glm(Sale_Price ~ ., data = ames, family = gaussian)
lm_caret <- train(Sale_Price ~ ., data = ames, method = "lm")

Table 1:

Table 1: Syntax for computing predicted class probabilities with direct engines.

Algorithm	Package	Code
Linear discriminant analysis	MASS	predict(obj)
Generalized linear model	stats	predict(obj, type = "response")
Mixture discriminant analysis	mda	predict(obj, type = "posterior")
Decision tree	rpart	predict(obj, type = "prob")
Random Forest	ranger	predict(obj)$predictions
Gradient boosting machine	gbm	predict(obj, type = "response", n.trees)

Resampling methods

k-fold cross validation

Figure 2.4:

knitr::include_graphics("images/cv.png")

Figure 2.5:

cv <- vfold_cv(mtcars, 10)
cv_plot <- cv$splits %>%
  purrr::map2_dfr(seq_along(cv$splits), ~ mtcars %>% mutate(
    Resample = paste0("Fold_", stringr::str_pad(.y, 2, pad = 0)),
    ID = row_number(),
    Data = ifelse(ID %in% .x$in_id, "Training", "Validation"))
    ) %>%
  ggplot(aes(Resample, ID, fill = Data)) +
  geom_tile() +
  scale_fill_manual(values = c("#f2f2f2", "#AAAAAA")) +
  scale_y_reverse("Observation ID", breaks = 1:nrow(mtcars), expand = c(0, 0)) +
  scale_x_discrete(NULL, expand = c(0, 0)) +
  theme_classic() +
  theme(legend.title=element_blank())

cv_plot

This code chunk is for illustrative purposes only; it is not designed to execute.

# Example using h2o
h2o.cv <- h2o.glm(
  x = x, 
  y = y, 
  training_frame = ames.h2o,
  nfolds = 10  # perform 10-fold CV
)

Bootstrapping

Figure 2.6:

knitr::include_graphics("images/bootstrap-scheme.png")

Figure 2.7:

boots <- rsample::bootstraps(mtcars, 10)
boots_plot <- boots$splits %>%
  purrr::map2_dfr(seq_along(boots$splits), ~ mtcars %>% 
             mutate(
               Resample = paste0("Bootstrap_", stringr::str_pad(.y, 2, pad = 0)),
               ID = row_number()
             ) %>%
             group_by(ID) %>%
             mutate(Replicates = factor(sum(ID == .x$in_id)))) %>%
  ggplot(aes(Resample, ID, fill = Replicates)) +
  geom_tile() +
  scale_fill_manual(values = c("#FFFFFF", "#F5F5F5", "#C8C8C8", "#A0A0A0", "#707070", "#505050", "#000000")) +
  scale_y_reverse("Observation ID", breaks = 1:nrow(mtcars), expand = c(0, 0)) +
  scale_x_discrete(NULL, expand = c(0, 0)) +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
  ggtitle("Bootstrap sampling") 

cv_plot <- cv_plot + 
  ggtitle("10-fold cross validation") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

cowplot::plot_grid(boots_plot, cv_plot, align = "h", nrow = 1)

# clean up
rm(boots, boots_plot, cv_plot)

Bias variance trade-off

Bias

Figure 2.8:

# Simulate some nonlinear monotonic data
set.seed(123)  # for reproducibility
x <- seq(from = 0, to = 2 * pi, length = 500)
y <- sin(x) + rnorm(length(x), sd = 0.3)
df <- data.frame(x, y) %>%
  filter(x < 4.5)

# Single model fit
bias_model <- lm(y ~ I(x^3), data = df)
df$predictions <- predict(bias_model, df)
p1 <- ggplot(df, aes(x, y)) +
  geom_point(alpha = .3) +
  geom_line(aes(x, predictions), size = 1.5, color = "dodgerblue") +
  scale_y_continuous("Response", limits = c(-1.75, 1.75), expand = c(0, 0)) +
  scale_x_continuous(limits = c(0, 4.5), expand = c(0, 0)) +
  ggtitle("Single biased model fit")

# Bootstrapped model fit
bootstrap_n <- 25
bootstrap_results <- NULL
for(i in seq_len(bootstrap_n)) {
  set.seed(i)  # for reproducibility
  index <- sample(seq_len(nrow(df)), nrow(df), replace = TRUE)
  df_sim <- df[index, ]
  fit <- lm(y ~ I(x^3), data = df_sim)
  df_sim$predictions <- predict(fit, df_sim)
  df_sim$model <- paste0("model", i)
  df_sim$ob <- index
  bootstrap_results <- rbind(bootstrap_results, df_sim)
}

p2 <- ggplot(bootstrap_results, aes(x, predictions, color = model)) +
  geom_line(show.legend = FALSE, size = .5) +
  scale_y_continuous("Response", limits = c(-1.75, 1.75), expand = c(0, 0)) +
  scale_x_continuous(limits = c(0, 4.5), expand = c(0, 0)) +
  ggtitle("25 biased models fit to bootstrap samples")

gridExtra::grid.arrange(p1, p2, nrow = 1)

Variance

Figure 2.9:

# Simulate some nonlinear monotonic data
set.seed(123)  # for reproducibility
x <- seq(from = 0, to = 2 * pi, length = 500)
y <- sin(x) + rnorm(length(x), sd = 0.3)
df <- data.frame(x, y) %>%
  filter(x < 4.5)

# Single model fit
variance_model <- knnreg(y ~ x, k = 3, data = df)
df$predictions <- predict(variance_model, df)
p1 <- ggplot(df, aes(x, y)) +
  geom_point(alpha = .3) +
  geom_line(aes(x, predictions), size = 1.5, color = "dodgerblue") +
  scale_y_continuous("Response", limits = c(-1.75, 1.75), expand = c(0, 0)) +
  scale_x_continuous(limits = c(0, 4.5), expand = c(0, 0)) +
  ggtitle("Single high variance model fit")

# Bootstrapped model fit
bootstrap_n <- 25
bootstrap_results <- NULL
for(i in seq_len(bootstrap_n)) {
  set.seed(i)  # for reproducibility
  index <- sample(seq_len(nrow(df)), nrow(df), replace = TRUE)
  df_sim <- df[index, ]
  fit <- knnreg(y ~ x, k = 3, data = df_sim)
  df_sim$predictions <- predict(fit, df_sim)
  df_sim$model <- paste0("model", i)
  df_sim$ob <- index
  bootstrap_results <- rbind(bootstrap_results, df_sim)
}

p2 <- ggplot(bootstrap_results, aes(x, predictions, color = model)) +
  geom_line(show.legend = FALSE) +
  scale_y_continuous("Response", limits = c(-1.75, 1.75), expand = c(0, 0)) +
  scale_x_continuous(limits = c(0, 4.5), expand = c(0, 0)) +
  ggtitle("25 high variance models fit to bootstrap samples")

gridExtra::grid.arrange(p1, p2, nrow = 1)

Hyperparameter tuning

Figure 2.10:

k_results <- NULL
k <- c(2, 5, 10, 20, 50, 150)

# Fit many different models
for(i in seq_along(k)) {
  df_sim <- df
  fit <- knnreg(y ~ x, k = k[i], data = df_sim)
  df_sim$predictions <- predict(fit, df_sim)
  df_sim$model <- paste0("k = ", stringr::str_pad(k[i], 3, pad = " "))
  k_results <- rbind(k_results, df_sim)
}

ggplot() +
  geom_point(data = df, aes(x, y), alpha = .3) +
  geom_line(data = k_results, aes(x, predictions), color = "dodgerblue", size = 1.5) +
  scale_y_continuous("Response", limits = c(-1.75, 1.75), expand = c(0, 0)) +
  scale_x_continuous(limits = c(0, 4.5), expand = c(0, 0)) +
  facet_wrap(~ model)

Figure 2.11:

cv <- trainControl(method = "repeatedcv", number = 10, repeats = 10, returnResamp = "all")
hyper_grid <- expand.grid(k = seq(2, 150, by = 2))
knn_fit <- train(x ~ y, data = df, method = "knn", trControl = cv, tuneGrid = hyper_grid)

ggplot() +
  geom_line(data = knn_fit$results, aes(k, RMSE)) +
  geom_point(data = knn_fit$results, aes(k, RMSE)) +
  geom_point(data = filter(knn_fit$results, k == as.numeric(knn_fit$bestTune)),
             aes(k, RMSE),
             shape = 21,
             fill = "yellow",
             color = "black",
             stroke = 1,
             size = 2) +
  scale_y_continuous("Error (RMSE)")

# clean up
rm(cv, hyper_grid, knn_fit)

Classification models

Figure 2.12

knitr::include_graphics("images/confusion-matrix.png")

Figure 2.13:

knitr::include_graphics("images/confusion-matrix2.png")

Figure 2.14:

library(plotROC)

# Generate data
set.seed(123)
response <- rbinom(200, size = 1, prob = .5)

set.seed(123)
curve1   <- rnorm(200, mean = response, sd = .40)

set.seed(123)
curve2   <- rnorm(200, mean = response, sd = .75)

set.seed(123)
curve3   <- rnorm(200, mean = response, sd = 2.0)

df <- tibble(response, curve1, curve2, curve3)

ggplot(df) + 
  geom_roc(aes(d = response, m = curve1), n.cuts = 0, size = .5, color = "#1E56F9") + 
  geom_roc(aes(d = response, m = curve2), n.cuts = 0, size = .5, color = "#7194F9") + 
  geom_roc(aes(d = response, m = curve3), n.cuts = 0, size = .5, color = "#B6C7F9") +
  geom_abline(lty = 'dashed') +
  annotate("text", x = .48, y = .46, label = c("No better than guessing"), 
           vjust = 1, angle = 34) +
  annotate("text", x = .3, y = .6, label = c("Ok"), 
           vjust = 1, angle = 33, color = "#B6C7F9") +
  annotate("text", x = .20, y = .75, label = c("Better"), 
           vjust = 1, angle = 33, color = "#7194F9") +
  annotate("text", x = .10, y = .96, label = c("Best"), 
           vjust = 1, angle = 33, color = "#1E56F9") +
  xlab("False positive rate") +
  ylab("True positive rate")

Putting the processes together

# Stratified sampling with the rsample package
set.seed(123)
split <- initial_split(ames, prop = 0.7, strata = "Sale_Price")
ames_train  <- training(split)
ames_test   <- testing(split)

This grid search takes approximately 3.5 minutes:

# Specify resampling strategy
cv <- trainControl(
  method = "repeatedcv", 
  number = 10, 
  repeats = 5
)

# Create grid of hyperparameter values
hyper_grid <- expand.grid(k = seq(2, 25, by = 1))

# Tune a knn model using grid search
knn_fit <- train(
  Sale_Price ~ ., 
  data = ames_train, 
  method = "knn", 
  trControl = cv, 
  tuneGrid = hyper_grid,
  metric = "RMSE"
)

Figure 2.15:

# Print and plot the CV results
knn_fit

k-Nearest Neighbors 

2053 samples
  80 predictor

No pre-processing
Resampling: Cross-Validated (10 fold, repeated 5 times) 
Summary of sample sizes: 1848, 1848, 1848, 1847, 1849, 1847, ... 
Resampling results across tuning parameters:

  k   RMSE      Rsquared   MAE     
   2  47844.53  0.6538046  31002.72
   3  45875.79  0.6769848  29784.69
   4  44529.50  0.6949240  28992.48
   5  43944.65  0.7026947  28738.66
   6  43645.76  0.7079683  28553.50
   7  43439.07  0.7129916  28617.80
   8  43658.35  0.7123254  28769.16
   9  43799.74  0.7128924  28905.50
  10  44058.76  0.7108900  29061.68
  11  44304.91  0.7091949  29197.78
  12  44565.82  0.7073437  29320.81
  13  44798.10  0.7056491  29475.33
  14  44966.27  0.7051474  29561.70
  15  45188.86  0.7036000  29731.56
  16  45376.09  0.7027152  29860.67
  17  45557.94  0.7016254  29974.44
  18  45666.30  0.7021351  30018.59
  19  45836.33  0.7013026  30105.50
  20  46044.44  0.6997198  30235.80
  21  46242.59  0.6983978  30367.95
  22  46441.87  0.6969620  30481.48
  23  46651.66  0.6953968  30611.48
  24  46788.22  0.6948738  30681.97
  25  46980.13  0.6928159  30777.25

RMSE was used to select the optimal model using the smallest value.
The final value used for the model was k = 7.

ggplot(knn_fit)

Shutdown H2O and clean up

h2o.shutdown(prompt = FALSE)

[1] TRUE

rm(list = ls())