如何计算广义线性混合模型中交互项的预测均值与标准误？

我用以下代码构建了包含交互项的广义线性混合模型：

be_model1 <- glmer(number_of_fledged_hatchlings ~ species*nestbox_quality + (1|habitat) + (1|year), family = poisson, data = be)

固定效应species包含两个水平：大山雀（GT）、蓝山雀（BT）；nestbox_quality包含两个水平："high"（高质量）、"low"（低质量）。

我需要计算交互项的拟合均值及标准误，比如低质量巢箱中大山雀（GT）的拟合平均出飞雏鸟数量及其标准误，但尝试三种方法后结果不一致，存在困惑：

尝试的方法及问题 #

1. allEffects()方法：仅输出拟合均值，未直接显示标准误

> allEffects(be_model1, se=TRUE)
 model: number_of_fledged_hatchlings ~ species * nestbox_quality

 species* effect
    nestbox_quality
species    high low
  GT 6.275469   6.869725
  BT 6.547347   5.819534

2. ggpredict()方法：输出95%置信区间，但无法直接获取标准误

library(ggeffects)
> (ggpredict(be_model1, c("nestbox_quality", "species")))
# Predicted counts of number_of_fledged_hatchlings

# species = BT

nestbox_quality  | Predicted |       95% CI
-------------------------------------
high             |      6.28 | [5.50, 7.16]
low              |      6.87 | [6.07, 7.77]

# species = GT

nestbox_quality  | Predicted |       95% CI
-------------------------------------
high             |      6.55 | [5.80, 7.40]
low              |      5.82 | [5.21, 6.50]

Adjusted for:
*   habitat = 0 (population-level)
*      year = 0 (population-level)

3. predict()+tapply()方法：得到标准误，但拟合均值和前两种方法差异明显

> pred_model2 <- predict(be_model1, type="response")
> 
> tapply(pred_model2, list(species, nestbox_quality), mean)
         high low
BT 6.457035   6.838857
GT 6.498190   5.792777
> tapply(pred_model2, list(species, nestbox_quality), std.error)
         high low
BT 0.10072955 0.10153544
GT 0.08011282 0.05989635

1. 结果差异的核心原因 #

• allEffects()和ggpredict()：计算的是**群体水平（边际均值）**的拟合值，即把随机效应（habitat和year）设为0，仅考虑固定效应组合的预测结果，是排除随机效应后的平均水平。
• predict()默认参数：计算的是条件水平的预测值，即包含了每个观测对应的具体habitat和year的随机偏移，之后分组求平均是对带随机效应的个体预测值取组内平均，和边际均值的统计意义完全不同，所以结果必然存在差异。

2. 获取边际均值及标准误的正确方法 #

方法1：从allEffects()提取标准误 #

allEffects()的结果对象中已经包含标准误，只是默认打印未显示，可直接提取：

eff <- allEffects(be_model1, se=TRUE)
# 提取低质量巢箱中大山雀的拟合均值和标准误
eff$species:nestbox_quality$fit["GT", "low"]  # 拟合均值
eff$species:nestbox_quality$se["GT", "low"]   # 标准误

方法2：从ggpredict()结果推导标准误 #

泊松模型的95%置信区间基于均值 ± 1.96×标准误（正态近似），可通过区间宽度反推标准误：

gg_res <- ggpredict(be_model1, c("nestbox_quality", "species"))
# 筛选目标组合的行
gt_low_row <- subset(gg_res, species == "GT" & nestbox_quality == "low")
# 计算标准误
se_gt_low <- (gt_low_row$conf.high - gt_low_row$conf.low)/(2*1.96)
# 拟合均值
mean_gt_low <- gt_low_row$predicted

方法3：用emmeans()直接计算（推荐） #

emmeans是计算边际均值的专用工具，支持混合模型，可直接输出均值、标准误和置信区间：

library(emmeans)
# 针对交互项计算边际均值，type="response"转换到原始计数尺度
emm <- emmeans(be_model1, ~ species*nestbox_quality, type="response")
# 查看所有组合的结果
emm
# 单独提取低质量巢箱中大山雀的结果
subset(emm, species == "GT" & nestbox_quality == "low")

输出会直接包含emmean（拟合均值）、SE（标准误）、lower.CL/upper.CL（置信区间），结果和allEffects()/ggpredict()完全一致。

3. 关于predict()方法的说明 #

若你需要的是包含随机效应的条件水平均值，predict()的结果是合理的，但此时的组内平均是带特定habitat和year偏移的条件均值，和前两种方法的边际均值统计意义不同，需根据研究需求选择，避免混淆。

内容的提问来源于stack exchange，提问作者Marle Mül