course
/
r4ds
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Restructuring model chapters

This commit is contained in:
hadley 2016-06-13 09:50:55 -05:00
parent f84f26a7ec
commit ae14d3a3cc
7 changed files with 459 additions and 1090 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
_bookdown.yml
View File

@ -25,7 +25,8 @@ rmd_files: [
"model.Rmd",
"model-basics.Rmd",
"model-vis.Rmd",
"model-building.Rmd",
"model-many.Rmd",
"model-assess.Rmd",
"communicate.Rmd",

10
model-assess.Rmd
View File

@ -10,16 +10,6 @@ options(digits = 3)
* Bootstrapping to understand uncertainty in parameters.
* Cross-validation to understand predictive quality.
## Purrr + dplyr
i.e. how do dplyr and purrr intersect.
* Why use a data frame?
* List columns in a data frame
* Mutate & filter.
* Creating list columns with `group_by()` and `do()`.
## Multiple models
A natural application of `map2()` is handling test-training pairs when doing model evaluation. This is an important modelling technique: you should never evaluate a model on the same data it was fit to because it's going to make you overconfident. Instead, it's better to divide the data up and use one piece to fit the model and the other piece to evaluate it. A popular technique for this is called k-fold cross validation. You randomly hold out x% of the data and fit the model to the rest. You need to repeat this a few times because of random variation.

62
model-basics.Rmd
View File

@ -2,9 +2,6 @@
A model is a function that summarizes how the values of one variable vary in relation to the values of other variables. Models play a large role in hypothesis testing and prediction, but for the moment you should think of models just like you think of statistics. A statistic summarizes a *distribution* in a way that is easy to understand; and a model summarizes *covariation* in a way that is easy to understand. In other words, a model is just another way to describe data.
In this book, we'll focus on deep tools for rather simple models. We want to give you tools to help you build your intuition for what models do. There's little mathematical formalism, and only a relatively small overlap with how modelling is normally presented in statistics. Modelling is a huge topic and we can only scratch the surface here. You'll definitely need to consult other resources as the complexity of your models grow.
We're going to give you a basic strategy, and point you to places to learn more. The key is to think about data generated from your model as regular data - you're going to want to manipulate it and visualise it in many different ways. Being good at modelling is a mixture of having some good general principles and having a big toolbox of techniques. Here we'll focus on general techniques to help you undertand what your model is telling you.
This chapter will explain how to build useful models with R.
@ -406,6 +403,28 @@ hes2 <- lm(income ~ height + education * sex, data = heights)
tidy(hes2)
```
### Nested variables
Another case that occassionally crops up is nested variables: you have an identifier that is locally unique, not globally unique. For example you might have this data about students in schools:
```{r}
students <- tibble::frame_data(
~student_id, ~school_id,
1, 1,
2, 1,
1, 2,
1, 3,
2, 3,
3, 3
)
```
The student id only makes sense in the context of the school: it doesn't make sense to generate every combination of student and school. You can use `nesting()` for this case:
```{r}
students %>% expand(nesting(school_id, student_id))
```
### Continuous
There appears to be a relationship between a person's education and how poorly the model predicts their income. If we graph the model residuals against `education` above, we see that the more a person is educated, the worse the model underestimates their income:
@ -545,6 +564,43 @@ data_frame(education = seq(5, 25)) %>%
geom_line()
```
Other useful arguments to `seq_range()`:
* `pretty = TRUE` will generate a "pretty" sequence, i.e. something that looks
nice to the human eye. This is useful if you want to produce tables of
output:
```{r}
seq_range(c(0.0123, 0.923423), n = 5)
seq_range(c(0.0123, 0.923423), n = 5, pretty = TRUE)
```
* `trim = 0.1` will trim off 10% of the tail values. This is useful if the
variables has an long tailed distribution and you want to focus on generating
values near the center:
```{r}
x <- rcauchy(100)
seq_range(x, n = 5)
seq_range(x, n = 5, trim = 0.10)
seq_range(x, n = 5, trim = 0.25)
seq_range(x, n = 5, trim = 0.50)
```
### Interpolation vs extrapolation
One danger with prediction plots is that it's easy to make predictions that are far away from the original data. This is dangerous because it's quite possible that the model (which is a simplification of reality) will no longer apply far away from observed values.
As the number of variables in your model grows ... "the curse of dimensionality": as the number of variables increases the average distance between points increases. That means most of the space is very sparse, and you have to rely on strong assumptions.
To help avoid this problem, it's good practice to include "nearby" observed data points in any prediction plot. These help you see if you're interpolating, making prediction "in between" existing data points, or extrapolating, making predictions about preivously unobserved slices of the data.
One way to do this is to use `condvis::visualweight()`.
<https://cran.rstudio.com/web/packages/condvis/>
## Response variations
### Transformations

357
model-building.Rmd Normal file
View File

@ -0,0 +1,357 @@
```{r setup, include = FALSE}
library(broom)
library(ggplot2)
library(dplyr)
library(lubridate)
library(tidyr)
library(nycflights13)
library(modelr)
```
# Model building
In the previous chapter you learned how some basic models worked, and learned some basic tools for understanding what a model is telling you about your data. In this chapter, we're going talk more about the model building process: how you start from nothing, and end up with a good model.
* Regression modelling strategies.
The way we're going to work is to subtract patterns from the data, while adding them to 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 is going to be a lot of work. There are certainly alternative approaches - a more machine learning approach is simply to focus on improving the predictive ability of the model, being careful to fairly assess it (i.e. not assessing the model on the data that was used to train it). These approaches tend to produce black boxes - i.e. the model does a really good job, but you don't know why. This is fine, but the main problem is that you can't apply your real world knowledge to the model to think about whether or not it's likely to 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 ML model building approach. If prediction is important, get to a good point, and then use visulisation to understand the most important parts of the model.
> 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."
-- Reddit user Broseidon241, https://www.reddit.com/r/datascience/comments/4irajq/mistakes_made_by_beginningaspiring_data_scientists/
### Prerequisites
```{r}
library(modelr)
options(
contrasts = c("contr.treatment", "contr.treatment"),
na.option = na.exclude
)
library(nycflights13)
library(lubridate)
library(dplyr)
```
## What affects the number of daily flights?
We're going to start by building a model to help us understand the number of flights per day that leave NYC. 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.
We'll start by using dplyr to generate the data of interest:
```{r}
daily <- flights %>%
mutate(date = make_datetime(year, month, day)) %>%
group_by(date) %>%
summarise(n = n())
daily
```
And then we'll plot it with ggplot2:
```{r}
ggplot(daily, aes(date, n)) +
geom_line()
```
### Day of week
This is a small dataset, but there's a lot of pattern to explore. Understanding the long-term trend is challenging because there's a very strong day-of-week effect that dominates the subtler patterns:
```{r}
daily <- daily %>%
mutate(wday = wday(date, label = TRUE))
ggplot(daily, aes(wday, n)) +
geom_boxplot()
```
There are fewer flights on weekends because a very large proportion of travel is for business. You might sometimes have to less on Sunday for an early flight, but it's very rare that you'd leave on Saturday: you'd much rather be home with your family.
One way to remove this strong pattern is to fit a model that "explains" (i.e. attempts to predict) the day of week effect, and then look at the residuals. Another way of thinking about this is that we're capturing the day-of-week effect, moving it from the data, into a model.
```{r}
mod <- lm(n ~ wday, data = daily)
daily <- daily %>% add_residuals(n_resid = mod)
daily %>%
ggplot(aes(date, n_resid)) +
geom_hline(yintercept = 0, size = 2, colour = "white") +
geom_line()
daily %>%
expand(wday) %>%
add_predictions(n_pred = mod) %>%
ggplot(aes(wday, n_pred)) +
geom_point()
```
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 interesting 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, n_resid, colour = wday)) +
geom_hline(yintercept = 0, size = 2, colour = "white") +
geom_line()
```
The problem appears to be Saturdays: it seems like during summer there are
more flights on Saturdays than we expect, and during Fall there are fewer.
I suspect this is because of summer holidays: many people go on holiday
in the summer, and people don't mind travelling on Saturdays for vacation.
(This doesn't, however, explain why there are more Satruday flights in
spring than fall).
1. There are some days with far fewer flights than expected:
```{r}
daily %>% filter(n_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, n_resid)) +
geom_hline(yintercept = 0, size = 2, colour = "white") +
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 more with this trend than brainstorm possible
explanations because we only have a single year's worth of data.
### Seasonal Saturday effect
We'll tackle the day of week effect first. Let's zoom in on Saturdays, going back to raw numbers:
```{r}
daily %>%
filter(wday == "Sat") %>%
ggplot(aes(date, n)) +
geom_line() +
geom_point(alpha = 1/3) +
scale_x_datetime(NULL, date_breaks = "1 month", date_labels = "%b")
```
So it looks like 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 is Jun 26 - Sep 9. Few families travel in the fall because of the big Thanksgiving and Christmas holidays. So lets add a "term" variable to attemp to control for that.
```{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}
mod2 <- lm(n ~ wday * term, data = daily)
daily$n_resid2 <- resid(mod2)
daily %>%
gather(model, prediction, n_resid, n_resid2) %>%
ggplot(aes(date, prediction, colour = model)) +
geom_line(alpha = 0.75)
```
That's because this model is basically calculating an average for each combination of wday and school term. (How many observations do we have for each day of week in each term?) 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: `rlm`. This greatly reduces the impact of the outliers on our estimates, and gives a result that does a good job of removing the day of week pattern:
```{r, warn = FALSE}
mod3 <- MASS::rlm(n ~ wday * term, data = daily)
daily <- daily %>% add_residuals(n_resid3 = mod3)
ggplot(daily, aes(date, n_resid3)) +
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.
Very common to use residual plots when figuring out if a model is ok. But it's easy to get the impression that there's just one type of residual plot you should do, when in fact there are infinite.
### Time of year: an alternative approach
In the previous section we used our knowledge of phenomonen 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.
When you have a continuous variable in the model, rather than using the unique values that you've seen, it's often more useful to generate an evenly spaced grid. One convenient way to do this is with `modelr::seq_range()` which takes a continuous variable, calculates its range, and then generates an evenly spaced points between the minimum and maximum.
```{r, warn = FALSE}
mod <- MASS::rlm(n ~ wday * yday(date), data = daily)
grid <- daily %>%
tidyr::expand(wday, date = seq_range(date, n = 13)) %>%
add_predictions(mod = mod)
ggplot(grid, aes(date, mod, colour = wday)) +
geom_line() +
geom_point()
```
(Why use `yday(date)` instead of `date`? That's saying we think that the pattern depends only the day of the year, so we might expect it to be the same in other years. Again, since we only have a single year of data we can't test that hypothesis.)
We know that this pattern doesn't do a good job of capturing the variation in the data. There isn't a simple linear trend across the entire year, 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 = mod) %>%
ggplot(aes(date, mod, colour = wday)) +
geom_line() +
geom_point()
```
Particularly, we see the strongly pattern in Saturdays that we identified when coming in the opposite direction. It's always a good sign when you see the same signal from multiple approaches. (But note our previous model was explanatory - this is just predictatory.)
How many degrees of freedom to use? Either pick manually to extract the shape of the data, or you can use one of the model assessment techniques in the following chapter to pick algorithmically. Here we're most interested in explanation, so picking by hand (with a little though and plenty of scepticism) is typically fine.
### Public holidays
### 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.
```{r}
compute_vars <- function(data) {
data %>% mutate(
term = term(date),
wday = wday(date, label = TRUE)
)
}
```
Another option is to wrap it ito the model formula:
```{r}
wday2 <- function(x) wday(x, label = TRUE)
mod3 <- lm(n ~ wday2(date) * term(date), data = daily)
daily %>%
expand(date) %>%
add_predictions(pred = mod3)
```
I think this is fine to do provided that you've carefully checked that the functions do what you think they do (i.e. with a visualisation). There are two disadvantages:
1. You may need to add the variables back in anyway if you want to use
them in a visualsiation.
1. When looking at the coefficients, their values are longer and harder to
read. (But this is a general problem with the way that linear models report
categorical coefficients in R, not a specific problem with this case.)
### 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(n_resid2 > 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?
Why is this not very helpful?
1. Above we made the hypothesis 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:
if it's true, you'd expect to see more Sunday 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 levels of the factor so
that the week starts on Monday.
1. Compare the predictions for each `wday` combined with `term` for the
`lm` and `rlm`
## Case study: predicting flight delays
Finish off with a somewhat more realistic case study where we combine the techniques of visualising predictions and residuals to attack the problem of predicting flight delays.
Can't predict delays for next year. Why not? Instead we'll focus on predicting the amount that your flight will be delayed if it's leaving soon.
To tackle a problem like this, it's worthwhile to start with some brainstorming. More ideas will come to you as you proceed but it's a good idea to start by seeding the pot
* time of day
* weather
* plane?
* number of flights
We'll start with some exploratory analysis, and then work on the model:
```{r}
delays <- flights %>%
mutate(date = make_datetime(year, month, day)) %>%
group_by(date) %>%
summarise(delay = mean(arr_delay, na.rm = TRUE), cancelled = mean(is.na(dep_time)), n = n())
# delays %>%
# ggplot(aes(wday(date, label = TRUE), delay)) +
# geom_boxplot()
delays %>%
ggplot(aes(n, delay)) +
geom_point() +
geom_smooth(se = F)
```

