I often need to run multiple sets of regressions on the same or similar datasets. This is usually for some set of robustness checks, either to help me better understand the stability of some result or to respond to referee requests. For years, this has been a mostly ad-hoc process for me. First, I would just copy-paste my regressions, modifying one variable or filter on the dataset with each paste. When this got to be too manual, I turned to nested loops and/or apply functions. This was an improvement in terms of running the regressions in a more systematic way, but extracting results I wanted to look at or highlight easily wasn’t straightforward. However, the purrr package (part of the tidyverse) provides tools that can make all of this easier.
The following code, along with a couple functions I’ve added to baylisR (call it a vanity package), allows me to facilitate a few common tasks:
I can easily build a set of regressions I want to run by combining different possible variables and datasets.
The output can be saved compactly as a tibble with a list-column containing either a stripped-down version of the fitted felm object, or a tidied data.frame of the same.
I can easily select specific statistical specifications for display in tables or figures.
There’s a bit more beneath the hood of this code. You’ll want to refer to the code of reg_felm to see how to call it. If you’re going to implement this yourself, I recommend you don’t rely on baylisR but instead extract the code for reg_felm (and strip_felm, which it calls) and modify it for your own purposes.
pacman::p_load(tidyverse, huxtable, stargazer, plm)# Load datadata("Wages", package ="plm")# Convenience function, fits model use lfe::felm following inputsreg_felm <-function(lhs, rhs, fe ="0", iv ="0", clus ="0", data_str) { data <-eval(parse(text = data_str)) fmla <-sprintf("%s ~ %s | %s | %s | %s", lhs, rhs, fe, iv, clus) fit <- lfe::felm(as.formula(fmla), data = data)}reg_tibble <-as_tibble(expand.grid(lhs ="lwage", rhs =c("exp", "exp + wks"), fe ="married", clus ="0", data_str =c("Wages", "Wages %>% filter(bluecol == 'no')"),stringsAsFactors = F))reg_tibble$fit <- purrr::pmap(reg_tibble, reg_felm)# These can be used directly in stargazer or huxtable...stargazer::stargazer(reg_tibble$fit)