R中结合Location因素从模拟后验分布计算拟合值

1. 构造包含所有变量组合的预测数据集 #

首先生成覆盖sex所有水平、size连续范围、location所有水平的全交叉数据集，确保覆盖所有需要预测的场景：

# 假设原始数据框为df，fit是lme4拟合的LMM对象
new_dat <- expand.grid(
  sex = unique(df$sex),  # 保留原数据中sex的全部水平
  location = unique(df$location),  # 保留原数据中location的全部水平
  size = seq(min(df$size), max(df$size), length.out = 50)  # 生成平滑的size连续序列
)

# 确保分类变量类型与原模型一致（比如因子类型）
new_dat$sex <- factor(new_dat$sex, levels = levels(df$sex))
new_dat$location <- factor(new_dat$location, levels = levels(df$location))

2. 基于后验模拟计算拟合值 #

利用arm::sim()得到的后验参数对象（记为sim_obj），结合预测数据集的设计矩阵计算每个模拟下的拟合值：

情况1：Location为固定效应

如果模型中location是固定效应（例如公式为y ~ sex + size + location）：

# 生成预测数据集的设计矩阵
X <- model.matrix(formula(fit), new_dat)
# 计算1000组模拟下的拟合值：维度为1000个模拟 × 预测点数量
fit_vals <- t(X %*% t(sim_obj@fixef))

情况2：Location为随机效应

如果模型中location是随机截距/斜率（例如公式为y ~ sex + size + (1|location)）：

# 计算固定效应部分的预测值
X <- model.matrix(formula(fit), new_dat)
fixef_vals <- t(X %*% t(sim_obj@fixef))

# 提取每个模拟下的location随机效应值
ranef_sim <- sim_obj@ranef$location[, , 1]  # 维度：1000个模拟 × location的数量
# 匹配预测数据集中每个location对应的随机效应
ranef_vals <- ranef_sim[, match(new_dat$location, rownames(ranef_sim))]

# 总拟合值 = 固定效应预测值 + 随机效应值
fit_vals <- fixef_vals + ranef_vals

注：若要预测新的location（边际预测），则无需添加随机效应，直接使用固定效应部分的结果即可。

3. 计算均值与可信区间 #

对每个预测点的1000组模拟值进行汇总，得到拟合均值和95%可信区间：

new_dat$fit_mean <- apply(fit_vals, 2, mean)
new_dat$lci <- apply(fit_vals, 2, quantile, probs = 0.025)
new_dat$uci <- apply(fit_vals, 2, quantile, probs = 0.975)

4. 可视化结果 #

用ggplot2实现两种常见的可视化方案：

方案1：按Location分面板展示

适合location水平较多的场景，每个location单独成图便于观察：

library(ggplot2)
ggplot(new_dat, aes(x = size, y = fit_mean)) +
  geom_line(aes(color = sex), linewidth = 1) +
  geom_ribbon(aes(ymin = lci, ymax = uci, fill = sex), alpha = 0.2, color = NA) +
  facet_wrap(~location) +
  labs(x = "Size", y = "Fitted Value", color = "Sex", fill = "Sex") +
  theme_bw()

方案2：同一图中用线型区分Location

适合location水平较少的场景，便于直接对比不同location的结果：

ggplot(new_dat, aes(x = size, y = fit_mean)) +
  geom_line(aes(color = sex, linetype = location), linewidth = 1) +
  geom_ribbon(aes(ymin = lci, ymax = uci, fill = sex), alpha = 0.2, color = NA) +
  labs(x = "Size", y = "Fitted Value", color = "Sex", fill = "Sex", linetype = "Location") +
  theme_bw()

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