Hi all,

I have a model where there are different groups. I want to know whether these different groups are equal or not.

My question is based on part 12.2.2 and beyond in *Doing Bayesian Data Analysis* of J.K. Krushke.

Initially , I estimated the parameters of each group and compared the posterior differences between the parameter in question.

However, I would like to explore a model-comparison approach where I compare a model with a single shared parameter and a model where each group has his own distinct parameter.

Krushke says the following regarding the approach:

The model specification begins with saying that each subject has an individual ability

theta[s]from

a condition-specific beta distribution:

```
model {
for ( s in 1:nSubj ) {
nCorrOfSubj[s] ˜ dbin( theta[s] , nTrlOfSubj[s] )
theta[s] ˜ dbeta( aBeta[CondOfSubj[s]] , bBeta[CondOfSubj[s]] )
}
```

The shape parameters of the beta distribution are then re-written in terms of the mode and concentration. Model 1 uses condition-specific

omega[j], while Model 2 uses the sameomega0for all conditions. The JAGS functionequals(mdlIdx,…)is used to select the appropriate model for index mdlIdx:

```
for ( j in 1:nCond ) {
# Use omega[j] for model index 1, omega0 for model index 2:
aBeta[j] <- ( equals(mdlIdx,1)*omega[j] + equals(mdlIdx,2)*omega0 ) * (kappa[j]-2)+1
bBeta[j] <-(1-( equals(mdlIdx,1)*omega[j] + equals(mdlIdx,2)*omega0 ) ) * (kappa[j]-2)+1
omega[j] ˜ dbeta( a[j,mdlIdx] , b[j,mdlIdx] )
}
omega0 ˜ dbeta( a0[mdlIdx] , b0[mdlIdx] )
```

I would like to reproduce this in Stan, but I don’t know how I can achieve this.

Inspired by the JAGS example I write

```
data {
....
}
parameters {
vector<lower=0, upper=1>[n_trials] p;
....
}
model {
for (trial in 1:n_trials) {
int<lower=0, upper=1> model_idx;
model_idx ~ categorial(p[trial]);
for (group in 1:n_groups) {
if (model_idx == 0) {
theta[group] ~ beta(phi[group] * shared_kappa, (1 - phi[group] * shared_kappa);
y[group] ~ binomial(trials[group], theta[group]);
} else {
theta[group] ~ beta(phi[group] * kappa[group], (1 - phi[group] * kappa[group]);
y[group] ~ binomial(trials[group], theta[group]);
}
}
}
}
```

Is this the correct way? My goal is to reproduce the below chart.

Thank you for your time