## References

``````library(magrittr)
library(rstan)
library(lme4)
library(dplyr)``````

## Linear Mixed Effects Model (intercept only)

childid is nested within classid, which in turn is nested within schoolid.

``````classroom <- read.csv("http://www-personal.umich.edu/~bwest/classroom.csv")

## Sort by class ID for easier handling in Stan
classroom <- arrange(classroom, classid, schoolid)

## Create a vector of school IDs where j-th element gives school ID for class ID j
schoolLookupVec <- unique(classroom[c("classid","schoolid")])[,"schoolid"]

## Combine as a stan dataset
dat <- with(classroom,
list(Ni           = length(unique(childid)),
Nj           = length(unique(classid)),
Nk           = length(unique(schoolid)),
classid      = classid,
schoolid     = schoolid,
schoolLookup = schoolLookupVec,
mathgain     = mathgain))``````
``````## Load Stan file
fileName <- "./multilevel3.stan"
cat(stan_code)``````
``````## data {
##   // Define variables in data
##   // Number of level-1 observations (an integer)
##   int<lower=0> Ni;
##   // Number of level-2 clusters
##   int<lower=0> Nj;
##   // Number of level-3 clusters
##   int<lower=0> Nk;
##
##   // Cluster IDs
##   int<lower=1> classid[Ni];
##   int<lower=1> schoolid[Ni];
##
##   // Level 3 look up vector for level 2
##   int<lower=1> schoolLookup[Nj];
##
##   // Continuous outcome
##   real mathgain[Ni];
##
##   // Continuous predictor
##   // real X_1ijk[Ni];
## }
##
## parameters {
##   // Define parameters to estimate
##   // Population intercept (a real number)
##   real beta_0;
##   // Population slope
##   // real beta_1;
##
##   // Level-1 errors
##   real<lower=0> sigma_e0;
##
##   // Level-2 random effect
##   real u_0jk[Nj];
##   real<lower=0> sigma_u0jk;
##
##   // Level-3 random effect
##   real u_0k[Nk];
##   real<lower=0> sigma_u0k;
## }
##
## transformed parameters  {
##   // Varying intercepts
##   real beta_0jk[Nj];
##   real beta_0k[Nk];
##
##   // Individual mean
##   real mu[Ni];
##
##   // Varying intercepts definition
##   // Level-3 (10 level-3 random intercepts)
##   for (k in 1:Nk) {
##     beta_0k[k] <- beta_0 + u_0k[k];
##   }
##   // Level-2 (100 level-2 random intercepts)
##   for (j in 1:Nj) {
##     beta_0jk[j] <- beta_0k[schoolLookup[j]] + u_0jk[j];
##   }
##   // Individual mean
##   for (i in 1:Ni) {
##     mu[i] <- beta_0jk[classid[i]];
##   }
## }
##
## model {
##   // Prior part of Bayesian inference
##   // Flat prior for mu (no need to specify if non-informative)
##
##   // Random effects distribution
##   u_0k  ~ normal(0, sigma_u0k);
##   u_0jk ~ normal(0, sigma_u0jk);
##
##   // Likelihood part of Bayesian inference
##   // Outcome model N(mu, sigma^2) (use SD rather than Var)
##   for (i in 1:Ni) {
##     mathgain[i] ~ normal(mu[i], sigma_e0);
##   }
## }``````
``````## Run Stan
resStan <- stan(model_code = stan_code, data = dat,
chains = 4, iter = 10000, warmup = 1000, thin = 10)``````
``traceplot(resStan, pars = c("beta_0","sigma_e0","sigma_u0jk","sigma_u0k"), inc_warmup = FALSE)``