R中多重插补数据执行exp(lincom)计算GEE模型百分比变化报错求助

lincom仅支持单个统计模型的线性组合计算，无法直接适配mira（多重插补模型集合）和mipo（合并后的插补结果）对象，以下是针对两种对象的解决方案：

一、处理mira类型对象 #

mira本质是多个插补后GEE模型的列表，需逐个计算每个模型的结果再合并：

library(mice)
library(geepack)
library(lmtest)

# 遍历mira中的每个GEE模型，计算v1的百分比变化及标准误
mira_calc <- lapply(gee_imputed_data$analyses, function(mod) {
  lc_result <- lincom(mod, "v1")
  perc_change <- (exp(lc_result$estimate) - 1) * 100
  # 用链式法则计算百分比变化的标准误
  perc_se <- lc_result$se * exp(lc_result$estimate)
  data.frame(estimate = perc_change, se = perc_se, df = lc_result$df)
})

# 合并多个插补结果
pooled_result <- pool(mira_calc)
# 查看合并后的统计量
summary(pooled_result)

二、处理mipo类型对象 #

mipo已经是合并后的系数结果，直接提取对应变量的统计量转换即可：

# 提取v1的合并系数、标准误和置信区间
v1_coef <- pool.data$estimate[pool.data$term == "v1"]
v1_se <- pool.data$std.error[pool.data$term == "v1"]
v1_ci_low <- pool.data$conf.low[pool.data$term == "v1"]
v1_ci_high <- pool.data$conf.high[pool.data$term == "v1"]

# 转换为百分比变化
perc_change <- (exp(v1_coef) - 1) * 100
perc_se <- v1_se * exp(v1_coef)
perc_ci_low <- (exp(v1_ci_low) - 1) * 100
perc_ci_high <- (exp(v1_ci_high) - 1) * 100

# 输出结果
cat(sprintf("v1的百分比变化：%.2f%%，标准误：%.2f\n95%%置信区间：(%.2f%%, %.2f%%)\n",
            perc_change, perc_se, perc_ci_low, perc_ci_high))

报错原因说明 #

lincom需要调用模型的vcov方法提取方差-协方差矩阵，但mira和mipo并非单个模型对象，没有对应vcov实现，因此直接调用会触发方法不匹配的错误。

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