基于mice插补数据用rms包lrm函数合并OR并绘图求助

核心解决思路 #

针对长格式多插补数据，采用遍历单个插补模型+结果合并的方式替代直接调用fit.mult.impute（规避自定义变换的属性传递报错），最终生成合并后的OR曲线。

步骤1：拟合单个插补模型并提取OR值 #

library(mice)
library(rms)
library(dplyr)
library(tidyr)
library(ggplot2)

# 重新初始化datadist（覆盖长格式所有变量）
dd <- datadist(imputed_df_l)
options(datadist = "dd")

# 定义统一模型公式
model_formula <- y ~ child_sex + bmi + gTrans(tiktok_ban, h)

# 生成tiktok_ban的取值序列（覆盖1-90，重点包含跳跃点20）
tiktok_seq <- seq(1, 90, by = 1)

# 固定协变量为参考水平（保证OR解释一致）
ref_covariates <- data.frame(
  child_sex = "Female",
  bmi = mean(imputed_df_l$bmi, na.rm = TRUE)
)

# 遍历所有插补数据集，计算OR
pred_results <- list()
for (imp in 1:20) {
  # 拟合单个插补数据集的lrm模型
  single_model <- lrm(model_formula, data = subset(imputed_df_l, .imp == imp))
  
  # 预测OR（以tiktok_ban=1为基准，转换为指数形式）
  pred_or <- predict(single_model,
                     newdata = cbind(ref_covariates, tiktok_ban = tiktok_seq),
                     type = "fitted",
                     fun = exp)
  
  pred_results[[imp]] <- data.frame(
    .imp = imp,
    tiktok_ban = tiktok_seq,
    OR = pred_or
  )
}

# 合并所有插补结果
combined_pred <- bind_rows(pred_results)

步骤2：计算合并后的OR及置信区间 #

# 按tiktok_ban分组，计算均值与95%分位数置信区间
final_or <- combined_pred %>%
  group_by(tiktok_ban) %>%
  summarise(
    OR_mean = mean(OR),
    OR_lower = quantile(OR, 0.025),
    OR_upper = quantile(OR, 0.975)
  )

步骤3：绘制OR与tiktok_ban的关系图 #

ggplot(final_or, aes(x = tiktok_ban, y = OR_mean)) +
  geom_line(linewidth = 1, color = "#2c3e50") +
  geom_ribbon(aes(ymin = OR_lower, ymax = OR_upper), alpha = 0.2, fill = "#3498db") +
  # 标注跳跃点20天
  geom_vline(xintercept = 20, linetype = "dashed", color = "#e74c3c", linewidth = 1) +
  labs(
    x = "tiktok_ban 天数",
    y = "比值比 (OR)",
    title = "tiktok_ban与事件y的OR关系（多插补合并结果）",
    subtitle = "红色虚线为20天跳跃点"
  ) +
  theme_minimal() +
  scale_y_continuous(limits = c(0, max(final_or$OR_upper) + 0.5))

关键注意事项 #

• 若尝试fit.mult.impute直接拟合报错，本质是自定义变换函数h的nonlinear属性无法在多插补拟合中正确传递，遍历单个模型是更稳定的替代方案。
• 固定协变量时需选择有意义的参考值（如分类变量取基准水平，连续变量取均值），确保OR的可比性。
• 多插补结果合并采用分位数置信区间，符合多插补结果的统计推断规范。

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