14
model-many.Rmd Normal file
View File

@ -0,0 +1,14 @@
# Working with many models
```{r setup-model, include=FALSE}
library(purrr)
set.seed(1014)
options(digits = 3)
```
i.e. how do tidyr, dplyr and purrr intersect.
* Why use a data frame?
* List columns in a data frame
* Mutate & filter.
* Creating list columns with `group_by()` and `do()`.

1076
model-vis.Rmd
View File

File diff suppressed because it is too large Load Diff

27
model.Rmd
View File

@ -1,5 +1,32 @@
```{r include=FALSE, cache=FALSE}
set.seed(1014)
options(digits = 3)
knitr::opts_chunk$set(
comment = "#>",
collapse = TRUE,
cache = TRUE
)
options(dplyr.print_min = 6, dplyr.print_max = 6)
```
# (PART) Model {-}
# Introduction
The scientific method guides data science. Data science solves known problems with the scientific method.
In this book, we'll focus on deep tools for rather simple models. We want to give you tools to help you build your intuition for what models do. There's little mathematical formalism, and only a relatively small overlap with how modelling is normally presented in statistics. Modelling is a huge topic and we can only scratch the surface here. You'll definitely need to consult other resources as the complexity of your models grow.
We're going to give you a basic strategy, and point you to places to learn more. The key is to think about data generated from your model as regular data - you're going to want to manipulate it and visualise it in many different ways. Being good at modelling is a mixture of having some good general principles and having a big toolbox of techniques. Here we'll focus on general techniques to help you undertand what your model is telling you.
In the course of modelling, you'll often discover data quality problems. Maybe a missing value is recorded as 999. Whenever you discover a problem like this, you'll need to review an update your import scripts. You'll often discover a problem with one variable, but you'll need to think about it for all variables. This is often frustrating, but it's typical.
<https://blog.engineyard.com/2014/pets-vs-cattle>.
<https://en.wikipedia.org/wiki/R/K_selection_theory>
In this book we are going to focus on models primarily as tools for description. This is rather non-standard because we're normally interested in models for their inferential power: their ability to make accurate predictions for observations that we haven't seen yet.
In other words, in this book, we're typically going to think about a good model as a model that well captures the patterns that we see in the data. For now, a good model captures the majority of the patterns that are generated by the underlying mechanism of interest, and captures few patterns that are not generated by that mechanism. When you go on from this book and learn other ways of thinking about models this will stand you in good stead: if you can't capture patterns in the data that you can see, it's unlikely you'll be able to make good predictions about data that you haven't seen.
Generally, however, this will tend to make us over-optimistic about the quality of our model. Chapter XXX you'll start to learn more about how we can judge the quality of a model on data that it was 't fit it. But you have to beware of overfitting the data - in the next section we'll discuss some formal methods. But a healthy dose of scepticism is also as powerful as precise quantitative methods: do you believe that a pattern you see in your sample is going to generalise to a wider population?