如何在gtsummary包的tbl_uvregression中添加聚类标准误？

可以通过结合sandwich和lmtest包，为tbl_uvregression输出的单变量逻辑回归结果替换聚类稳健标准误、置信区间和p值，具体实现步骤如下：

1. 加载所需包 #

library(gtsummary)
library(sandwich)
library(lmtest)

2. 构造示例数据 #

df <- data.frame(
  cluster = c(1,2,3,4,5),
  y = c(0,0,1,1,1),
  x1 = c(10,20,30,40,50),
  x2 = c(0,1,0,1,1)
)

3. 生成初始单变量回归表 #

先运行基础的tbl_uvregression得到初始结果：

x <- tbl_uvregression(
  df[c("x1", "x2", "y", "cluster")],
  method = glm,
  y = y,
  method.args = list(family = binomial),
  include = -cluster,
  exponentiate = TRUE
)

4. 替换为聚类稳健标准误 #

使用modify_table_body函数，针对每个变量的模型计算聚类稳健的统计量，替换原表中的对应值：

x_cluster <- x %>%
  modify_table_body(
    ~ .x %>%
      dplyr::rowwise() %>%
      dplyr::mutate(
        # 重新拟合单变量模型
        model = list(glm(formula = formula, data = df, family = binomial)),
        # 计算聚类稳健方差协方差矩阵
        vcov_cl = list(vcovCL(model, cluster = ~cluster)),
        # 提取聚类稳健的系数、标准误、p值
        coef_test = list(coeftest(model, vcov = vcov_cl)),
        estimate = exp(coef_test[2, 1]),
        std.error = coef_test[2, 2],
        p.value = coef_test[2, 4],
        # 计算聚类稳健置信区间
        conf.low = exp(coef_test[2, 1] - qnorm(0.975) * coef_test[2, 2]),
        conf.high = exp(coef_test[2, 1] + qnorm(0.975) * coef_test[2, 2])
      ) %>%
      dplyr::ungroup() %>%
      dplyr::select(-model, -vcov_cl, -coef_test)
  )

# 查看结果
x_cluster

说明 #

• 针对每个单变量模型重新拟合后，用vcovCL()指定聚类变量计算稳健方差协方差矩阵；
• 通过coeftest()提取调整后的标准误和p值，因已设置exponentiate = TRUE，需对对数尺度的置信区间取指数得到最终结果；
• 最终输出的表格将展示聚类稳健的标准误、置信区间和p值。

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