多插补与CBPS后Firth逻辑回归的pool函数适配问题求助

核心原因 #

pool()函数（来自mice包）依赖broom包的tidy()方法提取模型结果，但logistf对象没有内置的tidy适配方法，导致报错。以下是三种可行的解决路径：

方法一：手动提取结果并应用Rubin规则合并 #

直接从每个插补的logistf模型中提取系数、标准误等信息，手动调用pool.scalar()实现Rubin规则合并：

# 提取所有插补数据集的logistf模型
models <- weighted_pool_firth$analyses

# 从每个模型中提取系数、标准误、自由度
estimates <- lapply(models, function(x) coef(x))
ses <- lapply(models, function(x) sqrt(diag(vcov(x))))
df_resid <- models[[1]]$df.residual  # 假设各插补数据集自由度一致

# 逐个系数应用Rubin规则合并
pooled_list <- lapply(seq_along(estimates[[1]]), function(i) {
  est_vec <- sapply(estimates, `[`, i)
  se_vec <- sapply(ses, `[`, i)
  pool.scalar(est_vec, se_vec, df_resid)
})

# 整理为标准结果数据框
names(pooled_list) <- names(estimates[[1]])
pooled_results <- do.call(rbind, pooled_list)
pooled_df <- cbind(
  term = rownames(pooled_results),
  estimate = pooled_results[, "qbar"],
  std.error = sqrt(pooled_results[, "t"]),
  df = pooled_results[, "df"],
  conf.low = pooled_results[, "qbar"] - qt(0.975, pooled_results[, "df"]) * sqrt(pooled_results[, "t"]),
  conf.high = pooled_results[, "qbar"] + qt(0.975, pooled_results[, "df"]) * sqrt(pooled_results[, "t"]),
  p.value = sapply(models, function(x) x$prob[i])  # 提取logistf的p值
)
rownames(pooled_df) <- NULL

方法二：为logistf对象自定义tidy方法 #

给logistf类编写适配broom的tidy方法，让pool()函数能正常识别：

library(broom)
library(logistf)

# 自定义tidy方法
tidy.logistf <- function(x, conf.int = FALSE, conf.level = 0.95, ...) {
  res <- data.frame(
    term = names(coef(x)),
    estimate = coef(x),
    std.error = sqrt(diag(vcov(x))),
    statistic = coef(x) / sqrt(diag(vcov(x))),
    p.value = x$prob,
    df = x$df.residual,
    stringsAsFactors = FALSE
  )
  
  if (conf.int) {
    crit_val <- qt((1 + conf.level)/2, x$df.residual)
    res$conf.low <- res$estimate - crit_val * res$std.error
    res$conf.high <- res$estimate + crit_val * res$std.error
  }
  return(res)
}

# 现在可以正常运行pool()
weighted_pool_results <- pool(weighted_pool_firth)
summary(weighted_pool_results, conf.int = TRUE)

方法三：替代方案——用brms实现兼容多重插补的Firth惩罚回归 #

brms支持Firth惩罚，且原生兼容mice的多重插补对象，无需额外适配：

library(brms)

# 拟合每个插补数据集的Firth惩罚logistic回归
weighted_pool_firth_brms <- with(weighted_pool_data,
  brm(
    outcome ~ mainpredictor + abunchofcovariates,
    family = bernoulli(link = "logit"),
    prior = prior(firth(), class = b),  # 启用Firth惩罚
    chains = 2, iter = 2000, warmup = 1000, cores = 2  # 按需调整计算参数
  )
)

# 合并多重插补结果
pooled_brms <- pool(weighted_pool_firth_brms)
summary(pooled_brms)

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