R: Robust SE's and model diagnostics in stargazer table
Here's one way to do what you want:
require(lmtest)
rob.fit1 <- coeftest(fit1, function(x) vcovHC(x, type="HC0"))
rob.fit2 <- coeftest(fit2, function(x) vcovHC(x, type="HC0"))
summ.fit1 <- summary(fit1, vcov. = function(x) vcovHC(x, type="HC0"), diagnostics=T)
summ.fit2 <- summary(fit2, vcov. = function(x) vcovHC(x, type="HC0"), diagnostics=T)
stargazer(fit1, fit2, type = "text",
se = list(rob.fit1[,"Std. Error"], rob.fit2[,"Std. Error"]),
add.lines = list(c(rownames(summ.fit1$diagnostics)[1],
round(summ.fit1$diagnostics[1, "p-value"], 2),
round(summ.fit2$diagnostics[1, "p-value"], 2)),
c(rownames(summ.fit1$diagnostics)[2],
round(summ.fit1$diagnostics[2, "p-value"], 2),
round(summ.fit2$diagnostics[2, "p-value"], 2)) ))
Which will yield:
==========================================================
Dependent variable:
----------------------------
y
(1) (2)
----------------------------------------------------------
x -1.222 -0.912
(1.672) (1.002)
a -0.240 -0.208
(0.301) (0.243)
Constant 9.662 8.450**
(6.912) (4.222)
----------------------------------------------------------
Weak instruments 0.45 0.56
Wu-Hausman 0.11 0.18
Observations 100 100
R2 -4.414 -2.458
Adjusted R2 -4.526 -2.529
Residual Std. Error (df = 97) 22.075 17.641
==========================================================
Note: *p<0.1; **p<0.05; ***p<0.01
As you can see, this allows manually including the diagnostics in the respective models.
You could automate this approach by creating a function that takes in a list of models (e.g. list(summ.fit1, summ.fit2)) and outputs the objects required by se or add.lines arguments.
gaze.coeft <- function(x, col="Std. Error"){
stopifnot(is.list(x))
out <- lapply(x, function(y){
y[ , col]
})
return(out)
}
gaze.coeft(list(rob.fit1, rob.fit2))
gaze.coeft(list(rob.fit1, rob.fit2), col=2)
Will both take in a list of coeftest objects, and yield the SEs vector as expected by se:
[[1]]
(Intercept) x a
6.9124587 1.6716076 0.3011226
[[2]]
(Intercept) x a
4.2221491 1.0016012 0.2434801
Same can be done for the diagnostics:
gaze.lines.ivreg.diagn <- function(x, col="p-value", row=1:3, digits=2){
stopifnot(is.list(x))
out <- lapply(x, function(y){
stopifnot(class(y)=="summary.ivreg")
y$diagnostics[row, col, drop=FALSE]
})
out <- as.list(data.frame(t(as.data.frame(out)), check.names = FALSE))
for(i in 1:length(out)){
out[[i]] <- c(names(out)[i], round(out[[i]], digits=digits))
}
return(out)
}
gaze.lines.ivreg.diagn(list(summ.fit1, summ.fit2), row=1:2)
gaze.lines.ivreg.diagn(list(summ.fit1, summ.fit2), col=4, row=1:2, digits=2)
Both calls will yield:
$`Weak instruments`
[1] "Weak instruments" "0.45" "0.56"
$`Wu-Hausman`
[1] "Wu-Hausman" "0.11" "0.18"
Now the stargazer() call becomes as simple as this, yielding identical output as above:
stargazer(fit1, fit2, type = "text",
se = gaze.coeft(list(rob.fit1, rob.fit2)),
add.lines = gaze.lines.ivreg.diagn(list(summ.fit1, summ.fit2), row=1:2))
How to add additional GOFs such as ivreg diagnostics to regression tables?
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))
Result 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])))
Result 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))
How to report a standardized model in stargazer package?
You would need to enter the additional values by hand, list-wise for each model, as you started with the coefficients. Standardized se=, the p= values (for the stars), ... as well as the the GOFs (R2, R2adj., ...), read options in help page: ?stargazer.
However, lm.beta appears to add nothing but the standardized coefficients, and none are yet calculated to report them.
Standardized standard errors are calculated using the formula SE*beta_star/beta.
So you could wrap a function, and calculate them, in order to fill them in the stargazer table:
std_se <- \(x) x[, 'Std. Error']*x[, 'Standardized']/x[, 'Estimate']
std_se(summary(mod_std)$coefficients)
# (Intercept) cyl disp
# 0.0000000 0.2109356 0.2109356
However, it might definitely be easier to calculate a actual standardized model
mod_std2 <- lm(mpg ~ cyl + disp, as.data.frame(scale(mtcars)))
summary(mod_std2) |> getElement('coefficients')
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 34.66099474 2.54700388 13.608536 4.022869e-14
# cyl -1.58727681 0.71184427 -2.229809 3.366495e-02
# disp -0.02058363 0.01025748 -2.006696 5.418572e-02
and put that in:
stargazer(mod, mod_std2, type='text')
# ==========================================================
# Dependent variable:
# ----------------------------
# mpg
# (1) (2)
# ----------------------------------------------------------
# cyl -1.587** -0.470**
# (0.712) (0.211)
# disp -0.021* -0.423*
# (0.010) (0.211)
# Constant 34.661*** -0.000
# (2.547) (0.090)
# ----------------------------------------------------------
# Observations 32 32
# R2 0.760 0.760
# Adjusted R2 0.743 0.743
# Residual Std. Error (df = 29) 3.055 0.507
# F Statistic (df = 2; 29) 45.808*** 45.808***
# ==========================================================
# Note: *p<0.1; **p<0.05; ***p<0.01
Related Topics
Initialize an Empty Tibble with Column Names and 0 Rows
Checking Cran Incoming Feasibility ... Note Maintainer
Error with Ggplot2 Mapping Variable to Y and Using Stat="Bin"
Format a Date Column in a Data Frame
How to Neatly Clean My R Workspace While Preserving Certain Objects
Combining Elements of List of Lists by Index
R: How to Display Clustered Matrix Heatmap (Similar Color Patterns Are Grouped)
How to Add an External Legend to Ggpairs()
Split a String Vector at Whitespace
Forcing R Output to Be Scientific Notation with at Most Two Decimals
R: Find and Add Missing (/Non Existing) Rows in Time Related Data Frame
Package 'Stringi' Does Not Work After Updating to R3.2.1
How to Read CSV Data with Unknown Encoding in R
Convert a Matrix with Dimnames into a Long Format Data.Frame