References

Load packages

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.

Load data

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"
stan_code <- readChar(fileName, file.info(fileName)$size)
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)