Dear all

I am new to Stan, and as a learning exercise I wanted to estimate a Bayesian multilevel linear model for the `sleepstudy`

dataset (in the `lme4`

package). The model have subject-specific intercepts and slopes, assumed to have a normal distribution across subjects, with (possibly) non-zero correlation. The model seems to run correctly, and I get estimates that are sufficiently close to what I obtain with `lme4`

.

I wanted to simulate the model to do a posterior predictive check (estimating p \left( y^{\text{rep}} \mid y \right)). In order to do that I am simulating new observations, from new simulated subjects, in the `generated quantities`

block. Is that the best way to proceed? It seems to work correctly (see plot; the blue bands represents the 95% percentile intervals of the distribution of simulated y^{\text{rep}})

Moreover, the only way I found to generate the random intercepts and slopes is with a loop. Is that the only way, or there is a more efficient way to generate directly a matrix of normal random variates with given variance-covariance matrix?

Thanks!

```
data {
int<lower=1> N; // number of observations
real RT[N]; // reaction time (transformed in seconds)
int<lower=0,upper=9> Days[N]; // predictor
int<lower=1> J; // number of subjects
int<lower=1,upper=J> Subject[N]; // subject id
}
parameters {
vector[2] beta; // fixed-effects parameters
real<lower=0> sigma_e; // residual std
vector<lower=0>[2] sigma_u; // random effects standard deviations
cholesky_factor_corr[2] L_u; // L_u is the Choleski factor of the correlation matrix
matrix[2,J] z_u; // random effect matrix
}
transformed parameters {
matrix[2,J] u;
u = diag_pre_multiply(sigma_u, L_u) * z_u; // use Cholesky to set correlation
}
model {
real mu; // conditional mean of the dependent variable
//priors
L_u ~ lkj_corr_cholesky(2); // LKJ prior for the correlation matrix
to_vector(z_u) ~ normal(0,1); // before Cholesky, random effects are normal variates with SD=1
sigma_u ~ cauchy(0, 0.5); // SD of random effects (vectorized)
sigma_e ~ cauchy(0, 0.5); // prior for residual standard deviation
beta[1] ~ normal(0.3, 1); // prior for fixed-effect intercept
beta[2] ~ normal(0.01, 0.5); // prior for fixed-effect slope
//likelihood
for (i in 1:N){
mu = beta[1] + u[1,Subject[i]] + (beta[2] + u[2,Subject[i]])*Days[i];
RT[i] ~ normal(mu, sigma_e);
}
}
generated quantities {
// simulate model for posterior predictive check
real y_rep[N];
matrix[2,J] u_rep; // simulate 'random effects'
for (j in 1:J) {
u_rep[1:2,j] = multi_normal_rng(rep_vector(0,2), diag_matrix(sigma_u) * (L_u*L_u') * diag_matrix(sigma_u));
}
for (i in 1:N){
y_rep[i] = normal_rng(beta[1] + u_rep[1,Subject[i]] + (beta[2] + u_rep[2,Subject[i]])*Days[i], sigma_e);
}
}
```

ppc_sleepstudy.pdf (15.0 KB)