R modelsummary output: robust standard errors from coxph model

Question

Im running a series of coxph models in R and compiling the output into latex tables using the modelsummary package and command. The coxph provides SE and robust se as outputs and the p-value is based on the robust se. Here is a quick example to illustrate the output:

test_data <- list(time=c(4,3,1,1,2,2,3)
                  , status=c(1,1,1,0,1,1,0)
                  , x=c(0,2,1,1,1,0,0)
                  , sex=c(0,0,0,0,1,1,1)) 
model <- coxph(Surv(time, status) ~ x + cluster(sex), test_data)

Call:
coxph(formula = Surv(time, status) ~ x, data = test1, cluster = sex)

      coef exp(coef) se(coef) robust se     z      p
x 0.460778  1.585306 0.562800  0.001052 437.9 <2e-16

Likelihood ratio test=0.66  on 1 df, p=0.4176
n= 7, number of events= 5

Next, Im trying to create latex tables from these models, displaying the robust se.

modelsummary(model, output = "markdown", fmt = 3, estimate = "{estimate}{stars}", statistic = "std.error")

|                      |  Model 1   |
|:---------------------|:----------:|
|x                     |  0.461***  |
|                      |  (0.563)   |
|Num.Obs.              |     7      |
|R2                    |   0.090    |
|AIC                   |    13.4    |

As we can see, only the non-adjusted se is displayed. I could not find any alternative for this statistic = "std.error" parameter that fits and also something like vcov="robust" does not work.

How can I display any kind of robust standard errors using modelsummary for coxph models?

Thanks for reading and any help is appreciated.

Answer 1

The next version of modelsummary (now available on Github) will produce a more informative error message in those cases:

library(survival)
library(modelsummary)
test_data <- list(time=c(4,3,1,1,2,2,3)
                  , status=c(1,1,1,0,1,1,0)
                  , x=c(0,2,1,1,1,0,0)
                  , sex=c(0,0,0,0,1,1,1)) 
model <- coxph(Surv(time, status) ~ x + cluster(sex), test_data)
modelsummary(model, vcov = "robust")

# Error: Unable to extract a variance-covariance matrix for model object of class
# `coxph`. Different values of the `vcov` argument trigger calls to the `sandwich`
# or `clubSandwich` packages in order to extract the matrix (see
# `?insight::get_varcov`). Your model or the requested estimation type may not be
# supported by one or both of those packages, or you were missing one or more
# required arguments in `vcov_args` (like `cluster`).

As you can see, the problem is that modelsummary does not compute robust standard errors itself. Instead, it delegates this task to the sandwich or clubSandwich packages. Unfortunately, this coxph model does not appear appear to be supported by those packages:

sandwich::vcovHC(model)

#> Error in apply(abs(ef) < .Machine$double.eps, 1L, all): dim(X) must have a positive length

sandwich is the main package in the R ecosystem to compute robust standard errors. AFAICT, all the other table-making packages available (e.g., stargazer, texreg) also use sandwich, so you are unlikely to have success by looking at those. If you find another package which can compute robust standard errors for Cox models, please file a report on the modelsummary Github repository. I will investigate to see if it’s possible to add support then.

If the info you want is available in the summary object, you can add this information by following the instructions here:

https://vincentarelbundock.github.io/modelsummary/articles/modelsummary.html#adding-new-information-to-existing-models

tidy_custom.coxph <- function(x, ...) {
    s <- summary(x)$coefficients
    data.frame(
        term = row.names(s),
        robust.se = s[, "robust se", drop = FALSE])
}

modelsummary(model, statistic = "robust.se")

	Model 1
x	0.461
	(0.001)
Num.Obs.	7
AIC	13.4
BIC	13.4
RMSE	0.61