Use purrr's map to estimate linear regressions by group when some groups are empty for some models

user_334755 0
Viewed 142

I would like to estimate many linear regressions for many groups. I use the tidyverse, so I tried purrr's map and broom's tidy.

However, not all groups have observations for all variables in all models. In my example below, z is missing when t is 1, so lm cannot estimate y ~ x + z when t is 1. I thought a filter statement in each lm would solve the problem. However, sometimes filter provides an empty data set, and lm throws an error.

I thought that map would have an option to handle these cases, but I do not see one in the help file.

Is there a best practice here? FWIW, I only want the coefficient estimates. If you swap the commented mutate the code works as expected.

library(tidyverse)


df <- tibble(t = rep(1:2, each = 10),
             y = runif(20),
             x = runif(20)) %>%
    mutate(z = ifelse(t == 1, NA, runif(20)))
    # mutate(z = runif(20))


results <- df %>%
    nest(dat = -t) %>%
    mutate(
        model_1 = map(dat, ~ lm(y ~ x, data = .x %>% drop_na(y, x))),
        model_2 = map(dat, ~ lm(y ~ x + z, data = .x %>% drop_na(x, z))),
        coef_1 = map(model_1, tidy),
        coef_2 = map(model_2, tidy)
    ) %>%
    select(t, starts_with('coef')) %>% 
    pivot_longer(
        cols = starts_with('coef')
    ) %>% 
    unnest(value)
#> Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...): 0 (non-NA) cases

Created on 2020-05-15 by the reprex package (v0.3.0)

purrr r
1 Answers

In unnest, the default is keep_empty = FALSE, but here the issue is at the model creation as all the values in 'z' are NA. We can create a condition with if/else

library(dplyr)
library(purrr)
library(tidyr)
df %>%
     nest(dat = -t) %>%
     mutate(
         model_1 = map(dat, ~ lm(y ~ x, data = .x %>% drop_na(y, x))),
         model_2 = map(dat, ~ if(all(is.na(.x$z))) NULL else 
               lm(y ~ x + z, data = .x %>% drop_na(x, z))),
         coef_1 = map(model_1, tidy),
         coef_2 = map(model_2, tidy)) %>% 
      select(t, starts_with('coef')) %>% 
      pivot_longer(
         cols = starts_with('coef')
     ) %>% 
     unnest(c(value))
# A tibble: 7 x 7
#      t name   term        estimate std.error statistic p.value
#  <int> <chr>  <chr>          <dbl>     <dbl>     <dbl>   <dbl>
#1     1 coef_1 (Intercept)    0.333     0.196     1.70  0.127  
#2     1 coef_1 x              0.444     0.386     1.15  0.283  
#3     2 coef_1 (Intercept)    0.852     0.195     4.37  0.00239
#4     2 coef_1 x             -0.854     0.363    -2.35  0.0463 
#5     2 coef_2 (Intercept)    0.597     0.438     1.36  0.215  
#6     2 coef_2 x             -0.659     0.479    -1.37  0.212  
#7     2 coef_2 z              0.280     0.425     0.658 0.532
user_37322710
Related