Three-level nested linear model in Stan (Classroom empirical data)
kaz_yos
2015-03-07
References
Harvard Chan School SBS263 Multilevel analysis, Lecture 6 Random intercepts
Harvard Chan School SBS263 Multilevel analysis, Lecture 15 Three-level model and Review
Classroom example data: http://www-personal.umich.edu/~bwest/chapter4.html
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)
print(resStan, pars = c("beta_0","sigma_e0","sigma_u0jk","sigma_u0k"))## Inference for Stan model: stan_code.
## 4 chains, each with iter=10000; warmup=1000; thin=10;
## post-warmup draws per chain=900, total post-warmup draws=3600.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## beta_0 57.44 0.03 1.46 54.58 56.45 57.43 58.41 60.29 3338 1.00
## sigma_e0 32.14 0.02 0.77 30.69 31.61 32.13 32.65 33.71 2111 1.00
## sigma_u0jk 9.93 0.10 2.19 5.36 8.53 10.01 11.41 14.09 463 1.01
## sigma_u0k 8.71 0.06 2.01 4.58 7.41 8.80 10.05 12.54 1179 1.00
##
## Samples were drawn using NUTS(diag_e) at Sat Mar 7 23:08:56 2015.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).lmer
resLmer <- lmer(mathgain ~ (1 | schoolid/classid), data = classroom)
summary(resLmer)## Linear mixed model fit by REML ['lmerMod']
## Formula: mathgain ~ (1 | schoolid/classid)
## Data: classroom
##
## REML criterion at convergence: 11768.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.6441 -0.5984 -0.0336 0.5334 5.6335
##
## Random effects:
## Groups Name Variance Std.Dev.
## classid:schoolid (Intercept) 99.23 9.961
## schoolid (Intercept) 77.49 8.803
## Residual 1028.23 32.066
## Number of obs: 1190, groups: classid:schoolid, 312; schoolid, 107
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 57.427 1.443 39.79Linear Mixed Effects Model (random intercepts with fixed effect covariates)
## Design matrix for model 4.4
desMat <- model.matrix(object = ~ 1 + mathkind + sex + minority + ses + housepov, data = classroom)
## Combine as a stan dataset
dat2 <- with(classroom,
list(Ni = length(unique(childid)),
Nj = length(unique(classid)),
Nk = length(unique(schoolid)),
p = ncol(desMat),
desMat = desMat,
classid = classid,
schoolid = schoolid,
schoolLookup = schoolLookupVec,
mathgain = mathgain))## Load Stan file
fileName <- "./multilevel3.2.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;
## // Number of fixed effect parameters
## int<lower=0> p;
##
## // Design matrix
## real desMat[Ni,p];
##
## // 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
## // Fixed effects
## real beta[p];
##
## // 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[1] + 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]] +
## desMat[i,2]*beta[2] + desMat[i,3]*beta[3] + desMat[i,4]*beta[4] +
## desMat[i,5]*beta[5] + desMat[i,6]*beta[6];
## }
## }
##
## 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 = dat2,
chains = 4, iter = 10000, warmup = 1000, thin = 10)traceplot(resStan, pars = c("beta","sigma_e0","sigma_u0jk","sigma_u0k"), inc_warmup = FALSE)
print(resStan, pars = c("beta","sigma_e0","sigma_u0jk","sigma_u0k"))## Inference for Stan model: stan_code.
## 4 chains, each with iter=10000; warmup=1000; thin=10;
## post-warmup draws per chain=900, total post-warmup draws=3600.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## beta[1] 284.99 0.18 11.03 262.90 277.74 285.35 292.52 306.06 3600 1.00
## beta[2] -0.47 0.00 0.02 -0.51 -0.49 -0.47 -0.46 -0.43 3600 1.00
## beta[3] -1.22 0.03 1.67 -4.44 -2.40 -1.22 -0.12 2.16 3600 1.00
## beta[4] -7.75 0.04 2.37 -12.32 -9.34 -7.80 -6.22 -3.17 3339 1.00
## beta[5] 5.23 0.02 1.25 2.82 4.37 5.21 6.10 7.68 3042 1.00
## beta[6] -11.51 0.17 10.05 -31.02 -18.30 -11.44 -4.68 8.19 3600 1.00
## sigma_e0 27.17 0.02 0.66 25.95 26.73 27.14 27.60 28.51 902 1.00
## sigma_u0jk 8.83 0.16 2.05 3.88 7.78 8.99 10.18 12.23 155 1.02
## sigma_u0k 8.84 0.04 1.57 5.70 7.83 8.86 9.86 11.89 1463 1.00
##
## Samples were drawn using NUTS(diag_e) at Sun Mar 8 03:21:51 2015.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).lmer
resLmer <- lmer(mathgain ~ mathkind + sex + minority + ses + housepov + (1 | schoolid/classid), data = classroom)
summary(resLmer)## Linear mixed model fit by REML ['lmerMod']
## Formula: mathgain ~ mathkind + sex + minority + ses + housepov + (1 | schoolid/classid)
## Data: classroom
##
## REML criterion at convergence: 11378.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.7926 -0.6052 -0.0327 0.5524 4.2727
##
## Random effects:
## Groups Name Variance Std.Dev.
## classid:schoolid (Intercept) 81.56 9.031
## schoolid (Intercept) 77.76 8.818
## Residual 734.42 27.100
## Number of obs: 1190, groups: classid:schoolid, 312; schoolid, 107
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 285.05797 11.02077 25.866
## mathkind -0.47086 0.02228 -21.133
## sex -1.23459 1.65743 -0.745
## minority -7.75587 2.38499 -3.252
## ses 5.23971 1.24497 4.209
## housepov -11.43920 9.93736 -1.151
##
## Correlation of Fixed Effects:
## (Intr) mthknd sex minrty ses
## mathkind -0.969
## sex -0.042 -0.032
## minority -0.265 0.153 -0.015
## ses 0.123 -0.165 0.019 0.144
## housepov -0.173 0.035 -0.009 -0.184 0.078- Top Page: http://rpubs.com/kaz_yos/
- github: https://github.com/kaz-yos