Chapter 24: Model Building

Overview

This is a reproduction of the workflow for model building, which comes from Chapter 24 of R for Data Science. The purpose was to learn the concepts. I stick pretty close to the book, but you may find me taking a few creative liberties.

The diamond analysis walks through the process of building models with patterns that are hidden or confounded by other more prevalent patterns. For the diamond prices, carat weight confounds many of the other trends, and these trends are pulled out through the model building process.

The flight analysis walks through the process of building a model to predict flight counts by day over a time series. Two approaches are used. The first develops the model using day-of-the-week and seasonality trends to fit a model. The second (alternative) approach uses day-of-the-week and splines for seasonality. We use MASS::rlm() for the final model because it is much more resistant to outliers than lm().

---

Prerequisites

library(modelr)
options(na.action = na.warn)
library(tidyverse) # ggplot2, purrr, dplyr, tidyr, readr, tibble
library(lubridate)
library(nycflights13)

---

Diamond Analysis

Answer question: Why are low quality diamonds more expensive?

ggplot(diamonds, aes(cut, price)) + geom_boxplot()

ggplot(diamonds, aes(color, price)) + geom_boxplot()

ggplot(diamonds, aes(clarity, price)) + geom_boxplot()

Price and Carat

Explore price and relationship to carat. Use hex scatter plot to show concentration.

ggplot(diamonds, aes(carat, price)) + 
  geom_hex(bins = 50)

Need to transform data to get linear trend:

A log transformation will create a linear pattern that can be easily modelled. Hadley uses log base 2 and filters to remove carat size > 2.5. Log base 2 is needed for graph to turn out when performing the re-transform. Filtering out diamonds above size 2.5 is needed because the predictions become way off beyond 2.5 carats.

# remove very large diamonds (outliers), log transform price and carat
diamonds2 <- diamonds %>%
    filter(carat <= 2.5) %>%
    mutate(lprice = log2(price), 
           lcarat = log2(carat))
# Visualize linear pattern of log transform
ggplot(diamonds2, aes(lcarat, lprice)) +
    geom_hex(bins = 50)

Fit the model and visualize. Note that to visualize the data set, we re-transform to get price after the predictions to log2(price) are made.

# Model
mod_diamond <- lm(lprice ~ lcarat, data = diamonds2)
# Create data_grid() along range of carat, log transform and add_predicitons()
grid <- diamonds2 %>%
    data_grid(carat = seq_range(carat, 20)) %>%
    mutate(lcarat = log2(carat)) %>%
    add_predictions(mod_diamond, "lprice") %>%
    mutate(price = 2 ^ lprice)
# Visualize re-transformed prediction
ggplot(diamonds2, aes(carat, price)) + 
    geom_hex(bins = 50) + 
    geom_line(data = grid, color = "red", size = 1)

Check residuals to make sure linear pattern is removed that was due to carat size.

# add_residuals()
diamonds2 <- diamonds2 %>%
    add_residuals(mod_diamond, 'lresid')
# Visulize residuals
ggplot(diamonds2, aes(lcarat, lresid)) + 
    geom_hex(bins = 50)

Plot the residuals as a function of cut, color, and clarity. Plotting against residuals removes the effect of carat size, which dominated cut, color, and clarity effects.

We now see the effect we expected:

ggplot(diamonds2, aes(cut, lresid)) + geom_boxplot()

ggplot(diamonds2, aes(color, lresid)) + geom_boxplot()

ggplot(diamonds2, aes(clarity, lresid)) + geom_boxplot()

A more complicated model

We can add in effects due to color, cut, and clarity.

mod_diamond2 <- lm(lprice ~ lcarat + color + cut + clarity, data = diamonds2)

Create grid for new model. Use .model argument to data_grid() so predictors not specified are filled in with typical values. See documentation ?data_gird.

# Show output of data_grid() and add_predictions()
grid <- diamonds2 %>% 
  data_grid(cut, .model = mod_diamond2) %>% 
  add_predictions(mod_diamond2)
grid 

Visualize predictions (pred) of lprice versus cut:

grid %>%
    ggplot(aes(cut, pred)) +
    geom_point()

Check residuals of new model:

# Add residuals
diamonds2 <- diamonds2 %>% 
  add_residuals(mod_diamond2, "lresid2")
# Visualize residuals
ggplot(diamonds2, aes(lcarat, lresid2)) + 
  geom_hex(bins = 50)

Some diamonds have very large residuals. A residual of 2 indicates value is 2 ^ 2 = 4X our prediction. Can review diamonds with large residuals greater than abs(lresid2) > 1.

diamonds2 %>% 
    filter(abs(lresid2) > 1) %>% 
    add_predictions(mod_diamond2) %>% 
    mutate(pred = round(2 ^ pred)) %>% 
    select(price, pred, carat:table, x:z) %>% 
    arrange(price) 

Exercises

Exercise 4

Does the final model, mod_diamonds2, do a good job of predicting diamond prices? Would you trust it to tell you how much to spend if you were buying a diamond?

One way to determine the accuracy of the model is to examine the difference between the actual and predicted price. We’ll take a look at a histogram first.

# Add error calculation
diamonds2 <- diamonds2 %>% 
    add_predictions(mod_diamond2) %>% 
    mutate(pred = round(2 ^ pred),
           err = pred - price) 
# Visualize error
diamonds2 %>%
    ggplot(aes(err)) +
    geom_histogram(bins = 50) 

Th predictions appear to be centered around zero, but it’s difficult to see the range of errors. To better characterize the error, we can review the distribution using quantile().

probs <- c(0.005, 0.025, 0.25, 0.5, 0.75, 0.975, 0.995)
diamonds2$err %>% quantile(probs = probs)

     0.5%      2.5%       25%       50%       75%     97.5%     99.5% 
