r4ds/model-building.Rmd

# Model building

## Introduction

In the previous chapter you learned how linear models worked, and learned some basic tools for understanding what a model is telling you about your data. The previous chapter focussed on simulated datasets to help you learn about how models work. This chapter will focus on real data, showing you how you can progressively build up a model to aid your understanding of the data.

We will take advantage of the fact that you can think about a model partition your data into pattern and residuals. We'll find patterns with visualisation, then make the concrete and precise with a model. We'll them repeat the process, replace the old response variable with the residuals from the model. The goal is to transition from implicit knowledge in the data and your head to explicit knowledge in a quantitative model. This makes it easier to apply to new domains, and easier for others to use. 

For very large and complex datasets this will a lot of work. There are certainly alternative approaches - a more machine learning approach is simply to focus on the predictive ability of the model. These approaches tend to produce black boxes: the model does a really good job at generating predictions, but you don't know why. This is a totally reasonable approach, but it does make it hard to apply your real world knowledge to the model. That, in turn, makes it difficult to assess whether or not the model will continue work in the long-term, as fundamentals change. For most real models, I'd expect you to use some combination of this approach and a more classic automated approach.

It's a challenge to know when to stop. You need to figure out when your model is good enough, and when additional investment is unlikely to pay off. I particularly this quote from reddit user Broseidon241: 

> A long time ago in art class, my teacher told me "An artist needs to know 
> when a piece is done. You can't tweak something into perfection - wrap it up. 
> If you don't like it, do it over again. Otherwise begin something new". Later
> in life, I heard "A poor seamstress makes many mistakes. A good seamstress 
> works hard to correct those mistakes. A great seamstress isn't afraid to 
> throw out the garment and start over."

-- Broseidon241, <https://www.reddit.com/r/datascience/comments/4irajq>

### Prerequisites

We'll start with modelling and EDA tools we needed in the last chapter. Then we'll add in some real datasets: `diamonds` from ggplot2, and `flights` from nycflights13. We'll also need lubridate to extract useful components of datetimes.

```{r setup, message = FALSE}
# Modelling functions
library(modelr)
library(broom)

# EDA tools
library(ggplot2)
library(dplyr)
library(tidyr)

# Data
library(nycflights13)
library(lubridate)
```

## Why are low quality diamonds more expensive?

In previous chapters we've seen a surprising relationship between the quality of diamonds and their price: low quality diamonds (poor cuts, bad colours, and inferior clarity) have higher prices.

```{r dev = "png"}
ggplot(diamonds, aes(cut, price)) + geom_boxplot()
ggplot(diamonds, aes(color, price)) + geom_boxplot()
ggplot(diamonds, aes(clarity, price)) + geom_boxplot()
```

### Price and carat

The basic reason we see this pattern is that the variable that is most predictive of diamond price is its size:

```{r, dev = "png"}
ggplot(diamonds, aes(carat, price)) + 
  geom_hex(bins = 50)
```

To explore this relationship, lets make a couple of tweaks to the diamonds dataset:

1. Remove all diamonds bigger than 2.5 carats
1. Log-transform the carat and price variables

```{r}
diamonds2 <- diamonds %>% 
  filter(carat <= 2.5) %>% 
  mutate(lprice = log2(price), lcarat = log2(carat))

ggplot(diamonds2, aes(lcarat, lprice)) + 
  geom_hex(bins = 50)
```

Log-transforming is very useful here because it make a linear relationship, and linear relationships are generally much easier to work with.  We could go one step further and use a linear model to remove the strong effect of `lcarat` on `lprice`. First, let's see what a linear model tells about the data on the original scale:

```{r}
mod_diamond <- lm(lprice ~ lcarat, data = diamonds2)

grid <- diamonds2 %>% 
  expand(carat = seq_range(carat, 20)) %>% 
  mutate(lcarat = log2(carat)) %>% 
  add_predictions(mod_diamond, "lprice") %>% 
  mutate(price = 2 ^ lprice)

ggplot(diamonds2, aes(carat, price)) + 
  geom_hex(bins = 50) + 
  geom_line(data = grid, colour = "red", size = 1)
```

That's interesting! If we believe our model, then it suggests that the large diamonds we have are much cheaper than expected. This may because no diamond in this dataset costs more than $19,000.

We can also look at the residuals from this model. This verifies that we have successfully removed the strong linear pattern:

```{r}
diamonds2 <- diamonds2 %>% 
  add_residuals(mod_diamond, "lresid")

ggplot(diamonds2, aes(lcarat, lresid)) + 
  geom_hex(bins = 50)
```

Importantly, we can now use those residuals in plots instead of `price`. 

```{r dev = "png"}
ggplot(diamonds2, aes(cut, lresid)) + geom_boxplot()
ggplot(diamonds2, aes(color, lresid)) + geom_boxplot()
ggplot(diamonds2, aes(clarity, lresid)) + geom_boxplot()
```

Here we see the relationship we'd expect. Now that we've removed the effect of size on price,  

To interpret the `y` axis, we need to think about what the residuals are telling us, and what scale they are on. A residual of -1 indicates that `lprice` was 1 unit lower than expected, based on the `carat` alone. $2^{-1}$ is 1/2, so that sugggests diamonds with colour I1 are half the price you'd expect.

