我正在对一个临床队列数据中的死亡原因进行插补。插补过程运行正常。然而，当我尝试使用Fine-Gray模型来估计风险比时，我得到一个自由度，这导致置信区间变得极端宽大，甚至趋向无穷大！！
不过，当我改用CoxPH模型时，这个问题似乎就没有出现，这意味着问题可能出在CRR模型的合并方式上。

下面给出一个示例，使用来自 tidycmprsk 包的 trial 数据框来复现这个问题。

x <- tidycmprsk::trial

set.seed(123) # Set seed for reproducibility
# 1. Randomly select the indices to be turned into NA
na_indices <- sample(which(x$death_cr %in% c("death from cancer", "death other causes")), 24) 

# 2. Create the new variable
x$death_cr_new <- x$death_cr
x$death_cr_new[na_indices] <- NA

# 3. Verification check
print(table(x$death_cr, useNA = "always"))
print(table(x$death_cr_new, useNA = "always"))

# 4. Filter for only those who died # We exclude 'Alive' because we don't want the model to learn from them
x1 <- x %>%
    filter(death_cr_new != "censor" | is.na(death_cr_new)) %>%
    #5. Set up predictors for imputation. I didn't use age and response because they have missings
    select(trt, stage, grade, ttdeath, death_cr_new) %>% 
    droplevels()

# 5. imputation process
init <- mice(x1, maxit = 0)
meth <- init$method # imputing death_cr_new using a logistic regression
imp <- mice(x1, method = meth, m = 20, seed = 123)

# 6. Merging the imputed cause of death data with the censorded cases
x_imputed <- map(1:20, ~{
    x2 <- x %>% filter(death_cr_new == "censor")
    x3 <- complete(imp, .x)

    bind_rows(x2, x3) %>%
        mutate(death_cr_new = factor(death_cr_new, level = c("censor","death from cancer", "death other causes")))
})

# 7. Running the fine-grey competing risks model
models_list <- map(x_imputed, ~ {
    crr(Surv(ttdeath, death_cr_new) ~ trt, 
        data = .x,
        failcode = "death from cancer",
        censor_codes = "censor") # Explicitly define the censoring string
})


mira_models <- as.mira(models_list) # 8. Convert the list of models into a 'mira' object so mice can pool it
pooled_model <- pool(mira_models) # 9. Pool the results using Rubin's Rules
summary(pooled_model, exponentiate = TRUE, conf.int = TRUE) # 10. pooled estimes show extremly wide CI

并给出最终的估计结果，显示自由度为0.4974，置信区间极其宽。

term estimate std.error statistic     df p.value     2.5 %    97.5 %  conf.low conf.high
trtDrug B    1.718    0.2959     1.829 0.4318  0.4974 2.542e-50 1.162e+50 2.542e-50 1.162e+50

看来在合并时没有正确地从每个模型导入自由度。
谢谢

解决方案

手动应用鲁宾规则。

q <- do.call('cbind', lapply(models_list, \(m) m$coef))
u <- do.call('cbind', lapply(models_list, \(m) diag(m$var)))

m <- length(models_list)
qbar <- rowMeans(q)
ubar <- rowMeans(u)
b <- rowSums((q - qbar)^2)/(m - 1L)
t <- ubar + (1 + 1/m)*b
se <- sqrt(t)

out <- data.frame(
  estimate=exp(qbar),
  se=se,
  conf_low=exp(qbar - 1.96*se),
  conf_high=exp(qbar + 1.96*se),
  p=2*pnorm(abs(qbar/se), lower.tail=FALSE)
)

得到

> out
          estimate        se  conf_low conf_high          p
trtDrug B 1.718318 0.2959237 0.9620735  3.069014 0.06734821

注：Death Status 是一个三水平的无序因子。

> x$death_cr_new |> str()
 Factor w/ 3 levels "censor","death from cancer",..: 1 1 1 3 3 2 1 3 1 2 ...
 - attr(*, "label")= chr "Death Status"

插补应使用 meth["death_cr_new"] <- "polyreg"。"logreg" 适用于二元结局，在这里不合适。