如何将多重插补+PSM生成的Cox回归mimira对象绘制成森林图

尊敬的StackOverflow社区：
我是一名外科医生，已通过StackOverflow等网站自学R语言6个月，若我的问题较为基础，请各位多多包涵。

背景说明 #

我的研究目标是对癌症患者数据集进行Cox生存回归分析。由于是回顾性研究，我计划采用1:3倾向得分匹配（PSM），缺失数据使用mice包做多重插补处理，PSM通过MatchThem包实现。我使用survey包的svycoxph()函数通过with()函数池化生存模型，得到了mimira对象，该对象可以通过gtsummary包的tbl_regression()正常输出美观的结果表格。

遇到的问题 #

我通常会将Cox回归结果输出为风险比表格和森林图（使用survminer包的ggforest()），但本次我遇到了阻碍：ggforest无法识别mimira对象为coxph对象，抛出如下错误：

Error in ggforest(tbl_regression_object, data = mimira_object) : 
  inherits(model, "coxph") is not TRUE 

我推测问题出在多重插补之外新增的PSM步骤，因为之前仅使用多重插补得到的mira对象可以通过pool_and_tidy_mice()函数正常配合ggforest输出森林图。

复现代码 #

# 加载数据处理依赖
library(fabricatr)
library(simsurv)

# 模拟临床试验患者数据
participant_data <- fabricate(
  N = 2000,
  age = runif(N, min = 18, max = 85),
  is_female = draw_binary(prob = 0.5, N = N),
  is_smoker = draw_binary(prob = 0.2 + 0.2 * (age > 50), N = N),
  disease_stage = round(runif(N, min = 1 + 0.5 * (age > 65), max = 4)),
  treatment = draw_binary(prob = 0.5, N = N),
  kps = runif(N, min = 40, max = 100)
)

# 模拟生存结局数据
survival_data <- simsurv(
  lambdas = 0.1, gammas = 1.8,
  x = participant_data, 
  betas = c(is_female = -0.2, is_smoker = 1.2,
            treatment = -0.4, kps = -0.005,
            disease_stage = 0.2),
  maxt = 5)

# 合并数据集
library(dplyr)
mydata_complete <- bind_cols(survival_data, participant_data)

# 生成含缺失值的数据集
library(missMethods)
mydata_uncomp <- delete_MCAR(mydata_complete, 0.3)
mydata <- mydata_uncomp

# 1. 用mice包做多重插补
library(mice)
mydata$nelsonaalen <- nelsonaalen(mydata, eventtime, status)
mydata_mice_imp_m3 <- mice(mydata, maxit = 2, m = 3, seed = 20200801) # m=3为测试用参数


# 2. 用MatchThem包做1:3倾向得分匹配
library(MatchThem)
mydata_imp_m3_psm <- matchthem(treatment ~ age + is_female + disease_stage, data = mydata_mice_imp_m3, approach = "within" ,ratio= 1, method = "optimal")

# 3. 用survey包池化多重插补+匹配后的Cox模型
library(survey)
mimira_object <- with(data = mydata_imp_m3_psm, expr = svycoxph(Surv(eventtime, status) ~ age+ is_smoker + disease_stage))
pool_and_tidy_mice(mimira_object, exponentiate = TRUE, conf.int=TRUE) -> pooled_imp_m3_cph

# 运行上述代码后出现警告：In get.dfcom(object, dfcom) : Infinite sample size assumed.
# pooled_imp_m3_cph输出结果如下：
#            term  estimate   std.error  statistic p.value conf.low conf.high            b  df dfcom fmi    lambda m      riv         ubar
# 1           age 0.9995807 0.001961343 -0.2138208     NaN      NaN       NaN 1.489769e-06 NaN   Inf NaN 0.5163574 3 1.067643 1.860509e-06
# 2     is_smoker 2.8626952 0.093476026 11.2516931     NaN      NaN       NaN 4.182884e-03 NaN   Inf NaN 0.6382842 3 1.764601 3.160589e-03
# 3 disease_stage 1.2386947 0.044092483  4.8547535     NaN      NaN       NaN 8.995628e-04 NaN   Inf NaN 0.6169374 3 1.610540 7.447299e-04

# 4. 用gtsummary包生成汇总表格
library(gtsummary)
tbl_regression_object <- tbl_regression(mimira_object, exp=TRUE, conf.int = TRUE) 
# 上述代码无95%CI和p值输出，原因为mimira对象池化时Matchthem:::get.2dfcom函数返回dfcom = 999999


# 5. 预期输出效果（无插补、无匹配的普通Cox模型）
library(survival)
mydata.cox <- coxph(Surv(eventtime, status) ~ age+ is_smoker + disease_stage, mydata_uncomp)

# 用gtsummary绘制森林图
forestGT <- 
  mydata.cox %>%
  tbl_regression(exponentiate = TRUE,
                 add_estimate_to_reference_rows = TRUE) %>% 
  plot()
# 该输出基本符合预期，期望能补充样本量、95%CI、HR值、p值以及模型参数（AIC、事件数、一致性等）

# 用survminer绘制森林图
HRforest <- 
  survminer::ggforest(mydata.cox, data = mydata_uncomp)
# 该输出包含所有需要的Cox回归相关信息，是目标森林图样式

# 6. 插补+匹配后运行的实际效果
# 用gtsummary绘制森林图
forestGT_imp_psm <- 
  mimira_object %>%
  tbl_regression(exponentiate = TRUE,
                 add_estimate_to_reference_rows = TRUE) %>% 
  plot() # 警告：In get.dfcom(object, dfcom) : Infinite sample size assumed.
# 可输出图形但缺失95%CI

# 用survminer绘制森林图
HRforest_imp_psm <- 
  ggforest(mimira_object, data = mydata_imp_m3_psm) 
# 报错：in ggforest(mimira_object, data = mydata_imp_m3_psm) : inherits(model, "coxph") is not TRUE

非常感谢各位的帮助。
祝好
AK

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