Our data are longitudinal, with multiple observations per person i over time t.

First, we fit a standard multilevel logistic model:

E[y_{it} | x_{it}, \eta_i] = \text{logit}^{-1}(\beta_0 + \beta_1 x_{1,it} + \beta_2 x_{2,it} + \beta_3 x_{3,it} + \eta_i)

\eta_i \sim N(0,\sigma_\eta)

This model appeared to converge (no obvious fitting issues). But our posterior predictive check failed. We got advice to try a mixture model for \eta_i:

z_i \sim \text{Categorical}(\theta_1,\theta_2)

\eta_i | z_i \sim N(\mu_{z_i},\sigma_{z_i})

where we enforce:

E[\eta_i] = \theta_1\mu_1 + \theta_2\mu_2 = 0

However, the model (as coded below) does not converge (we used 4 chains, 2000 iterations). We cannot share the data, but will try simulated data next. Any advice just based on the code? Errors?

Thanks so much!

```
data {
int<lower=1> num_obs;
int<lower=1> num_people;
int y[num_obs];
int<lower=1,upper=num_people> ii[num_obs]; // person
real x1[num_obs];
real x2[num_obs];
real x3[num_obs];
}
parameters {
real person_effect[num_people];
simplex[2] theta; // mixing proportions
real mu_subset; // first location of mixture components
vector<lower=0>[2] sigma; // scales of mixture components
real intercept;
real beta_1;
real beta_2;
real beta_3;
}
model {
vector[2] log_theta = log(theta); // cache log calculation
vector[2] mu; // locations of mixture components
real y_mean_logit[num_obs];
mu[1] = mu_subset;
mu[2] = (- theta[1] * mu[1])/theta[2];
intercept~normal(0,4);
beta_1~normal(0,4);
beta_2~normal(0,4);
beta_3~normal(0,4);
mu~normal(0,4);
sigma~normal(0,4);
for (i in 1:num_people) {
vector[2] lps = log_theta;
for (k in 1:2)
lps[k] += normal_lpdf(person_effect[i] | mu[k], sigma[k]);
target += log_sum_exp(lps);
}
for(n in 1:num_obs)
y_mean_logit[n] = intercept+beta_1*x1[n]+beta_2*x2[n]+beta_3*x3[n]+person_effect[ii[n]];
y ~ bernoulli_logit(y_mean_logit);
}
```