```
::p_load(tidyverse, huxtable, stargazer, plm)
pacman
# Load data
data("Wages", package = "plm")
# Convenience function, fits model use lfe::felm following inputs
<- function(lhs, rhs, fe = "0", iv = "0", clus = "0", data_str) {
reg_felm <- eval(parse(text = data_str))
data <- sprintf("%s ~ %s | %s | %s | %s", lhs, rhs, fe, iv, clus)
fmla <- lfe::felm(as.formula(fmla), data = data)
fit
}
<- as_tibble(
reg_tibble expand.grid(lhs = "lwage", rhs = c("exp", "exp + wks"),
fe = "married", clus = "0", data_str = c("Wages", "Wages %>% filter(bluecol == 'no')"),
stringsAsFactors = F))
$fit <- purrr::pmap(reg_tibble, reg_felm)
reg_tibble
# These can be used directly in stargazer or huxtable...
::stargazer(reg_tibble$fit) stargazer
```

```
% 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}
```

`::huxreg(reg_tibble$fit) huxtable`

(1) | (2) | (3) | (4) | |
---|---|---|---|---|

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) | |||

N | 4165 | 4165 | 2036 | 2036 |

R2 | 0.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_id | term | estimate | std.error | statistic | p.value |
---|---|---|---|---|---|

1 | exp | 0.00705 | 0.000624 | 11.3 | 3.32e-29 |

2 | exp | 0.00714 | 0.000623 | 11.5 | 6.39e-30 |

2 | wks | 0.00433 | 0.00132 | 3.29 | 0.00102 |

3 | exp | 0.0116 | 0.000889 | 13.1 | 1.26e-37 |

4 | exp | 0.0117 | 0.000886 | 13.2 | 3.71e-38 |

4 | wks | 0.0079 | 0.00195 | 4.05 | 5.26e-05 |

```
# It's also straightforward to display only a subset of the regressions.
::huxreg(reg_tibble %>% filter(rhs == "exp + wks") %>% pull(fit)) huxtable
```

(1) | (2) | |
---|---|---|

exp | 0.007 *** | 0.012 *** |

(0.001) | (0.001) | |

wks | 0.004 ** | 0.008 *** |

(0.001) | (0.002) | |

N | 4165 | 2036 |

R2 | 0.112 | 0.166 |

*** p < 0.001; ** p < 0.01; * p < 0.05. |