Programatic way to create, name, store, and summarize models

Viewed 42

Suppose you have a data frame like this:

dat1 <- data.frame(Loc = rep(c("NY","MA","FL","GA"), each = 1000),
                   Region = rep(c("a","b","c","d"),each = 1000),
                   ID = rep(c(1:10), each=200),
                   var1 = rnorm(1000),
                   var2=rnorm(1000),
                   var3=rnorm(1000))

You start off by creating some empty lists to store things you are interested in, something like this:

models <- list()
fitvals <- list()
summaries <- list()
plots <- ()

Now for each combination of var# and ID you want to create a separate model, each named according to the Loc, Reg, ID, and var# combination it belongs to (e.g., the model name for Var1 / ID == 1 who is from NY, which is in Region1 would be something like: NYa1var1mod). Each of these models will be stored in the models list. Similarly, you want to store the fit values, standard errors, and confidence intervals for each model in a data frame with a similar naming convention and store these data frames in the list fitvals. Further, you want to produce a standard ggplot of some kind upon each (say, just a line plot), which has a title that follows the naming convention, and y-axis name that corresponds with the variable being plotted.

Obviously there are tons of ways to go about doing this in a programmatic fashion, I am just curious about the different approaches to tackling this situation and what the most computationally efficient methods would be.

1 Answers

I like to keep everything in data frames, e.g. like this:

library(dplyr)

# keeping all models stored
dat1 %>% group_by(Loc, Region, ID) %>% do(mod = lm(var1 ~ var2 + var3, data = .))

# A tibble: 20 x 4
   Loc   Region    ID mod     
 * <fct> <fct>  <int> <list>  
 1 FL    c          1 <S3: lm>
 2 FL    c          2 <S3: lm>
 3 FL    c          3 <S3: lm>

Here, you can access each model by simply getting the element of the dataframe.

# keeping the coefficients etc stored
library(broom)
dat1 %>% group_by(Loc, Region, ID) %>% do(tidy(lm(var1 ~ var2 + var3, data = .)))

# A tibble: 60 x 8
# Groups:   Loc, Region, ID [20]
   Loc   Region    ID term        estimate std.error statistic p.value
   <fct> <fct>  <int> <chr>          <dbl>     <dbl>     <dbl>   <dbl>
 1 FL    c          1 (Intercept)  0.0284     0.0748    0.380    0.704
 2 FL    c          1 var2         0.0786     0.0745    1.06     0.293
 3 FL    c          1 var3        -0.0994     0.0719   -1.38     0.169
 4 FL    c          2 (Intercept) -0.00290    0.0693   -0.0418   0.967
 5 FL    c          2 var2         0.0609     0.0759    0.803    0.423
 6 FL    c          2 var3        -0.116      0.0715   -1.62     0.106
Related