Consider the following simplistic example of fixed effects regressions:
# Load packages
packs <- list("texreg", "plm", "lmtest")
lapply(packs, require, character.only = T)
# Load mtcars data set
d <- mtcars
# Run fixed effects regressions
fe_results <- lapply(c("hp", "drat", "qsec", "vs"), function(x) {
plm_eq <- paste0("mpg ~ wt + ", x)
plm_outp <- plm(plm_eq, index = "cyl", data = d, model = "within")
return(plm_outp)
})
# Print via texreg
texreg(fe_results, stars = c(0.01, 0.05, 0.1), include.rsquared = F, include.rmse = F, include.adjrsq = T)
The last line produces this table:
\begin{table}
\begin{center}
\begin{tabular}{l c c c c }
\hline
& Model 1 & Model 2 & Model 3 & Model 4 \\
\hline
wt & $-3.18^{***}$ & $-3.24^{***}$ & $-3.89^{***}$ & $-3.30^{***}$ \\
& $(0.72)$ & $(0.83)$ & $(0.91)$ & $(0.78)$ \\
hp & $-0.02^{*}$ & & & \\
& $(0.01)$ & & & \\
drat & & $-0.14$ & & \\
& & $(1.32)$ & & \\
qsec & & & $0.50$ & \\
& & & $(0.38)$ & \\
vs & & & & $0.86$ \\
& & & & $(1.64)$ \\
\hline
Adj. R$^2$ & 0.39 & 0.30 & 0.34 & 0.31 \\
Num. obs. & 32 & 32 & 32 & 32 \\
\hline
\multicolumn{5}{l}{\scriptsize{$^{***}p<0.01$, $^{**}p<0.05$, $^*p<0.1$}}
\end{tabular}
\caption{Statistical models}
\label{table:coefficients}
\end{center}
\end{table}
However, I would like to change the standard errors into heteroskedasticity-robust standard errors. Therefore, I add coeftest to the given function:
fe_results <- lapply(c("hp", "drat", "qsec", "vs"), function(x) {
plm_eq <- paste0("mpg ~ wt + ", x)
plm_outp <- plm(plm_eq, index = "cyl", data = d, model = "within")
plm_outp <- coeftest(plm_outp, vcov = vcovHC(plm_outp, type = "HC1"))
return(plm_outp)
})
Unfortunately, coeftest has the undesired effect of dropping information on the Adj. R$^2$ and the Num. obs., which then also removes it from the texreg output:
\begin{table}
\begin{center}
\begin{tabular}{l c c c c }
\hline
& Model 1 & Model 2 & Model 3 & Model 4 \\
\hline
wt & $-3.18^{***}$ & $-3.24^{***}$ & $-3.89^{***}$ & $-3.30^{***}$ \\
& $(0.95)$ & $(0.83)$ & $(1.29)$ & $(1.08)$ \\
hp & $-0.02^{*}$ & & & \\
& $(0.01)$ & & & \\
drat & & $-0.14$ & & \\
& & $(0.72)$ & & \\
qsec & & & $0.50^{***}$ & \\
& & & $(0.10)$ & \\
vs & & & & $0.86$ \\
& & & & $(0.63)$ \\
\hline
\multicolumn{5}{l}{\scriptsize{$^{***}p<0.01$, $^{**}p<0.05$, $^*p<0.1$}}
\end{tabular}
\caption{Statistical models}
\label{table:coefficients}
\end{center}
\end{table}
One solution that found circumvents dropping Adj. R$^2$ and Num. obs. via:
fe_results <- lapply(c("hp", "drat", "qsec", "vs"), function(x) {
plm_eq <- paste0("mpg ~ wt + ", x)
plm_outp <- plm(plm_eq, index = "cyl", data = d, model = "within")
plm_outp_s <- summary(plm_outp)
plm_outp_s$coefficients <- unclass(coeftest(plm_outp, vcov = vcovHC(plm_outp, type = "HC1")))
return(plm_outp_s)
})
The disadvantage is that the returned output is now in summary.plm format and cannot be printed by texreg anymore (Error in (function (classes, fdef, mtable): unable to find an inherited method for function ‘extract’ for signature ‘"summary.plm"’). How do I fix this? How do I print regression output incl. Adj. R$^2$, Num. obs. and heteroskedasticity robust standard errors using texreg?
Furthermore, I would like to add an "intercept" in a Stata-like manner. I figured that the solution might rely on within_intercept() but am uncertain of how to use it in a way that it gets included into the texreg output.
I tried to keep this example as simple and coherent as possible and am looking forward to any comments and suggestions.
