Pairwise Mean Comparison BRMS

Interfaces brms
1

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.

2

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.

4 Likes
3

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

4

Thanks Peter for the very clear example and Matti for the very helpful answer (and blog post)!

I have a follow-up question about whether it’s possible to do the same hypothesis tests in a 2x2 model with sum contrast coded predictors? For example, in the dummy data below, I would like to quantify whether [cond1 -0.5, cond2 -0.5] is different to [cond1 -0.5, cond2 0.5]… is that possible?

I tried fitting a separate model with four factorised conditions similar to your example above which was fine, but then it’s a slightly different model to the one I’m basing inference on so it would be nice to be able to probe the model below.

test_data <- tibble(
  cond1 = rep(c(-.5,-.5,.5,.5), 25),
  cond2 = rep(c(-.5,.5,-.5,.5), 25),
  X = rnorm(100, 10, 5)
)

priors <- c(
  prior(normal(0, 5), class = Intercept),
  prior(normal(0, 1), class = b), 
  prior(normal(8, 2), class = sigma)
)

m <- brm(X ~ cond1 * cond2,
         prior = priors,
         family = gaussian(),
         data   = test_data)

Many thanks!
Kate