如何从glmmTMB对象中获取整合随机效应后的估算期望值？

你需要的是整合了随机效应的边际期望值（即给定固定效应设计矩阵X时Y的期望，随机效应已被积分出去），glmmTMB的predict()函数默认返回的是包含随机效应预测值的条件预测，re.form=NA则是将随机效应设为0的固定效应预测，这两者都不符合你的需求。以下是针对不同分布的解决方法：

有闭合形式的分布（如泊松/负二项（对数连接）、高斯（恒等连接）） #

泊松/负二项模型（对数连接） #

这类模型的边际期望可以通过解析计算得到：
边际期望 = exp(Xβ) * E[exp(Zb)]
其中：

• Xβ是固定效应部分的线性预测值（可通过predict(fit, re.form=NA, type="link")获取）
• b是服从多元正态分布的随机效应，E[exp(Zb)]的计算利用正态分布的性质：若Zb ~ N(0, V)，则E[exp(Zb)] = exp(Var(Zb)/2)

具体代码示例：

library(glmmTMB)
# 拟合泊松混合模型
fit <- glmmTMB(count ~ mined + (1|site), data = Salamanders, family = poisson)

# 提取固定效应线性预测值（Xβ）
fixef_link <- predict(fit, re.form = NA, type = "link")

# 提取随机效应协方差矩阵
vc <- VarCorr(fit)$cond$site
# 构建随机效应设计矩阵Z（此处为随机截距，对应每个观测的site指示变量）
Z <- model.matrix(~0 + site, data = Salamanders)
# 计算每个观测的Zb的方差
var_Zb <- diag(Z %*% vc %*% t(Z))

# 计算边际期望值
marginal_exp <- exp(fixef_link) * exp(var_Zb / 2)

高斯模型（恒等连接） #

高斯模型的边际期望值直接等于固定效应的预测值，此时用predict(fit, re.form=NA, type="response")即可得到正确结果，因为随机效应的期望为0，积分后不影响最终期望。

无闭合形式的分布（如二项（logit连接）、Gamma（对数连接）） #

这类模型无法通过解析计算得到边际期望，需要用蒙特卡洛模拟近似：
核心思路是多次模拟随机效应的取值，计算对应的条件期望后取平均。

以二项logit模型为例：

library(glmmTMB)
library(MASS)

# 拟合二项混合模型（示例数据需自行替换）
fit <- glmmTMB(cbind(success, failure) ~ treatment + (1|clinic), data = my_data, family = binomial)

# 提取固定效应系数和随机效应协方差
beta <- fixef(fit)$cond
sigma <- VarCorr(fit)$cond$clinic
# 固定效应线性预测值
fixef_link <- predict(fit, re.form = NA, type = "link")
# 随机效应设计矩阵Z
Z <- model.matrix(~0 + clinic, data = my_data)

# 蒙特卡洛模拟随机效应（模拟次数越多结果越准确）
n_sim <- 1000
b_sim <- mvrnorm(n_sim, mu = rep(0, nrow(sigma)), Sigma = sigma)

# 计算每次模拟的条件期望，再取平均得到边际期望
linear_pred <- matrix(fixef_link, nrow = nrow(my_data), ncol = n_sim) + Z %*% t(b_sim)
cond_exp <- plogis(linear_pred)  # logit连接的逆函数
marginal_exp <- rowMeans(cond_exp)

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