Stargazer dosen't show number of observations for robust s.e

Question

I'm trying to use stargazer to output some regression results with robust standard errors, but the lines in the bottom for values like nobs and f stat doesn't show. here are the code and the output.

r_FEB1_se <- coeftest(r_FEB1, vcov = vcovHC(r_FEB1, method = "arellano", type = "HC3", cluster = "group"))
r_FED1_se <- coeftest(r_FED1, vcov = vcovHC(r_FED1, method = "arellano", type = "HC3", cluster = "group"))
r_FEF1_se <- coeftest(r_FEF1, vcov = vcovHC(r_FEF1, method = "arellano", type = "HC3", cluster = "group"))
r_FEH1_se <- coeftest(r_FEH1, vcov = vcovHC(r_FEH1, method = "arellano", type = "HC3", cluster = "group"))

r_FEB2_se <- coeftest(r_FEB2, vcov = vcovHC(r_FEB2, method = "arellano", type = "HC3", cluster = "group"))
r_FED2_se <- coeftest(r_FED2, vcov = vcovHC(r_FED2, method = "arellano", type = "HC3", cluster = "group"))
r_FEF2_se <- coeftest(r_FEF2, vcov = vcovHC(r_FEF2, method = "arellano", type = "HC3", cluster = "group"))
r_FEH2_se <- coeftest(r_FEH2, vcov = vcovHC(r_FEH2, method = "arellano", type = "HC3", cluster = "group"))

stargazer::stargazer(r_FEB1_se,r_FEB2_se,r_FED1_se,r_FED2_se,r_FEF1_se,r_FEF2_se,r_FEH1_se,r_FEH2_se, type = "text")

=========================================================================================================================
                                                                          Dependent variable:                            
                                               --------------------------------------------------------------------------

                                                 (1)      (2)     (3)      (4)      (5)      (6)       (7)        (8)    
-------------------------------------------------------------------------------------------------------------------------
`Ratio Immigrants/pop t-1`                      0.594*  0.653*  1.193*** 1.368*** 0.973*** 1.031***  0.841***   0.904*** 
                                               (0.337)  (0.338) (0.450)  (0.458)  (0.295)  (0.306)   (0.294)    (0.299)  

`Disposable income in Thousand EUR t-1`         0.000*           0.000             0.000            0.00000***           
                                               (0.000)          (0.000)           (0.000)            (0.000)             

`Unemployment rate on all civilian income t-1` 0.010***         0.031***           -0.006            0.009**             
                                               (0.003)          (0.007)           (0.005)            (0.004)             

`Disposable income in Thousand EUR t-2`                  0.000            0.000             0.000              0.00000***
                                                        (0.000)          (0.000)           (0.000)              (0.000)  

`Unemployment rate on all civilian income t-2`          0.008**          0.025***          -0.011**              0.007   
                                                        (0.003)          (0.007)           (0.005)              (0.005)  

=========================================================================================================================
=========================================================================================================================
Note:                                                                                         *p<0.1; **p<0.05; ***p<0.01

Answer 1

Usually the "lm" objects are used to produce the tables. But you may override the values. Since I'm more familiar with the texreg package, I hope you don't mind that I provide an according solution.

You need both, the "lm" objects, denoted as fit* below, and the matrices with robust standard errors resulting from coeftest, denoted by rob.fit* below. Example:

fit1 <- lm(mpg ~ hp, mtcars)
fit2 <- lm(mpg ~ hp + am, mtcars)
library(sandwich);library(lmtest)
rob.fit1 <- coeftest(fit1, vcov.=vcovHC(fit1))
rob.fit2 <- coeftest(fit2, vcov.=vcovHC(fit2))

Then, in texreg use override.* to provide modified standard errors, as well as the p-values.

library(texreg)
screenreg(list(fit1, fit2),
          override.se=list(rob.fit1[,2], rob.fit2[,2]),
          override.pvalues=list(rob.fit1[,4], rob.fit2[,4]))

# =================================
#              Model 1    Model 2  
# ---------------------------------
# (Intercept)  30.10 ***  26.58 ***
#              (2.41)     (1.47)   
# hp           -0.07 ***  -0.06 ***
#              (0.02)     (0.01)   
# am                       5.28 ***
#                         (1.22)   
# ---------------------------------
# R^2           0.60       0.78    
# Adj. R^2      0.59       0.77    
# Num. obs.    32         32                        <-- THERE THEY ARE!!
# RMSE          3.86       2.91    
# =================================
# *** p < 0.001, ** p < 0.01, * p < 0.05

Just replace screenreg with texreg for LaTeX, or htmlreg for HTML, whatever you want.

I'm sure this might be somehow also possible with stargazer.

Answer 2

I agree with jay.sf that stargazer will not show test statistics if coeftest objects instead of regression models are supplied. But there is a solution without the need to switch to texreg.

You can still use coeftest to find robust standard errors. You need to save the standard errors from each model in one numeric vector, and then combine all vectors into a single list, which will be suppplied to the se option in stargazer. But be careful that stargazer must rely on the element names in those vectors to match each robust standard error to each independent variable. So, for each model, you should subset the coeftest object by using [,2], as the second column is for robust standard errors. The variable names will be retained in this way. Then, you can merge all numeric vectors into a list and supply it into the se option.

A sample code that does not fit into your dataset:

(I have eight models to run, from q4a1 to q4a8, and they are stored in the list q4a. I want to find two-way clustered standard errors instead of robust standard errors.)

se.q4a <- sapply(q4a, function(z) coeftest(z, vcov = vcovCL, cluster = ~ wbcode + cluster)[, 2])
stargazer(q4a[[1]], q4a[[2]], q4a[[3]], q4a[[4]], q4a[[5]], q4a[[6]], q4a[[7]], q4a[[8]], keep = c("rlaw", "corrupt"), type = "text", se = se.q4a)