Chapter 25: Many Models

Overview

This is a reproduction of the workflow for creating many models, which comes from Chapter 25 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 Gapminder Analysis illustrates a workflow for integrating modelling to multiple subsets of data using the modelr and broom packages. In addition, the purrr mapping functions are integrated into the analysis workflow. The ultimate goal is to showcase the concepts that enable analysts to visualize the models, extracting meaningful information on residuals, r-squared values, etc.

The details behind creating many models cover both creating and simplifying nested data frames. These skills are necessary for the modelling workflow.

---

Prerequisites

library(modelr)
options(na.action = na.warn)
library(broom)
library(tidyverse) # ggplot2, purrr, dplyr, tidyr, readr, tibble
library(gridExtra) # grid.arange()
library(gapminder) # gapminder data

---

Gapminder Analysis: Multiple Models Workflow

gapminder 

Question to answer: “How does life expectancy (lifeExp) change over time (year) for each country (country)?”

gapminder %>%
    ggplot(aes(year, lifeExp, group = country)) + 
    geom_line(alpha = 0.35)

Problems:

• Difficult to see what’s going on.
• General trend is positive, but some countries don’t follow the group pattern.
• Need to be able to make countries not following pattern easier to stand out.

Solution: Use models. Model captures trend and residuals capture remaining pattern (if any).

Process for Single Country:

nz <- filter(gapminder, country == "New Zealand")
p1 <- nz %>%
    ggplot(aes(year, lifeExp)) + 
    geom_line() + 
    ggtitle("Full data = ")
nz_mod <- lm(lifeExp ~ year, data = nz) # Model
p2 <- nz %>%
    add_predictions(nz_mod) %>%
    ggplot(aes(year, pred)) +
    geom_line() + 
    ggtitle("Linear trend + ")
p3 <- nz %>%
    add_residuals(nz_mod) %>%
    ggplot(aes(year, resid)) +
    geom_line() + 
    ggtitle("Remaining pattern")
    
grid.arrange(p1, p2, p3, ncol=3)

Nested data (nest)

Now we need to extend this analysis to many other countries. We do this using a map function from purrr by extracting out a common function and using map to apply to all elements of a list or data frame.

Problem: Data is structured differently. We need to repeat the map function one a subset of rows instead of each single row.

Solution: Nested data frame. Easiest way is to use group_by() then pass to nest().

by_country <- gapminder %>%
    group_by(country, continent) %>%
    nest()
by_country %>% head 

Can drill into the data column.

by_country$data[[1]] %>% head()

List-columns (use map on nested columns)

Next, need a model function to iterate over.

country_model <- function(df) {
    lm(lifeExp ~ year, data = df)
}

Then apply to each country using purrr::map(). The workflow stores the model as part of the nested data frame.

by_country <- by_country %>%
    mutate(model = map(.x = data, .f = country_model))
by_country %>% head()

Follow the same process to add predictions and residuals.

by_country <- by_country %>%
    mutate(
        preds  = map2(.x = data, .y = model, .f = add_predictions),
        resids = map2(.x = data, .y = model, .f = add_residuals)
    )
by_country %>% head()

Unnesting (unnest)

Now we need to plot, but this is very difficult in current structure.

Problem: The nested data frame structure makes plotting difficult. Need to be able to plot, but structure is a list of data frames.

Solution: Unnest with unnest() function.

preds <- unnest(data = by_country, preds)
preds %>% head()

resids <- unnest(data = by_country, resids)
resids %>% head()

First, take a look at trends (preds):

preds %>%
    ggplot(aes(year, pred, group = country)) +
    geom_line(alpha = 0.35) +
    facet_wrap(~ continent)

Second, take a look at resids (resids):

resids %>%
    ggplot(aes(year, resid, group = country)) +
    geom_line(alpha = 0.35) + 
    facet_wrap(~ continent)

Upon inspection of residuals, we can get the following takeaways:

1. It looks like we missed some mild pattern.

2. Africa has a weird wave motion in the residuals indicating something interesting that needs further investigation.

Assessing model quality (broom package)

Quick overview of broom package

The broom package provides three general tools for turning models into tidy data frames:

1. broom::glance(model) returns a row for each model. Each column gives a model summary: either a measure of model quality, or complexity, or a combination of the two.

