R BEST package how to set Bayesian priors by group

I'm using the R BEST package to test for a difference in means between two groups, and I'd like to set the prior beliefs by group.

In the R code below, I am able to set the prior means by group (see priors2), but not the prior standard deviations by group (see priors3).
Am I doing something wrong?

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)

priors1 <- list(muM = 6, muSD = 2)
BESTout1 <- BESTmcmc(y1, y2, priors=priors1, parallel=FALSE)

priors2 <- list(muM = c(6, 4), muSD = 2)
BESTout2 <- BESTmcmc(y1, y2, priors=priors2, parallel=FALSE)

priors3 <- list(muM = c(6, 4), muSD = c(2, 2))
BESTout3 <- BESTmcmc(y1, y2, priors=priors3, parallel=FALSE)

The R BEST package is a Bayesian version of a two-sample t-test (John K. Kruschke).

From the example code you have provided, it seems like you are trying to use different priors (ways to incorporate existing knowledge or beliefs about a parameter before the data is considered) for the two groups you are comparing.

However, the way the Bayesian t-test is implemented in the BEST package is that it uses a single shared prior for the standard deviation of the groups. This is because it assumes that the two groups being compared are likely to have similar variances.

In the code you have posted, the muM argument refers to the prior mean of the group means, and the muSD argument refers to the prior standard deviation of the group means. While you can set different prior means for each group (as you have done in priors2), the muSD parameter is shared between the groups.

Unfortunately, the BEST package does not support setting different prior standard deviations for each group.
You can see the details in mikemeredith/BEST issue 6

If you are interested in setting different priors for the standard deviation of each group, you will need to use a different package or perform the analysis manually.

One alternative could be to use the paul-buerkner/brms package (presented in this paper), which is a versatile package for Bayesian analysis in R. It uses Stan (a platform for Bayesian statistics) for its computations, and supports a wide variety of models and priors. However, note that using brms or Stan directly can be more complex than using the BEST package, as they require you to specify the model and priors in more detail.

The Bayesian t-test in brms could look like this:

library(brms)

# create a data frame
dat <- data.frame(
  y = c(y1, y2),
  group = factor(rep(c("y1", "y2"), each = 6))
)

# fit the model
fit <- brm(
  y ~ group,
  prior = c(
    set_prior("normal(6, 2)", class = "Intercept", coef = "groupy1"),
    set_prior("normal(4, 2)", class = "Intercept", coef = "groupy2")
  ),
  data = dat,
  family = gaussian()
)

This code creates a Bayesian linear regression model where the means of the two groups are modeled as intercepts, and different priors are set for each group.

This is just an example, and you will need to choose the appropriate priors for your specific case.

You can set priors in brms. For example:

set_prior("normal(0,5)", class = "b", coef = "x1") 
set_prior("student_t(10, 0, 1)", class = "b", coef = "x2")

Here, we are setting a normal prior with mean 0 and standard deviation 5 for the coefficient x1, and a Student's t prior with 10 degrees of freedom for x2. This code would need to be used in conjunction with a brm function call to run a model.

One key thing to note is that the coef argument in set_prior corresponds to the variable names in your data, not to the group names. If you have group-specific effects that you want to model, you might need to consider using group-level or random effects in your model.

To integrate this with the model code, you would do something like:

library(brms)
prior1 <- set_prior("normal(0,5)", class = "b", coef = "x1") 
prior2 <- set_prior("student_t(10, 0, 1)", class = "b", coef = "x2")
priors <- c(prior1, prior2)
mod_eqvar <- brm(
  IQ ~ Group,
  data = d,
  cores = 4,  # Use 4 cores for parallel processing
  prior = priors
)
summary(mod_eqvar)

(Do replace x1 and x2 with the actual names of the variables in your dataset for which you wish to set priors.)

This code sets the specified priors for x1 and x2 and uses them in the model.
It then performs a Bayesian t-test comparing IQ between different Groups in the dataframe d, using the specified priors.
The results of the model are then viewed with summary(mod_eqvar).

To compare groups, the initial code above would become:

# Load the library
library(brms)

# Set a seed for reproducibility
set.seed(42)

