ggeffects与marginaleffects预测概率标准误差差异问询

我尝试绘制欧洲社会调查（ESS）西班牙样本中，抑郁症状与抗议参与关系的预测概率图。使用survey包的svyglm拟合了包含抑郁症状二次项的加权调查逻辑回归模型——由于早期ESS轮次没有PSU/分层信息，采用ids=~1的简单抽样设计与后分层权重。运行代码后发现，ggeffects和marginaleffects两个包返回的预测概率标准误差存在细微差异。

以下是数据模拟、模型拟合及结果对比的完整代码：

1. 数据模拟与模型拟合 #

library(survey)
library(ggeffects)
library(marginaleffects)
library(dplyr)
library(ggplot2)

rm(list = ls())
set.seed(123)

# 模拟ESS数据
n <- 6000
df_full <- tibble(
  cntry       = sample(c("ES","FR","DE"), n, replace = TRUE, prob = c(0.4, 0.3, 0.3)),
  essround    = factor(sample(c("3","6","7"), n, replace = TRUE), levels = c("3","6","7")),
  depmeans_rs = runif(n, 0, 1),  # 抑郁症状变量
  pspwght     = runif(n, 0.5, 2.0)  # 调查权重
) %>%
  mutate(
    # 线性预测项：包含二次效应+轮次差异
    lp = -1.5 + 0.9*depmeans_rs - 0.6*(depmeans_rs^2) +
      if_else(essround == "6", 0.25, 0) +
      if_else(essround == "7", 0.15, 0),
    p = plogis(lp),  # 转换为概率
    protest = factor(rbinom(n(), 1, p), levels = c(0, 1))  # 抗议参与因变量
  )

# 构建调查设计对象
df_design <- svydesign(ids = ~1, weights = ~pspwght, data = df_full)

# 筛选西班牙样本
df_design_es <- subset(df_design, cntry == "ES")

# 拟合含二次项的svyglm模型
m2_es <- svyglm(
  protest ~ depmeans_rs + I(depmeans_rs^2) + essround,
  design = df_design_es,
  family = quasibinomial()
)

2. 生成两种方法的预测值 #

# 使用ggeffects生成响应尺度的边际预测值
gge <- predict_response(
  m2_es,
  terms = c("depmeans_rs [0:1 by=.10]", "essround")
)

# 使用marginaleffects生成相同网格的预测值
dep_seq <- seq(0, 1, by = 0.10)
grid_typical <- datagrid(
  model = m2_es,
  depmeans_rs = dep_seq,
  essround = levels(model.frame(m2_es)$essround)
)

me <- predictions(
  m2_es,
  newdata = grid_typical,
  type = "response"
)

# 整理对比用的数据框
gge_cmp <- as.data.frame(gge) %>%
  transmute(
    essround    = group,
    depmeans_rs = x,
    estimate    = predicted,
    conf.low, conf.high,
    method = "ggeffects"
  )

me_cmp <- as.data.frame(me) %>%
  transmute(
    essround, depmeans_rs,
    estimate, conf.low, conf.high,
    method = "marginaleffects"
  )

3. 对比置信区间差异 #

# 合并数据并计算置信区间不对称性
compare_ci <- bind_rows(gge_cmp, me_cmp) %>%
  mutate(
    lower_width = estimate - conf.low,
    upper_width = conf.high - estimate,
    asymmetry   = upper_width - lower_width
  )

# 汇总两种方法的不对称性统计量
compare_ci %>%
  group_by(method) %>%
  summarise(
    mean_asymmetry = mean(asymmetry),
    max_abs_asym   = max(abs(asymmetry))
  )

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