Search code examples
rregressionstargazer

Robust Standard Errors in lm() using stargazer() for multiple models


I am using stargazer() to display the results of my OLS including robust standard errors for two models m3adj and m4adj. This is my code:

model.lst <- list(m3adj, m4adj)

stargazer(model.lst,
          title = "Initial Results",
          type = "html",
          align = TRUE,
          no.space=TRUE,
          column.labels = c("Model 1 (Main effect)", "Model 2 (Main effect of sub dimensions)"), 
          omit.stat=c("LL","ser", "f"), 
          ord.intercepts = TRUE,
          notes.append=TRUE,
          se = c(list(NULL,robust_se3, NULL, robust_se4)), #include robust standard errors
          single.row = TRUE, 
          report = "vcs*", 
          notes="Heteroscedasticity-robust standard errors in parantheses.",
          add.lines=list(c('Year-quarter fixed effects', 'Yes','Yes')),
          star.cutoffs = c(0.05, 0.01, 0.001)
)

However, my robust standard errors which are calculated in robust_se3 and robust_se4 only show up partially. I only get the SEs in parantheses for Model 1 but for Model 2 I get random numbers in parantheses and only for a couple of my variables. Any ideas how to fix my code so I can include robust SEs for both models?


Solution

  • Found the answer myself: First calculate the robust standard errors for each model m3 and m4 using the code

    robust_se2 <- sqrt(diag(vcovHC(m2, type = "HC0")))

    and then adapting stargazer as follows:

    model.lst <- list(m3, m4)
    
    stargazer(model.lst,
              title = "Initial Results",
              type = "html",
              align = TRUE,
              no.space=TRUE,
              column.labels = c("Model 1 (Main effect)", "Model 2 (Main effect of sub dimensions)"), 
              omit.stat=c("LL","ser", "f"), 
              ord.intercepts = TRUE,
              notes.append=TRUE,
              se = list(robust_se3, NULL, robust_se4), #include robust standard errors
              single.row = TRUE, 
              report = "vcs*", 
              notes="Heteroscedasticity-robust standard errors in parantheses.",
              add.lines=list(c('Year-quarter fixed effects', 'Yes','Yes')),
              star.cutoffs = c(0.05, 0.01, 0.001)
    )