Search code examples
rmodelglmmodel-fitting

How to fit multiple interaction models in a loop?


Lets say i have 3 response variables A,C and M and i want to fit a model for all possible models ie fit Y ~ A, Y ~ C, Y ~ M, Y ~ A * C, Y ~ A * M, Y ~ C * M, etc. Is there a quick way to do this without manually specifiying the interactions each time?

i do not want to write

M1 = glm(Y ~ A , data = subs, family = "poisson")
M2 = glm(Y ~ C , data = subs, family = "poisson")
M3 = glm(Y ~ M , data = subs, family = "poisson")
M4 = glm(Y ~ A*C , data = subs, family = "poisson")
...

In reality i have more than 3 variables and would like some sort of loop, is this even possible. Thanks


Solution

  • Here is a sort of functional programming approach. You create your data, and as long as your Y is the first column, this code would take all the rest of the variables (no matter how many) and construct models on their combinations.

    Finally, since you've done it in this framework, you can call broom's tidy and confint_tidy to extract the results into an easy to filter dataset.

    DF <- data_frame(Y = rpois(100, 5),
               A = rnorm(100),
               C = rnorm(100),
               M = rnorm(100))
    
    formula_frame <- bind_rows(data_frame(V1 = names(DF[,-1])),
                               as_data_frame(t(combn(names(DF[,-1]),2)))) %>%
      rowwise() %>%
      mutate(formula_text = paste0("Y ~", if_else(is.na(V2),
                                                  V1, 
                                                  paste(V1,V2, sep = "*"))), 
             formula_obj = list(as.formula(formula_text))) %>%
      ungroup()
    
    formula_frame %>% 
    mutate(fits = map(formula_obj, ~glm(.x, family = "poisson", data = DF) %>%
    (function(X)bind_cols(broom::tidy(X),broom::confint_tidy((X)))))) %>%
     unnest(fits) %>%
     select(-formula_obj)
    
    # A tibble: 18 x 10
       V1    V2    formula_text term        estimate std.error statistic   p.value conf.low conf.high
       <chr> <chr> <chr>        <chr>          <dbl>     <dbl>     <dbl>     <dbl>    <dbl>     <dbl>
     1 A     NA    Y ~A         (Intercept)  1.63       0.0443    36.8   6.92e-297   1.54      1.72  
     2 A     NA    Y ~A         A            0.0268     0.0444     0.602 5.47e-  1  -0.0603    0.114 
     3 C     NA    Y ~C         (Intercept)  1.63       0.0443    36.8   5.52e-296   1.54      1.72  
     4 C     NA    Y ~C         C            0.0326     0.0466     0.699 4.84e-  1  -0.0587    0.124 
     5 M     NA    Y ~M         (Intercept)  1.63       0.0454    35.8   1.21e-280   1.53      1.71  
     6 M     NA    Y ~M         M           -0.0291     0.0460    -0.634 5.26e-  1  -0.119     0.0615
     7 A     C     Y ~A*C       (Intercept)  1.62       0.0446    36.4   5.64e-290   1.54      1.71  
     8 A     C     Y ~A*C       A            0.00814    0.0459     0.178 8.59e-  1  -0.0816    0.0982
     9 A     C     Y ~A*C       C            0.0410     0.0482     0.850 3.96e-  1  -0.0532    0.136 
    10 A     C     Y ~A*C       A:C          0.0650     0.0474     1.37  1.70e-  1  -0.0270    0.158 
    11 A     M     Y ~A*M       (Intercept)  1.62       0.0458    35.5   1.21e-275   1.53      1.71  
    12 A     M     Y ~A*M       A            0.0232     0.0451     0.514 6.07e-  1  -0.0653    0.112 
    13 A     M     Y ~A*M       M           -0.0260     0.0464    -0.561 5.75e-  1  -0.116     0.0655
    14 A     M     Y ~A*M       A:M         -0.00498    0.0480    -0.104 9.17e-  1  -0.0992    0.0887
    15 C     M     Y ~C*M       (Intercept)  1.60       0.0472    34.0   1.09e-253   1.51      1.70  
    16 C     M     Y ~C*M       C            0.0702     0.0506     1.39  1.65e-  1  -0.0291    0.169 
    17 C     M     Y ~C*M       M           -0.0333     0.0479    -0.695 4.87e-  1  -0.127     0.0611
    18 C     M     Y ~C*M       C:M          0.0652     0.0377     1.73  8.39e-  2  -0.0102    0.138