I want to give a quick tutorial on fitting Linear Mixed Models (hierarchical models) with a full variance-covariance matrix for random effects (what <a href="http://www.sciencedirect.com/science/article/pii/S0749596X12001180">Barr et al 2013</a> call a maximal model) using Stan. <br><br>For a longer version of this tutorial, see: <a href="http://www.ling.uni-potsdam.de/~vasishth/statistics/BayesLMMs.html">Sorensen, Hohenstein, Vasishth, 2016</a>. <br><br><strong>Prerequisites</strong>: You need to have <a href="https://cran.r-project.org/">R</a> and preferably <a href="https://www.rstudio.com/">RStudio</a> installed; RStudio is optional. You need to have rstan installed. See <a href="https://github.com/stan-dev/rstan/wiki/RStan-Getting-Started">here</a>. I am also assuming you have fit lmer models like these before: <pre><br>lmer(log(rt) ~ 1+RCType+dist+int+(1+RCType+dist+int|subj) + (1+RCType+dist+int|item), dat)<br></pre> If you don't know what the above code means, first read chapter 4 of <a href="https://github.com/vasishth/Statistics-lecture-notes-Potsdam/blob/master/IntroductoryStatistics/StatisticsNotesVasishth.pdf">my lecture notes</a>. <h2> The code and data format needed to fit LMMs in Stan</h2> <h3>The data</h3> I assume you have a 2x2 repeated measures design with some continuous measure like reading time (rt) data and want to do a main effects and interaction contrast coding. Let's say your main effects are RCType and dist, and the interaction is coded as int. All these contrast codings are $\pm 1$. If you don't know what contrast coding is, see <a href="https://github.com/vasishth/Statistics-lecture-notes-Potsdam/blob/master/IntroductoryStatistics/StatisticsNotesVasishth.pdf">these notes</a> and read section 4.3 (although it's best to read the whole chapter). I am using an excerpt of an example data-set from <a href="http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0100986">Husain et al. 2014</a>. <pre><br>&quot;subj&quot; &quot;item&quot; &quot;rt&quot;&quot;RCType&quot; &quot;dist&quot; &quot;int&quot;<br>1       14    438  -1        -1      1<br>1       16    531   1        -1     -1<br>1       15    422   1         1      1<br>1       18   1000  -1        -1      1 <br>...<br></pre> Assume that these data are stored in R as a data-frame with name rDat. <h3> The Stan code</h3> Copy the following Stan code into a text file and save it as the file matrixModel.stan. For continuous data like reading times or EEG, you never need to touch this file again. You will only ever specify the design matrix X and the structure of the data. The rest is all taken care of. <pre><br>data {<br>  int N;               //no trials<br>  int P;               //no fixefs<br>  int J;               //no subjects<br>  int n_u;             //no subj ranefs<br>  int K;               //no items<br>  int n_w;             //no item ranefs<br>  int subj[N]; //subject indicator<br>  int item[N]; //item indicator<br>  row_vector[P] X[N];           //fixef design matrix<br>  row_vector[n_u] Z_u[N];       //subj ranef design matrix<br>  row_vector[n_w] Z_w[N];       //item ranef design matrix<br>  vector[N] rt;                 //reading time<br>}<br><br>parameters {<br>  vector[P] beta;               //fixef coefs<br>  cholesky_factor_corr[n_u] L_u;  //cholesky factor of subj ranef corr matrix<br>  cholesky_factor_corr[n_w] L_w;  //cholesky factor of item ranef corr matrix<br>  vector[n_u] sigma_u; //subj ranef std<br>  vector[n_w] sigma_w; //item ranef std<br>  real sigma_e;        //residual std<br>  vector[n_u] z_u[J];           //spherical subj ranef<br>  vector[n_w] z_w[K];           //spherical item ranef<br>}<br><br>transformed parameters {<br>  vector[n_u] u[J];             //subj ranefs<br>  vector[n_w] w[K];             //item ranefs<br>  {<br>    matrix[n_u,n_u] Sigma_u;    //subj ranef cov matrix<br>    matrix[n_w,n_w] Sigma_w;    //item ranef cov matrix<br>    Sigma_u = diag_pre_multiply(sigma_u,L_u);<br>    Sigma_w = diag_pre_multiply(sigma_w,L_w);<br>    for(j in 1:J)<br>      u[j] = Sigma_u * z_u[j];<br>    for(k in 1:K)<br>      w[k] = Sigma_w * z_w[k];<br>  }<br>}<br><br>model {<br>  //priors<br>  L_u ~ lkj_corr_cholesky(2.0);<br>  L_w ~ lkj_corr_cholesky(2.0);<br>  for (j in 1:J)<br>    z_u[j] ~ normal(0,1);<br>  for (k in 1:K)<br>    z_w[k] ~ normal(0,1);<br>  //likelihood<br>  for (i in 1:N)<br>    rt[i] ~ lognormal(X[i] * beta + Z_u[i] * u[subj[i]] + Z_w[i] * w[item[i]], sigma_e);<br>}<br></pre> <h3>Define the design matrix</h3> Since we want to test the main effects coded as the columns RCType, dist, and int, our design matrix will look like this: <pre><br># Make design matrix<br>X attr(X, &quot;assign&quot;) </pre> <h4>Prepare data for Stan</h4> Stan expects the data in a list form, not as a data frame (unlike lmer). So we set it up as follows: <pre><br># Make Stan data<br>stanDat P = ncol(X),<br>n_u = ncol(X),<br>n_w = ncol(X),<br>X = X,<br>Z_u = X,<br>Z_w = X,<br>J = nlevels(rDat$subj),<br>K = nlevels(rDat$item),<br>rt = rDat$rt,<br>subj = as.integer(rDat$subj),<br>item = as.integer(rDat$item))<br></pre> <h3>Load library rstan and fit Stan model</h3> <pre><br>library(rstan) <br>rstan_options(auto_write = TRUE)<br>options(mc.cores = parallel::detectCores())<br><br># Fit the model<br>matrixFit iter = 2000, chains = 4)<br></pre> <h4>Examine posteriors</h4> <pre><br>print(matrixFit)<br></pre> This print output is overly verbose. I wrote a simple function to get the essential information quickly. <pre><br>stan_results  m_extr  par_names  means  quantiles                    function(x)quantile(x,probs=c(0.025,0.975)))<br>  means  quants  summry  colnames(summry)  summry<br>}<br></pre> For example, if I want to see only the posteriors of the four beta parameters, I can write: <pre><br>stan_results(matrixFit, params=c(&quot;beta[1]&quot;,&quot;beta[2]&quot;,&quot;beta[3]&quot;,&quot;beta[4]&quot;))<br></pre> For more details, such as interpreting the results and computing things like Bayes Factors, see <a href="http://onlinelibrary.wiley.com/doi/10.1111/lnc3.12207/abstract;jsessionid=EF48E0A9341AD02594F99C034EF41DDF.f03t01">Nicenboim and Vasishth 2016</a>. <h4>FAQ: What if I don't want to fit a lognormal?</h4> In the Stan code above, I assume a lognormal function for the reading times: <pre><br> rt[i] ~ lognormal(X[i] * beta + Z_u[i] * u[subj[i]] + Z_w[i] * w[item[i]], sigma_e);<br></pre> If this upsets you deeply and you want to use a normal distribution (and in fact, for EEG data this makes sense), go right ahead and change the lognormal to normal: <pre><br> rt[i] ~ normal(X[i] * beta + Z_u[i] * u[subj[i]] + Z_w[i] * w[item[i]], sigma_e);<br></pre> <h4>FAQ: What if I my dependent measure is binary (0,1) responses?</h4> Use this Stan code instead of the one shown above. Here, I assume that you have a column called response in the data, which has 0,1 values. These are the trial level binary responses. <pre><br>data {<br>  int N;               //no trials<br>  int P;               //no fixefs<br>  int J;               //no subjects<br>  int n_u;             //no subj ranefs<br>  int K;               //no items<br>  int n_w;             //no item ranefs<br>  int subj[N]; //subject indicator<br>  int item[N]; //item indicator<br>  row_vector[P] X[N];           //fixef design matrix<br>  row_vector[n_u] Z_u[N];       //subj ranef design matrix<br>  row_vector[n_w] Z_w[N];       //item ranef design matrix<br>  int response[N];                 //response<br>}<br><br>parameters {<br>  vector[P] beta;               //fixef coefs<br>  cholesky_factor_corr[n_u] L_u;  //cholesky factor of subj ranef corr matrix<br>  cholesky_factor_corr[n_w] L_w;  //cholesky factor of item ranef corr matrix<br>  vector[n_u] sigma_u; //subj ranef std<br>  vector[n_w] sigma_w; //item ranef std<br>  vector[n_u] z_u[J];           //spherical subj ranef<br>  vector[n_w] z_w[K];           //spherical item ranef<br>}<br><br>transformed parameters {<br>  vector[n_u] u[J];             //subj ranefs<br>  vector[n_w] w[K];             //item ranefs<br>  {<br>    matrix[n_u,n_u] Sigma_u;    //subj ranef cov matrix<br>    matrix[n_w,n_w] Sigma_w;    //item ranef cov matrix<br>    Sigma_u = diag_pre_multiply(sigma_u,L_u);<br>    Sigma_w = diag_pre_multiply(sigma_w,L_w);<br>    for(j in 1:J)<br>      u[j] = Sigma_u * z_u[j];<br>    for(k in 1:K)<br>      w[k] = Sigma_w * z_w[k];<br>  }<br>}<br><br>model {<br>  //priors<br>  beta ~ cauchy(0,2.5);<br>  sigma_u ~ cauchy(0,2.5);<br>  sigma_w ~ cauchy(0,2.5);<br>  L_u ~ lkj_corr_cholesky(2.0);<br>  L_w ~ lkj_corr_cholesky(2.0);<br>  for (j in 1:J)<br>    z_u[j] ~ normal(0,1);<br>  for (k in 1:K)<br>    z_w[k] ~ normal(0,1);<br>  //likelihood<br>  for (i in 1:N)<br>    response[i] ~ bernoulli_logit(X[i] * beta + Z_u[i] * u[subj[i]] + Z_w[i] * w[item[i]]);<br>}<br></pre> <h4>For reproducible example code</h4> See <a href="http://www.ling.uni-potsdam.de/~vasishth/statistics/BayesLMMs.html">here</a>.