BEST与brms包贝叶斯均值差异检验结果对比及技术疑问

针对你用BEST和brms做贝叶斯均值差异检验时遇到的三个问题，解答如下：

1. brms链接函数的作用、数据预处理及logm1链接说明 #

• 链接函数核心作用：建立模型线性预测器（即模型拟合的线性组合部分）与响应变量/分布参数（如sigma、nu）之间的映射，确保参数取值符合其固有约束（比如sigma必须为正、学生t分布的自由度nu必须大于1）。它不会直接修改原始数据，所有变换仅在模型内部计算流程中进行。
• logm1链接功能：公式为 g(x) = log(x - 1)，作用是将参数的取值范围从 (1, +∞) 映射到整个实数轴（(-∞, +∞)），让线性预测器可以取任意实数值，避免拟合时出现参数违反分布约束的情况。
• 为何是nu的默认选项：学生t分布的自由度nu必须大于1（否则分布无意义），logm1链接恰好完美适配这个参数的取值范围，是最合理的映射方式，因此被设为默认。

2. 结合链接函数设置先验分布的规则 #

先验分布的参数是基于链接函数变换后的线性预测器尺度设置的，具体分两种情况：

• 若使用identity链接（如你代码中响应变量y的链接）：线性预测器直接等于原始参数（比如mu），此时先验可直接基于原始数据尺度设置（比如你给groupy1的mu设normal(6, 2)，和BEST中的先验完全匹配）。
• 若使用其他链接（如sigma的log链接、nu的logm1链接）：先验是针对链接变换后的线性预测器设置的。例如你代码中sigma的链接是log，那么set_prior("cauchy(1.015332, ...)", dpar="sigma")里的1.015332，其实是原始sigma先验位置的对数（log(2.76)≈1.015），对应原始尺度下sigma的先验位置为2.76。

简言之：先验参数要对应「链接函数作用在原始参数上的结果」的分布，而非原始参数本身的分布（除非是identity链接）。

3. brms结果逆变换回原始尺度的方法 #

brms默认输出的参数值已经是原始尺度的，但如果需要从线性预测器（即链接变换后的数值）手动转换，对应规则如下：

• identity链接：无需变换，直接使用线性预测器的值（比如你代码中的b_groupy1就是原始尺度的mu）。
• log链接：用exp()逆变换，比如你代码中sigma1 <- mean(exp(brmout$b_sigma_groupy1))，就是把log变换后的线性预测器转换回原始sigma尺度。
• logm1链接：逆变换公式是exp(x) + 1。因为logm1链接是g(nu) = log(nu - 1)，解出原始nu就是nu = exp(g(nu)) + 1。注意：brms输出的nu列已经是逆变换后的原始尺度值，所以你代码中直接mean(brmout$nu)是正确的，无需额外处理。

附整理后的R代码 #

library(BEST)

y1 <- c(5.77, 5.33, 4.59, 4.33, 3.66, 4.48)
y2 <- c(3.88, 3.55, 3.29, 2.59, 2.33, 3.59)

best_priors <- list(muM = c(6, 4), muSD = 2)
BESTout <- BESTmcmc(y1, y2, priors=best_priors, parallel=FALSE, rnd.seed = 123)
print(BESTout)

library(brms)

y <- data.frame(y = c(y1, y2), group = c(rep("y1", 6), rep("y2", 6)))
with(y, mean(y))
with(y, sd(y))

brmout <- brm(
    formula = brmsformula(y ~ 0 + group, sigma ~ 0 + group),
    data    = y,
    family  = student(link = "identity", link_sigma = "log", link_nu = "logm1"),
    prior   = c(
        set_prior("normal(6, 2)",                 class = "b", coef = "groupy1"),
        set_prior("normal(4, 2)",                 class = "b", coef = "groupy2"),
        set_prior("cauchy(1.015332, 1.015332*5)", class = "b", coef = "groupy1", dpar = "sigma"),
        set_prior("cauchy(1.015332, 1.015332*5)", class = "b", coef = "groupy2", dpar = "sigma"),
        set_prior("gamma(30, 30)",                class = "nu")
    ),
    seed    = 123
)

brmout
brmout_df <- as.data.frame(brmout)

summary(BESTout)

(mu1 <- mean(brmout_df$b_groupy1))
(mu2 <- mean(brmout_df$b_groupy2))
(muDiff <- mean(brmout_df$b_groupy1 - brmout_df$b_groupy2))

(sigma1 <- mean(exp(brmout_df$b_sigma_groupy1)))
(sigma2 <- mean(exp(brmout_df$b_sigma_groupy2)))
(sigmaDiff <- sigma1 - sigma2)

(nu <- mean(brmout_df$nu))

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