如何用R的modelsummary展示coxph模型的稳健标准误？

当你在coxph中使用cluster()指定聚类变量后，modelsummary默认提取的是普通标准误，要展示聚类稳健标准误，可以通过以下两种方法实现：

方法1：利用vcov参数传入聚类稳健方差矩阵 #

通过sandwich包的vcovCL()函数计算聚类稳健方差矩阵，直接传给modelsummary的vcov参数，工具会自动基于该矩阵计算标准误和对应的稳健p值。

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

# 构建数据与模型
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), data = test_data)

# 生成带稳健标准误的LaTeX表格
modelsummary(
  model,
  output = "latex",
  fmt = 3,
  estimate = "{estimate}{stars}",
  statistic = "std.error",
  vcov = vcovCL(model, cluster = ~sex)  # 指定聚类稳健方差矩阵
)

方法2：手动提取稳健统计量并覆盖模型结果 #

从coxph模型的summary()结果中提取稳健标准误和p值，覆盖tidy()输出的对应字段后，再传给modelsummary：

library(survival)
library(modelsummary)
library(broom)

# 构建模型（同上）
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), data = test_data)

# 提取模型的tidy结果
tidy_model <- tidy(model, conf.int = TRUE)
# 从summary中提取稳健标准误和p值
sum_model <- summary(model)
tidy_model$std.error <- sum_model$coefficients[, "robust se"]
tidy_model$p.value <- sum_model$coefficients[, "p"]

# 生成表格
modelsummary(
  tidy_model,
  output = "latex",
  fmt = 3,
  estimate = "{estimate}{stars}",
  statistic = "std.error"
)

关键说明 #

• 直接使用vcov="robust"无效，因为该参数默认对应普通Huber-White稳健方差，而非聚类稳健方差，必须通过vcovCL()指定聚类变量来匹配你的cluster(sex)设置。
• 两种方法都会输出基于稳健标准误的p值，符合需求。

内容的提问来源于stack exchange，提问作者dmort