Making regressions purrr

programming
Author

Patrick Baylis

Published

June 11, 2019

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:

1. I can easily build a set of regressions I want to run by combining different possible variables and datasets.
2. 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.
3. 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)

data("Wages", package = "plm")

# Convenience function, fits model use lfe::felm following inputs
reg_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)

% Table created by stargazer v.5.2.3 by Marek Hlavac, Social Policy Institute. E-mail: marek.hlavac at gmail.com
% Date and time: Sat, Dec 10, 2022 - 23:35:58
\begin{table}[!htbp] \centering
\caption{}
\label{}
\begin{tabular}{@{\extracolsep{5pt}}lcccc}
\\[-1.8ex]\hline
\hline \\[-1.8ex]
& \multicolumn{4}{c}{\textit{Dependent variable:}} \\
\cline{2-5}
\\[-1.8ex] & \multicolumn{4}{c}{fmla} \\
\\[-1.8ex] & (1) & (2) & (3) & (4)\\
\hline \\[-1.8ex]
exp & 0.007$^{***}$ & 0.007$^{***}$ & 0.012$^{***}$ & 0.012$^{***}$ \\
& (0.001) & (0.001) & (0.001) & (0.001) \\
& & & & \\
wks &  & 0.004$^{***}$ &  & 0.008$^{***}$ \\
&  & (0.001) &  & (0.002) \\
& & & & \\
\hline \\[-1.8ex]
Observations & 4,165 & 4,165 & 2,036 & 2,036 \\
R$^{2}$ & 0.110 & 0.112 & 0.160 & 0.166 \\
Adjusted R$^{2}$ & 0.110 & 0.112 & 0.159 & 0.165 \\
Residual Std. Error & 0.435 (df = 4162) & 0.435 (df = 4161) & 0.415 (df = 2033) & 0.414 (df = 2032) \\
\hline
\hline \\[-1.8ex]
\textit{Note:}  & \multicolumn{4}{r}{$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01} \\
\end{tabular}
\end{table} 
huxtable::huxreg(reg_tibble$fit) (1)(2)(3)(4) exp0.007 ***0.007 ***0.012 ***0.012 *** (0.001) (0.001) (0.001) (0.001) wks 0.004 ** 0.008 *** (0.001) (0.002) N4165 4165 2036 2036 R20.110 0.112 0.160 0.166 *** p < 0.001; ** p < 0.01; * p < 0.05. # ...Or shown all together in a single table, after tidying. bind_rows(lapply(reg_tibble$fit, broom::tidy), .id = "reg_id")
reg_idtermestimatestd.errorstatisticp.value
1exp0.007050.00062411.3 3.32e-29
2exp0.007140.00062311.5 6.39e-30
2wks0.004330.00132 3.290.00102
3exp0.0116 0.00088913.1 1.26e-37
4exp0.0117 0.00088613.2 3.71e-38
4wks0.0079 0.00195 4.055.26e-05
# It's also straightforward to display only a subset of the regressions.
huxtable::huxreg(reg_tibble %>% filter(rhs == "exp + wks") %>% pull(fit))
(1)(2)
exp0.007 ***0.012 ***
(0.001)   (0.001)
wks0.004 ** 0.008 ***
(0.001)   (0.002)
N4165        2036
R20.112    0.166
*** p < 0.001; ** p < 0.01; * p < 0.05.