Pairwise Mean Comparison BRMS

Hi all,

I’m fairly new to Stan and BRMS, so maybe this will sound like a stupid question, but here goes…
I’m wondering how to model pairwise mean comparison between groups in brms.

For example, I have data on some variable X and the group of an observation (code in r):

require(rstan)
require(brms)
require(tidyverse)


test_data <- tibble(
  Group = rep(1:3, each = 30),
  X = c(rnorm(30, 10), rnorm(30, 20), rnorm(30, 30))
)

If I fit a simple brms model, like so:

fit_brm <- brm(X ~ factor(Group),
               data = test_data)
print(fit_brm)

This will give me group effects on the dependent variable.

I am, however, more interested in the comparing the means directly, similarly to a t-test, or in multivariate case, to an ANOVA.

something like this:

aov_fit <- aov(lm(X ~ factor(Group),
              data = test_data))
summary(aov_fit)
TukeyHSD(aov_fit)

I’ve also written (a very crude) Stan program to estimate means and sds from data, and was wondering whether it is possible to replicate it in brms:

data{
  int<lower=0> N;
  vector[N] x;
  int<lower=1> J;
  int<lower=1, upper=J> Subject[N];
}
parameters{
  real<lower=0> nu;
  vector<lower=0>[J] sigma;
  vector[J] mu;
}
model{
  for (i in 1:N) {
    target += normal_lpdf(nu | 30, 30);

    target += cauchy_lpdf(sigma[Subject[i]] | 1, 20);
    target += normal_lpdf(mu[Subject[i]] | 0, 30);

    target += student_t_lpdf(x[i] | nu, mu[Subject[i]], sigma[Subject[i]]);
  }
} 


data_list <- list(
  N = nrow(test_data),
  x = test_data$X,
  Subject = test_data$Group,
  J = length(unique(test_data$Group))
)


fit_stan <- sampling(my_model,
                     data = data_list)

Given that now the output is the distribution of parameters, i can directly subtract them and figure out how much of an overlap there is.

Hi @Peter, welcome to the Stan forums,

Your model has three important parameters, the mean of the reference group (Intercept), and differences from that to groups 2 and 3.

Then, similar to t-tests, the parameters Group2 and Group3 are differences to Group1, with standard errors and CIs. To get the group 2 vs group 3 “t-test” (often called a “contrast”), you can do:

hypothesis(fit_brm, "Intercept + Group2 = Intercept + Group3")

Which returns:

Hypothesis Tests for class b:
                Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
1 (Intercept+Group2... = 0    -9.93      0.25   -10.42    -9.44         NA        NA    *
---
'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
'*': For one-sided hypotheses, the posterior probability exceeds 95%;
for two-sided hypotheses, the value tested against lies outside the 95%-CI.
Posterior probabilities of point hypotheses assume equal prior probabilities.

See ?hypothesis, and my recent blog post about how to use the hypothesis() function.

As an alternative, you can estimate the model without an intercept:

fit_brm2 <- brm(
  X ~ 0 + Group,
  data = test_data,
  cores = 4
)

That model’s parameters will be the group means. Then you can use hypothesis() to calculate differences in them, e.g.:

hypothesis(fit_brm2, "Group1 = Group2")

You will see that the result is the same as the Group2 parameter in the first model.

3 Likes

Perfect, thanks! The hypothesis function is a useful tool to know!