Using gapminder to show how tidyr and broom and dplyr can fit a lot of models all at once.

nesting_models.R

library(tidyverse)
library(gapminder)
library(broom)
library(ggplot2)

life_exp_mods <- gapminder %>%
  
  #do this for each country
  group_by(country) %>%
  
  #nest the data frame
  nest(.key="country_data") %>%
  
  #now fit models and turn them into a list column
  mutate(mods = map(country_data, ~lm(lifeExp ~ year, data = .)))




#take a look at the object
life_exp_mods


#cool, right? There's a lot we can do with this. Get coefficients
life_exp_coefs <- life_exp_mods %>%
  
  #now get coefficients
  mutate(coefs = map(mods, tidy)) %>%
  
  #now pull it back out
  unnest(coefs) %>%
  
  #just the slopes
  filter(term=="year")

#who are the most "significant"?
arrange(life_exp_coefs, p.value)

#who are the most "biggest"?
arrange(life_exp_coefs, desc(estimate))


#Let's plot things! Volcano style!

ggplot(life_exp_coefs, 
       aes(x=estimate, y = log(p.value), label=country)) +
  geom_text() +
  geom_hline(yintercept=log(0.0005), lty=2, color="red")


### or plot fits!
life_exp_fits <- life_exp_mods %>%
  
  #get fitted values
  mutate(fits = map(mods, fitted)) %>%
  
  #extract fitted values and original data
  unnest(fits, country_data)

ggplot(life_exp_fits,
       aes(x=year, y=fits, group=country, color=country)) +
  geom_line() +
  geom_point(alpha=0.3, mapping=aes(y=lifeExp)) +
  facet_wrap(~continent) +
  guides(color = "none") +
  theme_bw()