### A model complicated model

We could continue this process, making our model complex:

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

add_predictions_trans <- function(df, mod) {
  df %>% 
    add_predictions(mod, "lpred") %>% 
    mutate(pred = 2 ^ lpred)
}

diamonds2 %>% 
  expand_model(mod_diamond2, cut) %>% 
  add_predictions_trans(mod_diamond2) %>% 
  ggplot(aes(cut, pred)) + 
    geom_point()
```


### Exercises

1.  In the plot of `lcarat` vs. `lprice`, there are some bright vertical
    strips. What do they represent?

1.  If `log(price) = a_0 + a_1 * log(carat)`, what does that say about 
    the relationship between `price` and `carat?
    
1.  Extract the diamonds that have very high and very low residuals. 
    Is there any unusual about these diamonds? Are the particularly bad
    or good, or do you think these are pricing errors?

## What affects the number of daily flights?

Let's explore the number of flights that leave NYC per day. We're not going to end up with a fully realised model, but as you'll see, the steps along the way will help us better understand the data. Let's get started by counting the number of flights per day and visualising it with ggplot2.

```{r}
daily <- flights %>% 
  mutate(date = make_datetime(year, month, day)) %>% 
  group_by(date) %>% 
  summarise(n = n())
daily

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

This is a really small dataset --- only 365 rows and 2 columns, but because as you'll see there's a rich set of interesting variables buried in the date.

### Day of week

Understanding the long-term trend is challenging because there's a very strong day-of-week effect that dominates the subtler patterns. Let's summarise the number of flights per day-of-week:

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

There are fewer flights on weekends because most travel is for business. The effect is particularly pronounced on Saturday: you might sometimes have to leave on Sunday for a Monday morning meeting, but it's very rare that you'd leave on Saturday as you'd much rather be at home with your family.

One way to remove this strong pattern is to use a model. First, we fit the model, and display its predictions overlaid on the original data:

```{r}
mod <- lm(n ~ wday, data = daily)

grid <- daily %>% 
  expand(wday) %>% 
  add_predictions(mod, "n")

ggplot(daily, aes(wday, n)) + 
  geom_boxplot() +
  geom_point(data = grid, colour = "red", size = 4)
```

Next we compute and visualise the residuals:

```{r}
daily <- daily %>% add_residuals(mod)
daily %>% 
  ggplot(aes(date, resid)) + 
  geom_ref_line(h = 0) + 
  geom_line()
```

Note the change in the y-axis: now we are seeing the deviation from the expected number of flights, given the day of week. This plot is useful because now that we've removed much of the large day-of-week effect, we can see some of the subtler patterns that remain:

1.  Our day of week adjustment seems to fail starting around June: you can
    still see a strong regular pattern that our model hasn't removed. Drawing
    a plot with one line for each day of the week makes the cause easier
    to see:

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

    Our model fails to accurately predict the number of flights on Saturday:
    during summer there are more flights than we expect, and during Fall there 
    are fewer. We'll see how we can do better in the next section.

1.  There are some days with far fewer flights than expected:

    ```{r}
    daily %>% filter(resid < -100)
    ```

    If you're familiar with American public holidays, you might spot New Year's 
    day, July 4th, Thanksgiving and Christmas. There are some others that don't 
    seem to correspond immediately to public holidays. You'll work on those 
    in the exercise below.
    
1.  There seems to be some smoother long term trend over the course of a year.
    We can highlight that trend with `geom_smooth()`:

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

    There are fewer flights in January (and December), and more in summer 
    (May-Sep). We can't do much with this pattern numerically, because we only
    have a single year of data. But we can use our domain knowledge to 
    brainstorm potential explanations.

### Seasonal Saturday effect

Let's first tackle our failure to accurately predict the number of flights on Saturday. A good place to start is to go back to the raw numbers, focussing on Saturdays:

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

(I've used both points and lines to make it more clear what is data and what is interpolation.)

I suspect pattern is caused by summer holidays: many people go on holiday in the summer, and people don't mind travelling on Saturdays for vacation. Looking at this plot, we might guess that summer holidays are from early June to late August. That seems to line up fairly well with the [state's school terms](http://schools.nyc.gov/Calendar/2013-2014+School+Year+Calendars.htm): summer break in 2013 was Jun 26--Sep 9. 

Why are there Saturday flights in the Fall than the Spring? I asked some American friends and they suggested that it's less common to plan family vacations during the Fall becuase of the big Thanksgiving and Christmas holidays. We can't tell if that's exactly the reason, but it seems like a plausible working hypothesis.

Lets create a "term" variable that roughly captures the three school terms, and check our work with a plot:

```{r}
term <- function(date) {
  cut(date, 
    breaks = as.POSIXct(ymd(20130101, 20130605, 20130825, 20140101)),
    labels = c("spring", "summer", "fall") 
  )
}

daily <- daily %>% mutate(term = term(date)) 

