Sure, I was trying to keep it simple, but here is the whole thing. If I figure out this multiple indexing thing, I can also make other improvements to the code.

Part of R run script:

```
stan_data <- list(
N = N,
I = I,
id = id,
group = group,
w = w,
x = x,
ECT = ECT,
MMSE = MMSE,
indices1 = as.integer(c(1,3,4,5,7,8)),
indices2 = as.integer(c(2,6,9,10,11,12))
)
model <- cmdstan_model(paste0(getwd(), '/', model_name,'.stan'))
fit <- model$sample(
data = stan_data,
seed = as.numeric(Sys.Date()),
chains = 8,
iter_warmup = 2000,
iter_sampling = 10000,
refresh = 100,
thin = 10,
adapt_delta = .95
)
```

Stan model:

```
data {
int<lower=1> I; // number of subjects
int<lower=1> N; // total observations; each n is ij
int<lower=1> id[N]; // id of nth observation
int<lower=1> group[I]; // diagnostic group of subject id i
vector[2] ECT[N]; //
int<lower=0, upper=30> MMSE[N];
matrix[N, 2] x; // ones column, t column (t is years since initial visit)
matrix[N, 4] w; // age_c, male, age_c_t, male_t
int indices1[6];
int indices2[6];
}
transformed data {}
parameters {
vector[12] alpha[3]; // fixed effects
real<lower=0> phi[3]; // scale parameter
vector[6] beta[I]; // random effects; (i, m) = (subject, index)
vector[6] mu[3]; // means of random effects; (l, m) = (group, index)
vector<lower=0>[6] tau[3]; // sd's of random effects; (l, m) = (group, index)
cholesky_factor_corr[6] L_Omega[3]; // random effects correlation matrix (l) = group
vector<lower=0>[2] sigma[3]; // error sd's for ECT
cholesky_factor_corr[2] L_Rho[3]; // error correlations for ECT (l) = group
}
transformed parameters {}
model {
// LIKELIHOOD
for(n in 1:N){
ECT[n] ~ multi_normal(
[ row(w, n) * alpha[group[id[n]]][1:4] + row(x, n) * beta[id[n]][1:2],
row(w, n) * alpha[group[id[n]]][5:8] + row(x, n) * beta[id[n]][3:4] ],
diag_pre_multiply(sigma[group[id[n]]], L_Rho[group[id[n]]]) *
diag_pre_multiply(sigma[group[id[n]]], L_Rho[group[id[n]]])'
);
MMSE[n] ~ neg_binomial_2(
exp(row(w, n) * alpha[group[id[n]]] + row(x, n) * beta[id[n]]),
phi[group[id[n]]]
);
}
for(i in 1:I){
beta[i] ~ multi_normal(
mu[group[i]],
diag_pre_multiply(tau[group[i]], L_Omega[group[i]]) *
diag_pre_multiply(tau[group[i]], L_Omega[group[i]])'
);
}
// PRIORS
for(l in 1:3){
L_Rho[l] ~ lkj_corr_cholesky(1);
L_Omega[l] ~ lkj_corr_cholesky(.8);
alpha[l][indices1] ~ normal(0,.25);
alpha[l][indices2] ~ normal(1, 2);
}
// FSLong intercepts
mu[1][1] ~ normal(6, 3);
mu[2][1] ~ normal(6, 3);
mu[3][1] ~ normal(6, 3);
tau[1][1] ~ exponential(2);
tau[2][1] ~ exponential(2);
tau[3][1] ~ exponential(2);
// FSLong slopes
mu[1][2] ~ normal(0, .5);
mu[2][2] ~ normal(0, .5);
mu[3][2] ~ normal(0, .5);
tau[1][2] ~ exponential(4);
tau[2][2] ~ exponential(4);
tau[3][2] ~ exponential(4);
// ANTs intercepts
mu[1][3] ~ normal(7, 3);
mu[2][3] ~ normal(7, 3);
mu[3][3] ~ normal(7, 3);
tau[1][3] ~ exponential(1);
tau[2][3] ~ exponential(1);
tau[3][3] ~ exponential(1);
// ANTs slopes
mu[1][4] ~ normal(0, .5);
mu[2][4] ~ normal(0, .5);
mu[3][4] ~ normal(0, .5);
tau[1][4] ~ exponential(10);
tau[2][4] ~ exponential(10);
tau[3][4] ~ exponential(10);
// MMSE intercepts
mu[1][5] ~ normal(0, 1);
mu[2][5] ~ normal(1, 1);
mu[3][5] ~ normal(1, 1);
tau[1][5] ~ exponential(2);
tau[2][5] ~ exponential(2);
tau[3][5] ~ exponential(2);
// MMSE slopes
mu[1][6] ~ normal(0, 1);
mu[2][6] ~ normal(0, 1);
mu[3][6] ~ normal(0, 1);
tau[1][6] ~ gamma(5, 50);
tau[2][6] ~ exponential(5);
tau[3][6] ~ exponential(5);
}
generated quantities {}
```