Can I use modelsummary with vcov="robust" for cox model in R?

Question

I have a few cox models created by coxph, and I want to summarise them with modelsummary. However, when I write robust=T in coxph script, and run modelsummary like below, I get an error message which says "Error in get_vcov.default(model, vcov = vcov, conf_level = conf_level, : Unable to extract a variance-covariance matrix of type robust from model of class coxph. The variance-covariance matrix is required to adjust the standard errors. The vcov argument accepts a variance-covariance matrix, a vector of standard errors, or a function that returns one of these, such as stats::vcov."

modelsummary(models,
             stars=T,
             vcov = 'robust'      
)

Also, when I use vcov = 'classical', I get the table but the stars disappear. They appear just when I write robust=F in my coxph script.

How can I use robust=T in coxph script and vcov = 'robust' in modelsummary together?

Try using the ‘coxrobust’ R-package – Ventrilocus Aug 28 '22 at 08:05 — Ventrilocus, Aug 28 '22 at 08:05

score 1 · Answer 1 · answered Aug 28 '22 at 18:02

The modelsummary function uses the sandwich package to compute robust standard errors. The vcov="robust" option calls sandwich::vcovHC() to compute heteroskedasticity-consistent standard error. It looks like this function does not support (at least some) coxph models.

However, some other functions from the sandwich package may support your model. For example, here is an example of a model where vcovHC does not work but vcovHAC does work (heteroskedasticity and autocorrelation-consistent standard errors).

library(modelsummary)
library(survival)
library(sandwich)

bladder1 <- bladder[bladder$enum < 5, ]

mod <- coxph(
    Surv(stop, event) ~ rx + size + number,
    cluster = id,
    data = bladder1,
    robust = TRUE)

sandwich::vcovHC(mod)
#> Error in meatHC(x, type = type, omega = omega): hatvalues() could not be extracted successfully but are needed for HC3

sandwich::vcovHAC(mod)
#>                rx        size      number
#> rx     0.49295246 0.021777986 0.020566583
#> size   0.02177799 0.022118330 0.002851453
#> number 0.02056658 0.002851453 0.020215100

modelsummary(mod, vcov = "HAC", output = "markdown")

	Model 1
rx	-0.540
	(0.200)
size	-0.055
	(0.070)
number	0.193
	(0.046)
Num.Obs.	340
AIC	1182.2
BIC	1193.7
RMSE	0.57
Std.Errors	HAC

Again, that may not work for all models. The key point to understand is that the model needs to be supported by sandwich for robust standard errors to be computed by modelsummary.

Can I use modelsummary with vcov="robust" for cox model in R?

1 Answers1