daily %>% 
  filter(wday == "Sat") %>% 
  ggplot(aes(date, n, colour = term)) +
  geom_point(alpha = 1/3) + 
  geom_line() +
  scale_x_datetime(NULL, date_breaks = "1 month", date_labels = "%b")
```

(I manually tweaked the dates to get nice breaks in the plot. Using a visualisation to help you understand what your function is doing is a really powerful and general technique.)

It's useful to see how this new variable affects the other days of the week:

```{r}
daily %>% 
  ggplot(aes(wday, n, colour = term)) +
    geom_boxplot()
```

It looks like there is significant variation across the terms, so fitting a separate day of week effect for each term is reasonable. This improves our model, but not as much as we might hope:

```{r}
mod1 <- lm(n ~ wday, data = daily)
mod2 <- lm(n ~ wday * term, data = daily)

daily %>% 
  gather_residuals(without_term = mod1, with_term = mod2) %>% 
  ggplot(aes(date, resid, colour = model)) +
    geom_line(alpha = 0.75)
```

We can see the problem by overlaying the predictions from the model on to the raw data:

```{r}
grid <- daily %>% 
  expand(wday, term) %>% 
  add_predictions(mod2, "n")

ggplot(daily, aes(wday, n)) +
  geom_boxplot() + 
  geom_point(data = grid, colour = "red") + 
  facet_wrap(~ term)
```

Our model is finding the _mean_ effect, but we have a lot of big outliers, so they tend to drag the mean far away from the typical value. We can alleviate this problem by using a model that is robust to the effect of outliers: `MASS::rlm()`. This greatly reduces the impact of the outliers on our estimates, and gives a model that does a good job of removing the day of week pattern:

```{r, warn = FALSE}
mod3 <- MASS::rlm(n ~ wday * term, data = daily)

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

It's now much easier to see the long-term trend, and the positive and negative outliers.

### Time of year: an alternative approach

In the previous section we used our domain knowledge (how the US school term affects travel) to improve the model. An alternative to using making our knowledge explicit in the model is to give the data more room to speak. We could use a more flexible model and allow that to capture the pattern we're interested in. We know that a simple linear trend isn't adeqaute, so instead we could use a natural spline to allow a smoothly varying trend across the year:

```{r}
library(splines)
mod <- MASS::rlm(n ~ wday * ns(date, 5), data = daily)

daily %>% 
  tidyr::expand(wday, date = seq_range(date, n = 13)) %>% 
  add_predictions(mod) %>% 
  ggplot(aes(date, pred, colour = wday)) + 
    geom_line() +
    geom_point()
```

We see a strong pattern in the numbers of Saturday flights. This is reassuring, because we also saw that pattern in the raw data. It's  a good sign when you see the same signal from multiple approaches.

How do you decide how many parameters to use for the spline? You can either either it pick by eye, or you could use automated techniques which you'll learn about in [model assessment]. For exploration, picking by eye to capture the most important patterns is fine.

### Computed variables

If you're experimenting with many models and many visualisations, it's a good idea to bundle the creation of variables up into a function so there's no chance of accidentally applying a different transformation in different places. For example, we could write:

```{r}
compute_vars <- function(data) {
  data %>% mutate(
    term = term(date), 
    wday = wday(date, label = TRUE)
  )
}
```

Another option is to put the transformations directly in the model formula:

```{r}
wday2 <- function(x) wday(x, label = TRUE)
mod3 <- lm(n ~ wday2(date) * term(date), data = daily)
```

Either approach is reasonable. Making the transformed variable explicit is useful if you want to check your work, or use them in a visualisation. But you can't easily use transformations (like splines) that return multiple columns. Including the transformations in the model function makes life a little easier when you're working with many different datasets because the model is self contained.

### Exercises

1.  Use your google sleuthing skills to brainstorm why there were fewer than
    expected flights on Jan 20, May 26, and Sep 9. (Hint: they all have the
    same explanation.) How would these days generalise to another year?

1.  What do the three days with high positive residuals represent?
    How would these days generalise to another year?

    ```{r}
    daily %>% filter(resid > 80)
    ```

1.  Create a new variable that splits the `wday` variable into terms, but only
    for Saturdays, i.e. it should have `Thurs`, `Fri`, but `Sat-summer`, 
    `Sat-spring`, `Sat-fall`. How does this model compare with the model with 
    every combination of `wday` and `term`?
    
1.  Create a new wday variable that combines the day of week, term 
    (for Saturdays), and public holidays. What do the residuals of 
    that model look like?

1.  What happens if you fit a day of week effect that varies by month 
    (i.e. `n ~ wday * month`)? Why is this not very helpful? 

1.  What would you expect the model `n ~ wday + ns(date, 5)` to look like?
    Knowing what you know about the data, why would you expect it to be
    not particularly effective?

1.  We hypothesised that people leaving on Sundays are more likely to be 
    business travellers who need to be somewhere on Monday. Explore that 
    hypothesis by seeing how it breaks down based on distance and itme: if 
    it's true, you'd expect to see more Sunday evening flights to places that 
    are far away.

1.  It's a little frustrating that Sunday and Saturday are on separate ends
    of the plot. Write a small function to set the manipulate the levels of the 
    factor so that the week starts on Monday.