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 same omega0 for all conditions. The JAGS function equals(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