The R packages "BEST" and "brms" can both calculate a Bayesian difference in means test. I'm trying to compare the results from both packages. My questions are: (1) what do the link functions do in the brms package? Do they preprocess the data somehow? What does the "logm1" link function do and why is that the default?, (2) In light of the brms 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) How do I back transform the brms results into the original scale? For example, I think the exp() function back transforms data from the "log" link function (maybe not though when taking means), but how do I back transform data from the "logm1" link function? 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)))
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 <- 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))
