Stan model of Kliegl et al 2011, version 3

gistfile1.txt

data {
int<lower=1> N;
real rt[N];                         //outcome
real c1[N];                         //predictor
real c2[N];                         //predictor
real c3[N];                         //predictor
int<lower=1> I;                     //number of subjects
int<lower=1, upper=I> id[N];        //subject id
vector[4] mu_prior;                 //vector of zeros passed in from R
}
parameters {
vector[4] beta;                     // intercept and slope
vector[4] u[I];                     // random intercept and slopes
real<lower=0> sigma_e;              // residual sd 
vector<lower=0>[4] sigma_u;         // subj sd
corr_matrix[4] Omega;               // correlation matrix for random intercepts and slopes
}
transformed parameters {
    matrix[4,4] D;
    D <- diag_matrix(sigma_u);
}
model {
matrix[4,4] L;
matrix[4,4] DL;
real mu[N];                         // mu for likelihood
real zero; 

//priors:
beta ~ normal(0,50);
sigma_e ~ normal(0,100);
sigma_u ~ normal(0,100);
Omega ~ lkj_corr(4.0);

L <- cholesky_decompose(Omega);
for (m in 1:4) 
 for (n in 1:m) 
     DL[m,n] <- L[m,n] * sigma_u[m]; 
zero <- 0; 
for (m in 1:4) 
  for (n in (m+1):4) 
       DL[m,n] <- zero; 

for (i in 1:I)                             // loop for subj random effects
    u[i] ~ multi_normal_cholesky(mu_prior, DL);

// can't do this yet:
//u ~ multi_normal_cholesky(mu_prior, DL);

for (n in 1:N) {
mu[n] <- beta[1] + beta[2]*c1[n] + beta[3]*c2[n] + beta[4]*c3[n] +  u[id[n], 1] + u[id[n], 2]*c1[n] + u[id[n], 3]*c2[n] + u[id[n], 4]*c3[n];
}
rt ~ normal(mu,sigma_e);           // likelihood
}
generated quantities {
matrix[4,4] L;
cov_matrix[4] Sigma;
real<lower=-1,upper=1> rho12;
real<lower=-1,upper=1> rho13;
real<lower=-1,upper=1> rho14;
real<lower=-1,upper=1> rho23;
real<lower=-1,upper=1> rho24;
real<lower=-1,upper=1> rho34;

Sigma <- L * L'; 

rho12 <- Sigma[1,2] / sqrt(Sigma[1,1] * Sigma[2,2]);
rho13 <- Sigma[1,3] / sqrt(Sigma[1,1] * Sigma[3,3]);
rho14 <- Sigma[1,4] / sqrt(Sigma[1,1] * Sigma[4,4]);
rho23 <- Sigma[2,3] / sqrt(Sigma[2,2] * Sigma[3,3]);
rho24 <- Sigma[2,4] / sqrt(Sigma[2,2] * Sigma[4,4]);
rho34 <- Sigma[3,4] / sqrt(Sigma[3,3] * Sigma[4,4]);
}

The paper related to this data is:

http://read.psych.uni-potsdam.de/pmr2/index.php?option=com_content&view=article&id=10:kliegt-type&catid=11:publications&Itemid=13