For my stargazer table I would like to include the diagnostics of weak instrument and Wu Hausman test.
I read the following answer: R: Robust SE's and model diagnostics in stargazer table.
Though, I don't really understand this. How can you apply it to a larger model with 9 or 10 columns?
How can I insert the diagnostics you get when you code summary(ivreg1, diagnostics=TRUE)?
I have the following code for the stargazer table:
ivreg1 <- ivreg(budget ~ policyfactor1 + control1 + control2 + control3 | iv + control1 + control2 + control3, data=df)
ivreg2 <- ivreg(budget ~ policyfactor2 + control1 + control2 + control3 | iv + control1 + control2 + control3, data=df)
ivreg3 <- ivreg(budget ~ policyfactor3 + control1 + control2 + control3 | iv + control1 + control2 + control3, data=df)
ivreg4 <- ivreg(budget ~ policyfactor4 + control1 + control2 + control3 | iv + control1 + control2 + control3, data=df)
ivreg5 <- ivreg(budget ~ policyfactor5 + control1 + control2 + control3 | iv + control1 + control2 + control3, data=df)
ivreg6 <- ivreg(budget ~ policyfactor6 + control1 + control2 + control3 | iv + control1 + control2 + control3, data=df)
ivreg7 <- ivreg(budget ~ policyfactor7 + control1 + control2 + control3 | iv + control1 + control2 + control3, data=df)
ivreg8 <- ivreg(budget ~ policyfactor8 + control1 + control2 + control3 | iv + control1 + control2 + control3, data=df)
ivreg9 <- ivreg(budget ~ policyfactor9 + control1 + control2 + control3 | iv + control1 + control2 + control3, data=df)
stargazer(ivreg1, ivreg2, ivreg3, ivreg4, ivreg5, ivreg6, ivreg7, ivreg8, ivreg9, title="Regression results", align=TRUE, column.sep.width = "-15pt", font.size = "small", type="latex")
I have the following dataframe example:
data.frame(
budget = c(10, 8, -7, 8, 3, 2, 0.5, 1.5, -2, -15, 30, -0.5, 12, 18),
policyfactor1 = c(1L, 1L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 1L, 0L, 1L, 1L),
policyfactor2 = c(0L, 1L, 0L, 1L, 0L, 0L, 1L, 1L, 0L, 0L, 1L, 0L, 1L, 0L),
policyfactor3 = c(0L, 1L, 1L, 0L, 1L, 1L, 1L, 0L, 0L, 1L, 1L, 1L, 0L, 0L),
policyfactor4 = c(1L, 0L, 1L, 0L, 0L, 1L, NA, NA, 1L, 1L, 0L, 1L, 1L, 0L),
policyfactor5 = c(1L, 0L, 1L, 0L, 0L, 1L, 1L, 0L, 0L, 1L, 1L, 0L, 1L, 1L),
policyfactor6 = c(0L, 1L, 1L, 0L, 1L, 1L, 0L, 1L, 1L, 1L, 0L, 1L, NA, 0L),
policyfactor7 = c(0L, 1L, 1L, 1L, 1L, 0L, 1L, 0L, 1L, 0L, 0L, 1L, 0L, 0L),
policyfactor8 = c(1L, 0L, 1L, 1L, NA, 1L, 1L, 0L, 0L, 1L, 0L, 0L, 1L, 0L),
policyfactor9 = c(0L, 0L, 1L, 0L, 1L, 0L, 1L, 1L, 0L, 1L, NA, 1L, 1L, 1L),
IV = c(17L,46L,18L,23L,35L,10L,0L,19L,
15L,12L,21L,6L,27L,10L),
control1 = c(29.4,27.4,33.2,27.1,27.4,34.2,
26.8,32.9,26,26.3,27.3,33.4,23.5,31.3),
control2 = c(1.3,20,-0.2,3,3.4,0.3,-1.1,1.9,
-2,-1.6,2.8,1.9,2,1.8)
)
#> budget policyfactor1 policyfactor2 policyfactor3 policyfactor4 policyfactor5
#> 1 10.0 1 0 0 1 1
#> 2 8.0 1 1 1 0 0
#> 3 -7.0 0 0 1 1 1
#> 4 8.0 0 1 0 0 0
#> 5 3.0 0 0 1 0 0
#> 6 2.0 1 0 1 1 1
#> 7 0.5 1 1 1 NA 1
#> 8 1.5 0 1 0 NA 0
#> 9 -2.0 0 0 0 1 0
#> 10 -15.0 0 0 1 1 1
#> 11 30.0 1 1 1 0 1
#> 12 -0.5 0 0 1 1 0
#> 13 12.0 1 1 0 1 1
#> 14 18.0 1 0 0 0 1
#> policyfactor6 policyfactor7 policyfactor8 policyfactor9 IV control1 control2
#> 1 0 0 1 0 17 29.4 1.3
#> 2 1 1 0 0 46 27.4 20.0
#> 3 1 1 1 1 18 33.2 -0.2
#> 4 0 1 1 0 23 27.1 3.0
#> 5 1 1 NA 1 35 27.4 3.4
#> 6 1 0 1 0 10 34.2 0.3
#> 7 0 1 1 1 0 26.8 -1.1
#> 8 1 0 0 1 19 32.9 1.9
#> 9 1 1 0 0 15 26.0 -2.0
#> 10 1 0 1 1 12 26.3 -1.6
#> 11 0 0 0 NA 21 27.3 2.8
#> 12 1 1 0 1 6 33.4 1.9
#> 13 NA 0 1 1 27 23.5 2.0
#> 14 0 0 0 1 10 31.3 1.8
Created on 2020-06-20 by the reprex package (v0.3.0)
You just need to list the model fits where you may use mget if you have the objects already in your global environment (see alternative, though, at bottom of the answer).
iv.fit <- mget(paste0("ivreg", 1:9))
Then in stargazer::stargazer, in the gaze.lines.ivreg.diagn function add the summaries with included diagnostics using lapply.
library(stargazer)
stargazer(iv.fit[1:3], title="Regression results", align=TRUE,
type="text", ## to get LaTeX output use `type="latex"`
column.sep.width="-15pt", font.size="small",
add.lines=
gaze.lines.ivreg.diagn(lapply(iv.fit[1:3], summary, diagnostics=TRUE), row=1:2))
stargazer# Regression results
# ===========================================================
# Dependent variable:
# -----------------------------
# budget
# (1) (2) (3)
# -----------------------------------------------------------
# policyfactor1 1.217
# (15.134)
#
# policyfactor2 3.209
# (39.816)
#
# policyfactor3 3.302
# (44.912)
#
# control1 -0.394 -0.237 -0.513
# (0.984) (2.407) (1.687)
#
# control2 0.517 0.439 0.486
# (0.721) (1.496) (1.079)
#
# Constant 14.493 9.342 16.727
# (31.255) (82.795) (33.982)
#
# -----------------------------------------------------------
# Weak instruments 0.17 0.57 0.63
# Wu-Hausman 0.36 0.91 0.83
# Observations 14 14 14
# R2 0.154 0.159 -0.011
# Adjusted R2 -0.099 -0.094 -0.315
# Residual Std. Error (df = 10) 11.465 11.437 12.539
# ===========================================================
# Note: *p<0.1; **p<0.05; ***p<0.01
This is also possible with texreg::texreg. Similarly we put the summaries into a list.
iv.gof2 <- lapply(iv.fit, function(x)
summary(x, diagnostics=T)$diagnostics[c("Weak instruments", "Wu-Hausman"), "p-value"])
Then we use texreg::extract to extract relevant model details to edit the GOF information. This gives list iv.fit.ex containing "texreg" objects.
library(texreg)
iv.fit.ex <- lapply(iv.fit, extract)
With a small function we may add additional GOFs using Map.
addGof<- function(x, y) {
x@gof.names <- c(x@gof.names, names(y))
x@gof <- c(x@gof, y)
x@gof.decimal <- c(x@gof.decimal, rep(TRUE, length(y)))
x
}
iv.fit.ex <- Map(addGof, iv.fit.ex, iv.gof2)
Then we run texreg::texreg on iv.fit.ex. (I show console version texreg::screenreg here).
screenreg(iv.fit.ex[1:3], digits=3, column.spacing=1,
custom.model.names=sprintf(" (%s)", 1:length(iv.fit.ex[1:3])))
texreg# ===========================================
# (1) (2) (3)
# -------------------------------------------
# (Intercept) 14.493 9.342 16.727
# (31.255) (82.795) (33.982)
# policyfactor1 1.217
# (15.134)
# control1 -0.394 -0.237 -0.513
# (0.984) (2.407) (1.687)
# control2 0.517 0.439 0.486
# (0.721) (1.496) (1.079)
# policyfactor2 3.209
# (39.816)
# policyfactor3 3.302
# (44.912)
# -------------------------------------------
# R^2 0.154 0.159 -0.011
# Adj. R^2 -0.099 -0.094 -0.315
# Num. obs. 14 14 14
# Weak instruments 0.168 0.567 0.628
# Wu-Hausman 0.357 0.911 0.825
# ===========================================
# *** p < 0.001; ** p < 0.01; * p < 0.05
By the way, you can also list your regression models in a less tedious way, by not having to code each regression individually:
fo <- sapply(paste0("policyfactor", 1:9), function(i)
as.formula(sprintf("budget ~ %s + control1 + control2 | iv + control1 + control2", i)))
iv.fit <- lapply(fo, function(fo) do.call(ivreg, list(formula=fo, data=quote(dat))))
Data:
dat <- structure(list(budget = c(10, 8, -7, 8, 3, 2, 0.5, 1.5, -2, -15,
30, -0.5, 12, 18), policyfactor1 = c(1L, 1L, 0L, 0L, 0L, 1L,
1L, 0L, 0L, 0L, 1L, 0L, 1L, 1L), policyfactor2 = c(0L, 1L, 0L,
1L, 0L, 0L, 1L, 1L, 0L, 0L, 1L, 0L, 1L, 0L), policyfactor3 = c(0L,
1L, 1L, 0L, 1L, 1L, 1L, 0L, 0L, 1L, 1L, 1L, 0L, 0L), policyfactor4 = c(1L,
0L, 1L, 0L, 0L, 1L, NA, NA, 1L, 1L, 0L, 1L, 1L, 0L), policyfactor5 = c(1L,
0L, 1L, 0L, 0L, 1L, 1L, 0L, 0L, 1L, 1L, 0L, 1L, 1L), policyfactor6 = c(0L,
1L, 1L, 0L, 1L, 1L, 0L, 1L, 1L, 1L, 0L, 1L, NA, 0L), policyfactor7 = c(0L,
1L, 1L, 1L, 1L, 0L, 1L, 0L, 1L, 0L, 0L, 1L, 0L, 0L), policyfactor8 = c(1L,
0L, 1L, 1L, NA, 1L, 1L, 0L, 0L, 1L, 0L, 0L, 1L, 0L), policyfactor9 = c(0L,
0L, 1L, 0L, 1L, 0L, 1L, 1L, 0L, 1L, NA, 1L, 1L, 1L), iv = c(17L,
46L, 18L, 23L, 35L, 10L, 0L, 19L, 15L, 12L, 21L, 6L, 27L, 10L
), control1 = c(29.4, 27.4, 33.2, 27.1, 27.4, 34.2, 26.8, 32.9,
26, 26.3, 27.3, 33.4, 23.5, 31.3), control2 = c(1.3, 20, -0.2,
3, 3.4, 0.3, -1.1, 1.9, -2, -1.6, 2.8, 1.9, 2, 1.8)), class = "data.frame", row.names = c(NA,
-14L))