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.