Multilevel logistic regression in STAN

multilevel_logistic_stan

// random effects covariances from
// https://github.com/vasishth/StanJAGSexamples/blob/master/FrankEtAlCogSci2015/FTVCogSciProficiency.Stan

data {
    int<lower=1> N;                         // Number of observations
    int<lower=1> J;                         // Number of subjects
    int<lower=1> K;                         // Number of predictors
    int<lower=1,upper=J> id[N];             // Subject ids
    real xInt[N];                           // ISI vector
    real<lower=-0.5,upper=0.5> xCond[N];    // Condition vector
    real xIntByCond[N];                     // Interaction vector
    int<lower=0, upper=1> y[N];             // Outcomes specified as 0s and 1s
}
parameters {
    vector[K] beta;                         // fixed effects
    vector<lower=0,upper=10>[K] sigma_u;    // SDs for raneffs
    
    // Correlation matrix for random intercepts and slopes
    cholesky_factor_corr[K] L_u; 
    
    matrix[K,J] z_u;                        // Random effects
    
}
transformed parameters{
    // These are only useful in transforming random effects to "random coefficients"
    matrix[J,K] u;                          // Random effects
    vector[J] u_int;
    vector[J] u_isi;
    vector[J] u_con;
    vector[J] u_iXc;
    u <- (diag_pre_multiply(sigma_u,L_u)*z_u)';
    for (i in 1:J){
        u_int[i] <- u[i,1] + beta[1];
        u_isi[i] <- u[i,2] + beta[2];
        u_con[i] <- u[i,3] + beta[3];
        u_iXc[i] <- u[i,4] + beta[4];
    }
}
model{
    real mu[N];                             // Mean of likelihood
    beta[1] ~ normal(0, 10);                // Fix intercept prior
    beta[2] ~ normal(0, 10);                // Fix ISI prior
    beta[3] ~ normal(0, 10);                // Fix Condition prior
    beta[4] ~ normal(0, 10);                // Fix interaction prior
    L_u ~ lkj_corr_cholesky(2.0);           // Ranef prior
    to_vector(z_u) ~ normal(0,1);
    
    for (n in 1:N)
        mu[n] <- (beta[1] + u[id[n],1]) +
                 (beta[2] + u[id[n],2]) * xInt[n] +
                 (beta[3] + u[id[n],3]) * xCond[n] +
                 (beta[4] + u[id[n],4]) * xIntByCond[n];
    y ~ bernoulli_logit(mu);
}