2. broom::tidy(model) returns a row for each coefficient in the model. Each column gives information about the estimate or its variability.

3. broom::augment(model, data) returns a row for each row in data, adding extra values like residuals, and influence statistics.

glance(nz_mod)

tidy(nz_mod)

augment(nz_mod, nz) %>% head()

Back to assessing the model quality

We want to use broom::glance() to get the lowest r.squared values. Two steps:

1. Use map() to apply glance() to the models.

2. Use unnest() to return the model statistics. Set .drop = TRUE to suppress unnested list-columns.

glance <- by_country %>% 
    mutate(glance = map(model, broom::glance)) %>% 
    unnest(glance, .drop = TRUE) %>%
    arrange(r.squared) 
glance %>% head()

Visualize the r-squared values:

Africa has some low r-squared values.

glance %>%
    ggplot(aes(continent, r.squared)) + 
    geom_jitter()

Review countries with low r-squared values:

bad_fit <- filter(glance, r.squared <= 0.25)
gapminder %>%
    filter(country %in% unique(bad_fit$country)) %>%
    ggplot(aes(year, lifeExp, color = country, group = country)) +
    geom_line()

Explanation is:

• Rwanda genocide in late 1980’s through early 2000’s.
• AIDS epidemic affecting African countries in 1990s through 2000s.

Exercises

I focus specifically on Exercise 1 as it pertains to extending the Gapminder analysis.

Exercise 1

A linear trend seems to be slightly too simple for the overall trend. Can you do better with a quadratic polynomial? How can you interpret the coefficients of the quadratic? (Hint you might want to transform year so that it has mean zero.)

Instead of jumping to 4th order polynomial, I first want to see which model appears to have the best fit on a subset of the data. First, I create a polynomial model with splines::ns() on a subset of the data. I use the New Zealand subset, nz. I loop through models with increasing polynomial complexity.

library(splines)
# Plot polynomial functions from order 1 to 5 on nz data
for (i in 1:5) {
    # Model
    nz_mod <- lm(lifeExp ~ ns(year, i), data = nz)
    
    # Visualize data
    p1 <- nz %>%
        ggplot(aes(year, lifeExp)) + 
        geom_line() + 
        ggtitle("Full data = ")
    
    # Visualize model
    p2 <- nz %>%
        add_predictions(nz_mod) %>%
        ggplot(aes(year, pred)) +
        geom_line() + 
        ggtitle("Linear trend + ")
    
    # Visualize residuals
    p3 <- nz %>%
        add_residuals(nz_mod) %>%
        ggplot(aes(year, resid)) +
        geom_line() + 
        ggtitle("Remaining pattern") +
        # Set limits to see improvement on same scale
        scale_y_continuous(limits = c(-1.2, 1.2)) 
    
    # Output in grid    
    grid.arrange(p1, p2, p3, ncol=3,
                 top = stringr::str_c("Polynomial: Order = ", i, sep = ""))
}

Fourth and fifth order polynomials have much lower residuals in the remaining pattern. Let’s go with the quadratic.

Next, we apply the model to the full data set as part of a nested data frame.

by_country <- gapminder %>%
    group_by(country, continent) %>%
    nest()
country_model_4 <- function(df) {
    lm(lifeExp ~ ns(year, 4), data = df)
}
by_country <- by_country %>%
    mutate(model = map(.x = data, .f = country_model_4))
by_country %>% head()

Assess the model quality:

glance <- by_country %>%
    mutate(glance = map(model, broom::glance)) %>%
    unnest(glance, .drop = TRUE) 
glance %>%
    ggplot(aes(continent, r.squared)) +
    geom_jitter() +
    scale_y_continuous(limits = c(0,1))

The r.squared improved considerably with a more flexible, quadratic model.

Let’s visualize outliers with a boxplot:

glance %>%
    ggplot(aes(continent, r.squared)) +
    geom_boxplot()

Use r.squared <= 0.8 as a limit to visualize “low” model fits.

bad_fit <- filter(glance, r.squared <= 0.8)
preds <- by_country %>%
    filter(country %in% unique(bad_fit$country)) %>%
    mutate(preds = map2(.x = data, .y = model, .f = add_predictions)) %>%
    unnest(preds, .drop = TRUE) %>%
    select(-(pop:gdpPercap)) 