# Generate some data
group1 <- rnorm(50, mean = 5, sd = 1)
group2 <- rnorm(50, mean = 6, sd = 1)

# Combine into a data frame
d <- data.frame(
  IQ = c(group1, group2),
  Group = factor(rep(c("Group1", "Group2"), each = 50))
)

# Set priors for the model. This prior represents a prior belief about the difference in means between Group2 and Group1
priors <- set_prior("normal(0, 5)", class = "b", coef = "Group2")

# Fit the model
mod_eqvar <- brm(
  IQ ~ Group,
  data = d,
  prior = priors,
  cores = 4,  # Use 4 cores for parallel processing
  file = "iqgroup"  # Save results into a file
)

# View the results
summary(mod_eqvar)

In this code, the model IQ ~ Group is comparing the variable IQ between different Groups in the dataframe d. In the model formula, IQ is the outcome variable and Group is the predictor variable. The predictor variable Group is categorical and it represents different groups in your data.

(See also "How to Compare Two Groups with Robust Bayesian Estimation in R")

In the model formula IQ ~ Group, IQ is the outcome variable and Group is the predictor variable. The predictor variable Group is categorical, and it represents different groups in your data.

When using factor variables in a regression model in R, the first level of the factor is taken as the reference group, and the estimates for the other levels represent the difference in the outcome variable (in this case, IQ) relative to the reference group. The estimate for "Group2" in the model represents the difference in means between Group2 and the reference group (Group1).

So, if you set a prior on the coefficient for "Group2" in the model, like in this line of code:

priors <- set_prior("normal(0, 5)", class = "b", coef = "Group2")

You are setting a prior belief about the difference in means between Group2 and the reference group, not about the overall mean for Group2. This prior distribution is then updated based on the data when you fit the model, resulting in a posterior distribution for the difference in means between the two groups.

The summary(mod_eqvar) function will give you a summary of the model, including the estimated difference in means between Group2 and Group1. This difference is represented by the "Group2" estimate in the output. In a Bayesian context, we typically look at the posterior distribution of the parameter of interest (in this case, the difference in means between the groups) and the 95% credible interval. If the 95% credible interval for the difference in means does not include 0, it suggests that there is a statistically significant difference in means between the groups.

The set_prior function is used to set a prior distribution for the difference in means between the groups in the model. The coef argument in set_prior corresponds to the variable names in your data, not to the group names. So when we use coef = "Group2" in set_prior, we are actually setting a prior for the difference in means between Group2 and the reference group (Group1), not for the overall mean of Group2.

(See also "Prior Definitions for brms Models")

In a Bayesian context, instead of p-values, we typically look at the posterior distribution of the parameter of interest (in this case, the effect of Group on IQ) and the 95% credible interval. If the credible interval does not include zero, it suggests that the effect is statistically significant at the 0.05 level.

About the error:

Error: The following priors do not correspond to any model parameter: 
  Intercept_groupy1 ~ normal(6, 2) 
  Intercept_groupy2 ~ normal(4, 2) 
Function 'get_prior' might be helpful to you. 

The error message indicates that the priors specified do not match any of the model parameters. This is likely because the group names do not match those in the data.

In the brms model syntax, the intercept term represents the overall mean for the reference group, and the coef terms represent the differences from this reference group mean for the other groups. The coef term should match the exact name of the group levels in your factor variable, which appears to be Group1 and Group2 in your case.

Let's correct the coef terms in the set_prior function to match the group names in your data.
You can try and replace Intercept_groupy1 and Intercept_groupy2 with Intercept_Group1 and Intercept_Group2 (or the exact names of the levels of your factor variable) in the set_prior function, like this:

# Set priors for the model
priors <- c(
  set_prior("normal(6, 2)", class = "Intercept", coef = "Group1"),
  set_prior("normal(4, 2)", class = "Intercept", coef = "Group2")
)

# Fit the model
mod_eqvar <- brm(
  IQ ~ Group,
  data = d,
  prior = priors,
  cores = 4,  # Use 4 cores for parallel processing
  file = "iqgroup"  # Save results into a file
)

Note: It is possible that this still will not work, because brms may not allow different priors for different levels of a factor variable within the same model. In this case, you may need to fit separate models for each group and then compare the results manually.