# 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.

4 Likes

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