preds %>%
    ggplot(aes(year, lifeExp, color = country, group = country)) +
    geom_line() + 
    geom_line(aes(x = year, y = pred, color = country, group = country), 
              linetype = "dashed", data = preds)

The deviations from the model are due to:

1. Cambodia: Impact of the Vietnam war in through the 1970s.

2. Rwanda: Impact of Rwanda genocide (as noted previously).

In summary, the quadratic model does a better job fitting the minor fluctuations in the data. Thus, outliers from the data significantly deviate from the normal pattern. The largest outliers tend to have dramatic effects due to significant events in the country’s history that caused life expectency to deviate from the norm.

---

Details behind Creating Many Models

This section is meant to go into greater detail on the functions used in the Gapminder analysis. It specifically focuses on creating list-columns within nested data frames and then simplifying list-columns into unnested (flat) data frames.

Creating list-columns (4 methods)

1. With nesting: group_by() + tidyr::nest(): to convert grouped data that has list-column

2. From vectorized functions: mutate() + vectorized function that returns a list

3. From multivalued summaries: summarize() + summary function that returns a list

4. From a named list: tibble::enframe()

With nesting

gapminder %>%
    group_by(country, continent) %>%
    nest() %>%
    head()

Alternatively, can use on ungrouped data frame by specifying data columns:

gapminder %>%
    nest(year:gdpPercap) %>%
    head()

From vectorized functions

df <- tibble(x1 = c("a,b,c", "d,e,f,g"))
df %>%
    mutate(x2 = stringr::str_split(x1, ","))

Use unnest() to handle this.

df %>%
    mutate(x2 = stringr::str_split(x1, ",")) %>%
    unnest()

Can also use tidyr::separate_rows():

df %>% separate_rows(x1, sep = ",")

Using map(), map2(), pmap(), or invoke_map():

sim <- tribble(
    ~f,      ~params,
    "runif", list(min = -1, max = -1),
    "rnorm", list(sd = 5),
    "rpois", list(lambda = 10)
)
sim %>%
    mutate(sims = invoke_map(f, params, n = 10))

From multivalued summaries

Problem: One restriction of summarise() is that is only works with summary functions that return a single value. For example, cannot be used with quantile().

mtcars %>%
    group_by(cyl) %>%
    summarise(q = quantile(mpg))

Above code results in error. Pass functions that return lists as a list to get nested data frame. Can then unnest.

probs <- c(0.01, 0.25, 0.5, 0.75, 0.99)
mtcars %>%
    group_by(cyl) %>%
    summarise(p = list(probs), 
              q = list(quantile(mpg))) %>%
    unnest()

From a named list

tibble::enframe() to create nested data frame from nested list.

x <- list(
    a = 1:5,
    b = 3:4, 
    c = 5:6
    )
x

$a
[1] 1 2 3 4 5

$b
[1] 3 4

$c
[1] 5 6

enframe(x) %>% unnest()

Can iterate over names and values using map2().

df <- enframe(x)
df

df %>%
    mutate(summary = map2_chr(name, value, ~ stringr::str_c(.x, ": ", .y[[1]])))

Simplifying list-columns

Need to be able to convert list columns back to atomic vector or set of atomic vectors (i.e. a flat data frame). Two methods:

1. If want single value: mutate() + (map_lgl() | map_int() | map_dbl() | map_chr()). These map functions create atomic vectors.

2. If want many values: unnest()

List to vector

df <- tribble(
  ~x,
  letters[1:5],
  1:3,
  runif(5)
)
  
df %>% mutate(
  type = map_chr(x, typeof),
  length = map_int(x, length)
)

Shortcuts: Use names of nested list items to drill into lists.

df <- tribble(
  ~x,
  list(a = 1, b = 2),
  list(a = 2, c = 4)
)
df %>% mutate(
  a = map_dbl(x, "a", .null = NA_real_),
  b = map_dbl(x, "b", .null = NA_real_) # Throws error if no NA handling
)

Unnesting

Unnest to get the entire sub-list items.

tibble(x = 1:2, y = list(1:4, 1))

tibble(x = 1:2, y = list(1:4, 1)) %>% unnest()