Applying multiple model formulas to groups of data

user_697065 0
Viewed 592

I'd like to apply 3 linear models to my data, and extract the residuals for each. I am wondering if there is a way to apply the same steps for each model using a combination of dplyr and purrr:

I want to keep:

  1. The lm object for each model
  2. The augment output for each model
  3. The residuals for each model

Here's a working example that analyzes the mpg dataset:

library(dplyr)
library(tidyr)
library(purrr)
library(broom)
library(ggplot2)

Here are the three different formulas I want to use for my lm 

f1 = hwy ~ cyl
f2 = hwy ~ displ
f3 = hwy ~ cyl + displ

lin_mod = function(formula) {
  function(data) {
    lm(formula, data = data)
  }
}

This is how I extract residuals for a single formula:

mpg %>% 
group_by(manufacturer) %>% 
nest() %>% 
mutate(model = map(data, lin_mod(f1)), 
       aug = map(model, augment), 
       res = map(aug, ".resid"))

However, this technique seems like a bad way to do it for all the formulas, since I rewrite a lot of code:

mpg %>% 
group_by(manufacturer) %>% 
nest() %>% 
mutate(model1 = map(data, lin_mod(f1)), 
       aug1 = map(model1, augment), 
       res1 = map(aug1, ".resid"),
       model2 = map(data, lin_mod(f2)), 
       aug2 = map(model2, augment), 
       res2 = map(aug2, ".resid"),
       model3 = map(data, lin_mod(f3)), 
       aug3 = map(model3, augment), 
       res3 = map(aug3, ".resid"))

How do I apply this function to each formula in an elegant way? I was thinking that mutate_all, or putting the formulas into a list might help in some way, but alas I'm stuck.

purrr dplyr r
3 Answers

I was looking into how to do a version of this — apply a list of model formulas to the same data — and came up with a more succinct way to do it with tidyverse version 1.3.2.

At a high level, the idea is to "map formulas to data" instead of "map data to formulas".

library("tidyverse")

formulas <- list(
  f1 = hwy ~ cyl,
  f2 = hwy ~ displ,
  f3 = hwy ~ cyl + displ
)

Once we apply the formulas to the data, we enframe the list of fitted lm models in a tibble.

formulas %>%
  map(
    ~ lm(., mpg)
  ) %>%
  enframe()
#> # A tibble: 3 × 2
#>   name  value 
#>   <chr> <list>
#> 1 f1    <lm>  
#> 2 f2    <lm>  
#> 3 f3    <lm>

Let's write a function to apply the list of formulas to any data.

fit_model_to_data <- function(data, ...) {
  formulas %>%
    map(
      ~ lm(., data)
    ) %>%
    enframe(
      name = "formula",
      value = "model"
    )
}

Once we have the helper function, we can group by manufacturer and fit the formulas to each split.

fits <- mpg %>%
  group_by(
    manufacturer
  ) %>%
  group_modify(
    fit_model_to_data
  )
fits
#> # A tibble: 45 × 3
#> # Groups:   manufacturer [15]
#>    manufacturer formula model 
#>    <chr>        <chr>   <list>
#>  1 audi         f1      <lm>  
#>  2 audi         f2      <lm>  
#>  3 audi         f3      <lm>  
#>  4 chevrolet    f1      <lm>  
#>  5 chevrolet    f2      <lm>  
#>  6 chevrolet    f3      <lm>  
#>  7 dodge        f1      <lm>  
#>  8 dodge        f2      <lm>  
#>  9 dodge        f3      <lm>  
#> 10 ford         f1      <lm>  
#> # … with 35 more rows

At this point we have the fitted models in a handy data frame. It's easy to extract the residuals.

fits %>%
  mutate(
    .resid = map(model, residuals)
  )
#> # A tibble: 45 × 4
#> # Groups:   manufacturer [15]
#>    manufacturer formula model  .resid    
#>    <chr>        <chr>   <list> <list>    
#>  1 audi         f1      <lm>   <dbl [18]>
#>  2 audi         f2      <lm>   <dbl [18]>
#>  3 audi         f3      <lm>   <dbl [18]>
#>  4 chevrolet    f1      <lm>   <dbl [19]>
#>  5 chevrolet    f2      <lm>   <dbl [19]>
#>  6 chevrolet    f3      <lm>   <dbl [19]>
#>  7 dodge        f1      <lm>   <dbl [37]>
#>  8 dodge        f2      <lm>   <dbl [37]>
#>  9 dodge        f3      <lm>   <dbl [37]>
#> 10 ford         f1      <lm>   <dbl [25]>
#> # … with 35 more rows


packageVersion("tidyverse")
#> [1] '1.3.2'

Created on 2022-09-06 with reprex v2.0.2

user_95664840
Related