The R package “BEST” calculates a simple Bayesian difference in means test. I’m trying to replicate the BEST package results using the brms package because the brms package has more flexibility and I want to learn how it works. My questions are: (1) what do the family link functions do in the brms package? Do they preprocess the data somehow?, (2) In light of the brms family link functions, how do I specify the parameters of the brms prior distributions? Do I specify the parameters as if the data were preprocessed or do I assume no preprocessing?, (3) The built-in priors for the BEST package for “sigma” and “nu” follow the Gamma distribution, which is what I’m trying to replicate. How do I correctly express those priors in the brms code? Thank you. Below is my R code.
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)))
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("gamma(sd(y1), sd(y1)*5)", class = "b", coef = "groupy1", dpar = "sigma"),
set_prior("gamma(sd(y2), sd(y2)*5)", class = "b", coef = "groupy2", dpar = "sigma"),
set_prior("gamma(30, 30)", class = "nu")
),
seed = 123
)
brmout
brmout <- as.data.frame(brmout)
summary(BESTout)
(mu1 <- mean(brmout$b_groupy1))
(mu2 <- mean(brmout$b_groupy2))
(muDiff <- mean(brmout$b_groupy1 - brmout$b_groupy2))
(sigma1 <- mean(exp(brmout$b_sigma_groupy1)))
(sigma2 <- mean(exp(brmout$b_sigma_groupy2)))
(sigmaDiff <- sigma1 - sigma2)
(nu <- mean(brmout$nu))