如何在R的brms包中指定随机效应间相关性并复现MCMCglmm模型？

架构老林

2026-5-19

当然可以用brms复现MCMCglmm里的多变量多层模型！结果差异大大概率是因为两边的默认设置（比如先验、MCMC参数）没对齐，咱们一步步来搞定：

1. 先锚定MCMCglmm的模型结构

假设原文的核心模型是多因变量的多层随机截距模型，重点估计不同因变量随机截距间的无约束协方差，代码大概长这样（如果你的原文代码有差异，直接替换对应部分就行）：

library(MCMCglmm)
# 加载数据
dat <- read.csv("appendix_example1_wide.SUPP.FINAL.csv")

# MCMCglmm多变量模型
model_mcmc <- MCMCglmm(
  # 固定效应：每个因变量的独立截距 + 固定协变量（假设x_fixed是固定效应）
  cbind(y1, y2, y3) ~ trait - 1 + trait:x_fixed,
  # 随机效应：id分组下，每个因变量的随机截距，无约束协方差结构（us=unstructured）
  random = ~ us(trait):id,
  # 残差协方差：无约束结构（不同因变量的残差允许相关）
  rcov = ~ us(trait):units,
  data = dat,
  # 原文指定的先验（这部分一定要对齐！）
  prior = list(
    R = list(V = diag(3), nu = 4),  # 残差协方差先验：逆Wishart，自由度4，尺度矩阵单位阵
    G = list(G1 = list(V = diag(3), nu = 4))  # 随机效应协方差先验：同残差设置
  ),
  # MCMC采样参数
  nitt = 10000, burnin = 2000, thin = 8
)

2. brms中对应的模型写法

brms处理多变量模型的逻辑和MCMCglmm略有不同，但完全可以做到1:1对齐，关键是协方差结构、先验、MCMC参数三者要严格匹配：

library(brms)
# 加载同一数据集
dat <- read.csv("appendix_example1_wide.SUPP.FINAL.csv")

# 构建多变量模型公式：每个因变量的固定效应和MCMCglmm一致
model_formula <- mvbind(y1, y2, y3) ~ 0 + trait + trait:x_fixed +
  # 随机效应：id分组下的多变量随机截距，无约束协方差（对应MCMCglmm的us）
  (0 + trait | id)

# 拟合模型
model_brms <- brm(
  formula = model_formula,
  data = dat,
  family = gaussian(),
  # 对齐原文的先验：完全匹配MCMCglmm的逆Wishart先验
  prior = c(
    # 随机效应协方差矩阵：逆Wishart先验，自由度4，尺度矩阵单位阵
    prior(inv_wishart(4, diag(3)), class = cov, coef = trait, group = id),
    # 残差协方差矩阵：开启残差相关（对应MCMCglmm的rcov=us），并设置相同先验
    prior(inv_wishart(4, diag(3)), class = cov, resp = y1),
    prior(inv_wishart(4, diag(3)), class = cov, resp = y2),
    prior(inv_wishart(4, diag(3)), class = cov, resp = y3)
  ),
  rescor = TRUE,  # 必须开这个！否则brms默认残差独立，和MCMCglmm的rcov设置矛盾
  # 对齐MCMC采样参数：和MCMCglmm的nitt/burnin/thin完全一致
  iter = 10000, warmup = 2000, thin = 8,
  chains = 1,  # MCMCglmm默认1链，brms默认3链，这里统一成1链方便对比
  control = list(adapt_delta = 0.95)  # 提升数值稳定性，避免采样发散
)

3. 结果差异大的常见排查点

如果还是差很多，优先检查这几个地方：

先验不一致：这是最常见的坑！MCMCglmm的默认先验非常弱（比如G矩阵默认nu=0.002），而brms默认的先验更具信息性，必须严格照搬原文的先验设置。
随机效应结构写错：比如你是不是写成了(1 | id)？这是单变量随机截距，和MCMCglmm的us(trait):id完全不是一回事，必须用(0 + trait | id)来指定每个因变量的独立随机截距。
残差协方差设置：brms默认不估计残差间的协方差，一定要加rescor = TRUE，否则和MCMCglmm的rcov=~us(trait):units差异极大。
MCMC收敛性：检查MCMCglmm的自相关系数（autocorr(model_mcmc)）和brms的Rhat值（rhat(model_brms)），确保两者都收敛良好（Rhat<1.01，自相关系数快速衰减）。

4. 验证一致性的小技巧

先对比固定效应估计值：如果固定效应都差很多，说明公式或者数据处理有问题；固定效应接近后，再看随机效应协方差。
提取协方差矩阵：MCMCglmm用model_mcmc$VCV查看随机效应协方差，brms用VarCorr(model_brms)或者posterior_samples(model_brms, "^cov_")提取对应后验样本。

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