在R中基于100组多重插补数据绘制Logit回归平均边际效应(AME)

我用R的mice包生成了100组多重插补数据，拟合Logit回归后通过Rubin法则得到了合并估计。目前已经能在单数据集上用sjPlot::plot_model()或marginaleffects::plot_slopes()绘制平均边际效应（AME），但不知道如何绘制这100组插补数据合并后的AME。

可复现代码 #

library(mice)
library(marginaleffects)

## 创建分析数据集
tmpdf<-iris
virginica<-ifelse(tmpdf$Species=="virginica", 1,0)
df2<-data.frame(Sepal.Length = tmpdf$Sepal.Length,
                Sepal.Width = tmpdf$Sepal.Width,
                virginica = virginica)

## 生成含缺失值的数据集
missdf<-mice::ampute(df2)$amp

## 生成100组插补数据
impdf<-mice::mice(missdf,m=100,maxit=1,seed=1234)

## 拟合Logit回归并合并结果
model<-with(impdf, 
            glm(virginica ~ Sepal.Length + Sepal.Width, 
                family=binomial("logit")))

## 计算并合并平均边际效应
marginal<-marginaleffects::avg_slopes(model)
marginal

单数据集绘图代码（参考） #

# 无缺失值的完整数据集
tmpdf<-iris
virginica<-ifelse(tmpdf$Species=="virginica", 1,0)
df2<-data.frame(Sepal.Length = tmpdf$Sepal.Length,
                Sepal.Width = tmpdf$Sepal.Width,
                virginica = virginica)

model2<- glm(virginica ~ Sepal.Length + Sepal.Width, 
             family=binomial("logit"), data=df2)

# 绘制Sepal.Length的斜率随Sepal.Width变化的边际效应
plot_slopes(model2,variables = "Sepal.Length",condition="Sepal.Width")

解决方案 #

1. 直接绘制合并后的平均边际效应

marginaleffects包支持直接对多重插补模型的合并边际效应结果绘图，自带的plot()函数会自动展示合并后的估计值与置信区间：

# 基于已计算的marginal对象绘图
plot(marginal)

如果需要自定义图形，可将合并结果转为数据框后用ggplot2绘制：

library(ggplot2)

# 提取合并后的边际效应数据
marginal_df <- as.data.frame(marginal)

# 绘制点图+误差棒
ggplot(marginal_df, aes(x = term, y = estimate)) +
    geom_point(size = 3, color = "#2E86AB") +
    geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.2, color = "#2E86AB") +
    labs(x = "自变量", y = "平均边际效应（AME）", title = "多重插补合并后的平均边际效应") +
    theme_minimal()

2. 绘制合并后的条件边际效应（如斜率随协变量变化）

要绘制类似单数据集plot_slopes()的条件效应，需先对每个插补数据集计算条件斜率，再合并结果后绘图：

library(dplyr)
library(ggplot2)

# 对每个插补数据集计算条件斜率（不直接绘图）
slopes_per_imp <- with(impdf, {
    mod <- glm(virginica ~ Sepal.Length + Sepal.Width, family = binomial("logit"))
    plot_slopes(mod, variables = "Sepal.Length", condition = "Sepal.Width", draw = FALSE)
})

# 合并所有插补数据集的结果，计算均值及95%分位数置信区间
combined_slopes <- bind_rows(lapply(slopes_per_imp$analyses, function(x) x$dat)) %>%
    group_by(condition, condition_value) %>%
    summarize(
        estimate = mean(estimate),
        conf.low = quantile(conf.low, 0.025),
        conf.high = quantile(conf.high, 0.975)
    )

# 绘制合并后的条件边际效应曲线
ggplot(combined_slopes, aes(x = condition_value, y = estimate)) +
    geom_line(color = "#2E86AB", linewidth = 1) +
    geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.2, fill = "#2E86AB") +
    labs(x = "Sepal.Width", y = "Sepal.Length对virginica的边际效应", title = "合并插补数据后的条件边际效应") +
    theme_minimal()

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