Can you get robust standard errors using the tbl_regression function in gtsummary?

Question

I have conducted several binominal logistic regression models using the standard glm() function in R and evaluated fit with tbl_regression() from the package gtsummary. I wanted output using robust standard errors similar to vcov(model, robust="HC1") from the sandwich-package.

Conducted the code tbl_regression (model, exponentiate=TRUE, robust="HC1") and got exact same output as when using jtools::summ(model, exp=TRUE, robust="HC!") which is based on sandwich. (no errors or warning-messages)

Does anyone know what method is used to get robust estimates in gtsummary()/ or if output is actually robust? Do not find any documentation that this is a functionality in the package

Grateful for all help!

gtsummary: https://github.com/ddsjoberg/gtsummary#readme

Answer 1

The tbl_regression() function uses the tidy function from the broom package to clean up model results, including the standard error and CI estimates. We can, however, pass a custom tidier in the tbl_regression(tidy_fun=) argument that calculates the robust estimates. In the example below, I wrote a tidier wrapping the jtools::summ() function.

library(gtsummary)
packageVersion("gtsummary")
#> [1] '1.4.2'
# build logistic regression model
mod <- glm(response ~ age + trt, data = trial, family = binomial)

# write a tidy function wrapping jtools::summ
tidy_robust <- function(x, 
                        exponentiate = FALSE, 
                        conf.level = 0.95, 
                        robust = c("HC0", "HC1", "HC2", "HC3", "HC4", "HC4m", "HC5"), 
                        ...) {
  robust <- match.arg(robust)
  
  jtools::summ(x, confint = TRUE, ci.width = conf.level, 
               exp = exponentiate, robust = robust) %>%
    purrr::pluck("coeftable") %>%
    as.data.frame() %>%
    tibble::rownames_to_column() %>%
    setNames(c("term", "estimate", "conf.low", "conf.high", 
               "statistic", "p.value"))
}

tidy_robust(mod)
#>          term   estimate     conf.low  conf.high  statistic     p.value
#> 1 (Intercept) -1.7424131 -2.918252621 -0.5665735 -2.9043646 0.003679993
#> 2         age  0.0189970 -0.003406793  0.0414008  1.6619255 0.096527704
#> 3   trtDrug B  0.1254909 -0.504805356  0.7557872  0.3902255 0.696369799
# use default robust SE, HC0
tbl1 <- 
  tbl_regression(
    mod, 
    exponentiate = TRUE,
    tidy_fun = tidy_robust
  )

# use robust SE, HC1
tbl2 <- 
  tbl_regression(
    mod, 
    exponentiate = TRUE,
    tidy_fun = purrr::partial(tidy_robust, robust = "HC1")
  )

^{Created on 2021-09-18 by the reprex package (v2.0.1)}

Robust SE reporting is something that will be made easier in a future release of gtsummary.