-2826.935 -1742.000  -195.000     1.000   156.000  1381.350  3196.480 

The model error (difference between actual and predicted price) is 95% of the time within [-$1742, $1381]. This is relatively high considering the median price is $2,396. Therefore, there is considerable variability in the predictions that should be weighed before using the model.

---

Flight Analysis

Question to answer: What affects the number of daily flights?

daily <- flights %>% 
    mutate(date = make_date(year, month, day)) %>% 
    group_by(date) %>% 
    summarise(n = n())
daily

ggplot(daily, aes(date, n)) + 
    geom_line()

Modelling the weekday effect

Problem: Understanding the trend is difficult due to the day of week variability.

Solution: Analyze the day of the week and remove the variability using a model.

daily <- daily %>% 
    mutate(wday = wday(date, label = TRUE))
ggplot(daily, aes(wday, n)) + 
    geom_boxplot()

Saturday is noticeably less due to less business travel.

We will apply a model to help reduce the variability.

# Model
mod <- lm(n ~ wday, data = daily)
# Create grid and add predictions
grid <- daily %>% 
    data_grid(wday) %>% 
    add_predictions(mod, "n")
# Visualize model
ggplot(daily, aes(wday, n)) + 
    geom_boxplot() +
    geom_point(data = grid, colour = "red", size = 4)

Next, visualize the residuals:

# add_residuals()
daily <- daily %>% 
    add_residuals(mod)
# Visualize residuals
daily %>% 
    ggplot(aes(date, resid)) + 
    geom_ref_line(h = 0) + 
    geom_line()

Day of the week effects have been removed, but there is still patterns that remain due to seasonality.

1. Viewing by weekday helps to show the large deviations. Saturday is noticeably different during summer (more flights than expected) and fall (less flights than expected).

ggplot(daily, aes(date, resid, colour = wday)) + 
    geom_ref_line(h = 0) + 
    geom_line()

2. Some days have far fewer flights than expected:

daily %>% filter(resid < -100)

3. There appears to be a long term trend over the course of the year.

daily %>% 
    ggplot(aes(date, resid)) + 
    geom_ref_line(h = 0) + 
    geom_line(colour = "grey50") + 
    geom_smooth(se = FALSE, span = 0.20)

Modelling the seasonal effect

Viewing Saturday highlights the seasonality effect:

daily %>% 
    filter(wday == "Sat") %>% 
    ggplot(aes(date, n)) + 
    geom_point() + 
    geom_line() +
    scale_x_date(NULL, date_breaks = "1 month", date_labels = "%b")

Create a term variable to split dates into spring, summer and fall seasons:

# Function to classify date as spring, summer, fall seasons
term <- function(date) {
    cut(date, 
        breaks = ymd(20130101, 20130605, 20130825, 20140101),
        labels = c("spring", "summer", "fall") 
    )
}
# Add season classification
daily <- daily %>% 
    mutate(term = term(date)) 
# Visualize season
daily %>% 
    filter(wday == "Sat") %>% 
    ggplot(aes(date, n, colour = term)) +
    geom_point(alpha = 1/3) + 
    geom_line() +
    scale_x_date(NULL, date_breaks = "1 month", date_labels = "%b")

We can visualize the term by day of the week:

daily %>% 
    ggplot(aes(wday, n, colour = term)) +
    geom_boxplot()

We can see summer tends to be higher, while fall and spring tend to be about the same with the exception of Saturday.

We can model using an interaction effect between weekday (wday) and season (term). Use gather_residuals() to compute residuals on both models for comparison.

# Models to compare
mod1 <- lm(n ~ wday, data = daily)
mod2 <- lm(n ~ wday * term, data = daily)
# gather_residuals(), visualize both models
daily %>% 
    gather_residuals(without_term = mod1, with_term = mod2) %>% 
    ggplot(aes(date, resid, colour = model)) +
    geom_line(alpha = 0.75)

Some improvement in summer, but not as much as we want. Visualizing predictions by term and day shows an issue.

# data_grid() by wday and term, add_predictions()
grid <- daily %>% 
    data_grid(wday, term) %>% 
    add_predictions(mod2, "n")
# Visualize predictions by term and wday
ggplot(daily, aes(wday, n)) +
    geom_boxplot() + 
    geom_point(data = grid, colour = "red") + 
    facet_wrap(~ term)

Problem: Our model is finding the mean effect, but large outliers are throwing it off.

Solution: Use a modelling algorithm that is affected less by outliers. MASS::rlm() is resistant to outliers.

# Model using rlm()
mod3 <- MASS::rlm(n ~ wday * term, data = daily)
# Visualize residuals
daily %>% 
    add_residuals(mod3, "resid") %>% 
    ggplot(aes(date, resid)) + 
    geom_hline(yintercept = 0, size = 2, colour = "white") + 
    geom_line()

Alternative approach

We can use splines::ns() to model by weekday throughout the year.

library(splines)
# Model using rlm() with 5th order polynomial for date
mod <- MASS::rlm(n ~ wday * ns(date, 5), data = daily)
# data_grid() by wday along 13 evenly spaced points, add_predictions() 
# and visualize
daily %>% 
    data_grid(wday, date = seq_range(date, n = 13)) %>% 
    add_predictions(mod) %>% 
    ggplot(aes(date, pred, colour = wday)) + 
    geom_line() +
    geom_point()

We can also check the residuals to see how this approach compares to the previous model.

daily %>% 
    add_residuals(mod) %>% 
    ggplot(aes(date, resid)) + 
    geom_hline(yintercept = 0, size = 2, colour = "white") + 
    geom_line()