I’m writing my first models in Stan and today I started with a multivariate model with 2 outcomes. To capture heterogeneity across individuals (my obervations are grouped in multiple individuals), next to the intercept alpha0 for each equation, I want to add an individual-specific intercept alpha to capture heterogeneity across individuals.

Here is my model code:

```
data {
int <lower = 1> N; // Number of observations
int<lower = 1> J; // number of outcomes
int<lower = 0> K; // number of covariates
vector[J] Y[N]; // vector of J outcome variables
vector[K] X[N]; // vector of K covariates
int <lower = 1> id[N]; // identifier of individuals
int <lower = 1> H; // number of individuals
}
parameters {
vector[J] alpha0; // intercept for each equation
vector[J] alpha[H]; // individual-specific intercept for each equation
vector<lower = 0>[J] sigma_alpha;
cholesky_factor_corr[J] L_Omega;
vector<lower=0>[J] L_sigma;
matrix[J, K] beta;
}
model {
vector[J] mu[N];
alpha[J] ~ normal(0, sigma_alpha);
for (n in 1:N)
mu[n] = alpha0 + alpha[id[n]] + beta * X[n];
to_vector(beta) ~ normal(0, 2);
L_Omega ~ lkj_corr_cholesky(1);
L_sigma ~ student_t(3, 0, 2);
Y ~ multi_normal_cholesky(mu, diag_pre_multiply(L_sigma, L_Omega));
}
```

Now I realized my data requires me to apply the non-cenetred parameterization. I’m struggling with the parameters now because I’m not confident which paramaters need to be individual-specific then.

Here is my code for a multivariate model capturing heterogeneity across individuals with a non-centered parameterization for the individual-specific intercept alpha.

```
data {
int <lower = 1> N; // Number of observations
int<lower = 0> J; // number of outcomes
int<lower = 0> K; // number of covariates
vector[J] Y[N]; // vector of J outcome variables
vector[K] X[N]; // vector of K covariates
int <lower = 1> id[N]; // identifier of individuals
int <lower = 1> H; // number of individuals
}
parameters {
real alpha0;
vector[H] alphai_raw;
real<lower = 0> sigma_alphai;
cholesky_factor_corr[J] L_Omega;
vector<lower=0>[J] L_sigma;
matrix[J, K] beta;
}
transformed parameters {
vector[H] alphai;
alphai = alpha0 + sigma_alphai*alphai_raw;
}
model {
vector[J] mu[N];
alphai_raw ~ cauchy(0,.2);
alpha0 ~ cauchy(0, 1);
for (n in 1:N)
mu[n] = alphai[id[n]] + beta * X[n];
to_vector(beta) ~ normal(0, 2);
L_Omega ~ lkj_corr_cholesky(1);
L_sigma ~ student_t(3, 0, 2);
Y ~ multi_normal_cholesky(mu, diag_pre_multiply(L_sigma, L_Omega));
}
```

- I’m wondering if it’s right that I can get rid of the non-individual-specific intercept alpha0 when defining mu[n] in the model section?
- The way I have coded it I get one intercept for each individual for the two equations together, but how can I get one intercept for each individual for each equation?