Three-level nested linear model in Stan
kaz_yos
2015-03-03
References
Harvard Chan School SBS263 Multilevel analysis, Lecture 6 Random intercepts
Harvard Chan School SBS263 Multilevel analysis, Lecture 15 Three-level model and Review
Multilevel Analysis for Applied Research. It’s just Regression!
Load packages
library(magrittr)
library(rstan)
library(lme4)Linear Mixed Effects Model
Random intercept model
Level 1 model: \(Y_{ijk} = \beta_{0jk} + \beta_1 X_{1ijk} + e_{0ijk}\)
Level 1 errors: \(e_{0ijk} \sim N(0, \sigma_{e_0}^2)\)
Level 2 model: \(\beta_{0jk} = \beta_{0k} + u_{0jk}\)
Level 2 errors: \(u_{0jk} \sim N(0, \sigma_{u_{0jk}}^2)\)
Level 3 model: \(\beta_{0k} = \beta_0 + u_{0k}\)
Level 3 errors: \(u_{0k} \sim N(0, \sigma_{u_{0k}}^2)\)
Data generation
## Data generation
set.seed(201503031)
## Population intercept
beta_0 <- 3
## Population slope
beta_1 <- 1.5
## Level-1 errors
sigma_e0 <- 3
e_0ijk <- rnorm(n = 10*10*10, mean = 0, sd = sigma_e0)
level1 <- seq_len(10*10*10)
## Level-2 errors
sigma_u0jk <- 2
u_0jk <- rnorm(n = 10*10, mean = 0, sd = sigma_u0jk)
level2 <- rep(seq_len(10*10), each = 10)
## Level-3 errors
sigma_u0k <- 1
u_0k <- rnorm(n = 10, mean = 0, sd = sigma_u0k)
level3 <- rep(seq_len(10), each = 100)
## Varying intercepts
## Level 3 (repeat for level 2)
beta_0k <- rep(beta_0 + u_0k, each = 10)
## Level 2 (repeat for level 1)
beta_0jk <- rep(beta_0k + u_0jk, each = 10)
## Predictor
X_1ijk <- runif(n = 1000, min = 0, max = 100)
## Outcome
Y_ijk <- beta_0jk + beta_1*X_1ijk + e_0ijk
## Combine
dat <- list(Ni = length(unique(level1)),
Nj = length(unique(level2)),
Nk = length(unique(level3)),
level2 = level2,
level3 = level3,
lev3ForLev2 = rep(seq_len(10), each = 10),
Y_ijk = Y_ijk,
X_1ijk = X_1ijk)Stan
## Load Stan file
fileName <- "./multilevel2.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> level2[Ni];
## int<lower=1> level3[Ni];
## int<lower=1> lev3ForLev2[Nj];
##
## // Continuous outcome
## real Y_ijk[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
## 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[lev3ForLev2[j]] + u_0jk[j];
## }
## // Individual mean
## for (i in 1:Ni) {
## mu[i] <- beta_0jk[level2[i]] + beta_1 * X_1ijk[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) {
## Y_ijk[i] ~ normal(mu[i], sigma_e0);
## }
## }## Run Stan
resStan <- stan(model_code = stan_code, data = dat,
chains = 3, iter = 3000, warmup = 500, thin = 10)resStanExt <- rstan::extract(resStan, permuted = TRUE)
rstan::traceplot(resStan, pars = c("beta_0","beta_1","sigma_e0","sigma_u0jk","sigma_u0k"), inc_warmup = FALSE)
print(resStan, pars = c("beta_0","beta_1","sigma_e0","sigma_u0jk","sigma_u0k"))## Inference for Stan model: stan_code.
## 3 chains, each with iter=3000; warmup=500; thin=10;
## post-warmup draws per chain=250, total post-warmup draws=750.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## beta_0 4.23 0.02 0.49 3.18 3.93 4.22 4.56 5.19 689 1
## beta_1 1.50 0.00 0.00 1.49 1.50 1.50 1.50 1.51 750 1
## sigma_e0 2.96 0.00 0.07 2.82 2.91 2.96 3.00 3.09 750 1
## sigma_u0jk 1.84 0.01 0.19 1.51 1.71 1.82 1.96 2.24 750 1
## sigma_u0k 1.20 0.02 0.46 0.55 0.89 1.12 1.39 2.37 684 1
##
## Samples were drawn using NUTS(diag_e) at Sat Mar 7 19:51:37 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(Y_ijk ~ X_1ijk + (1 | level3/level2))
varCorrMat <- as.data.frame(VarCorr(resLmer))
summary(resLmer)## Linear mixed model fit by REML ['lmerMod']
## Formula: Y_ijk ~ X_1ijk + (1 | level3/level2)
##
## REML criterion at convergence: 5180.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0604 -0.6424 0.0027 0.6755 3.5638
##
## Random effects:
## Groups Name Variance Std.Dev.
## level2:level3 (Intercept) 3.294 1.815
## level3 (Intercept) 1.087 1.043
## Residual 8.725 2.954
## Number of obs: 1000, groups: level2:level3, 100; level3, 10
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 4.212844 0.425965 9.9
## X_1ijk 1.500429 0.003429 437.6
##
## Correlation of Fixed Effects:
## (Intr)
## X_1ijk -0.414Comparison
## True
resDf <- data.frame(beta_0 = c(beta_0, median(resStanExt$beta_0), coef(summary(resLmer))[1,1]),
beta_1 = c(beta_1, median(resStanExt$beta_1), coef(summary(resLmer))[2,1]),
sigma_e0 = c(sigma_e0, median(resStanExt$sigma_e0), varCorrMat[3,"sdcor"]),
sigma_u0jk = c(sigma_u0jk, median(resStanExt$sigma_u0jk), varCorrMat[1,"sdcor"]),
sigma_u0k = c(sigma_u0k, median(resStanExt$sigma_u0k), varCorrMat[2,"sdcor"]))
rownames(resDf) <- c("true","stan","lmer")
print(resDf, digits = 3)## beta_0 beta_1 sigma_e0 sigma_u0jk sigma_u0k
## true 3.00 1.5 3.00 2.00 1.00
## stan 4.22 1.5 2.96 1.82 1.12
## lmer 4.21 1.5 2.95 1.81 1.04- Top Page: http://rpubs.com/kaz_yos/
- github: https://github.com/kaz-yos