CODEWITHSUNDEEP
rmodelsummary
hobsbawm
hobsbawm

Reputation: 21

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

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?

Upvotes: 0

Views: 917

Answers (1)

Vincent
Vincent

Reputation: 17715

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.

Upvotes: 1

Related Questions