References

Load packages

library(magrittr)
library(rstan)
library(lme4)

Linear Mixed Effects Model

Random slope model

\(Y_ij = \beta_{0j} + \beta_1 X_{1ij} + e_{0ij}\)

where

\(\beta_{0j} = \beta_0 + u_{0j}\)

\(u_{0j} \sim N(0, \sigma_{u_0}^2)\)

\(e_{0ij} \sim N(0, \sigma_{e_0}^2)\)

## Data generation
set.seed(201503031)
## Random intercepts
sigma_u0 <- 1
u_0j     <- rnorm(n = 10, mean = 0, sd = sigma_u0)
## Population intercept
beta_0   <- 3
## Varying intercepts
beta_0j  <- rep(beta_0 + u_0j, 10)
## Cluster indicator
cluster  <- rep(1:10, each = 10)
## Individual-level error
sigma    <- 2
e_0ij    <- rnorm(n = 100, mean = 0, sd = sigma)
##
X_1ij    <- runif(n = 100, min = 0, max = 20)
## Population slope
beta_1   <- 1.5
## Outcome
Y_ij     <- beta_0j + beta_1*X_1ij + e_0ij
## Combine
dat      <- list(Ni      = length(Y_ij),
                 Nj      = length(unique(cluster)),
                 cluster = cluster,
                 Y_ij    = Y_ij,
                 X_1ij   = X_1ij)
## Load Stan file
fileName <- "./multilevel1.stan"
stan_code <- readChar(fileName, file.info(fileName)$size)
cat(stan_code)
## data {
##   // Define variables in data
##   // Number of observations (an integer)
##   int<lower=0> Ni;
##   // Number of clusters (an integer)
##   int<lower=0> Nj;
##   // Cluster IDs
##   int<lower=1> cluster[Ni];
##   // Continuous outcome
##   real Y_ij[Ni];
##   // Continuous predictor
##   real X_1ij[Ni];
## }
## 
## parameters {
##   // Define parameters to estimate
##   // Population intercept (a real number)
##   real beta_0;
##   // Random effects
##   real u_0j[Nj];
##   // Standard deviation of random effects
##   real<lower=0> sigma_u0;
## 
##   // Population slope
##   real beta_1;
## 
##   // Population standard deviation
##   real<lower=0> sigma;
## }
## 
## transformed parameters  {
##   // Varying intercepts
##   real beta_0j[Nj];
##   // Individual mean
##   real mu[Ni];
## 
##   // Varying intercepts definition
##   for (j in 1:Nj) {
##     beta_0j[j] <- beta_0 + u_0j[j];
##   }
## 
##   // Individual mean
##   for (i in 1:Ni) {
##     mu[i] <- beta_0j[cluster[i]] + beta_1 * X_1ij[i];
##   }
## }
## 
## model {
##   // Prior part of Bayesian inference
##   // Flat prior for mu (no need to specify if non-informative)
## 
##   // Random effects distribution
##   u_0j ~ normal(0, sigma_u0);
## 
##   // Likelihood part of Bayesian inference
##   // Outcome model N(mu, sigma^2) (use SD rather than Var)
##   for (i in 1:Ni) {
##     Y_ij[i] ~ normal(mu[i], sigma);    
##   }
## }
## Run Stan
resStan <- stan(model_code = stan_code, data = dat,
                chains = 3, iter = 3000, warmup = 500, thin = 10)
traceplot(resStan, pars = c("beta_0","beta_1","sigma_u0","sigma"), inc_warmup = FALSE)

print(resStan, pars = c("beta_0","beta_1","sigma_u0","sigma"))
## 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   2.48    0.02 0.52 1.47 2.13 2.47 2.82  3.54   586 1.00
## beta_1   1.49    0.00 0.05 1.39 1.46 1.49 1.52  1.57   576 1.00
## sigma_u0 0.57    0.02 0.39 0.09 0.26 0.48 0.77  1.55   356 1.01
## sigma    2.59    0.01 0.19 2.24 2.46 2.58 2.72  2.99   750 1.00
## 
## Samples were drawn using NUTS(diag_e) at Tue Mar  3 21:15:20 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_ij ~ X_1ij + (1 | cluster))
summary(resLmer)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Y_ij ~ X_1ij + (1 | cluster)
## 
## REML criterion at convergence: 477.5
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.47891 -0.74770  0.03418  0.56005  2.56016 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  cluster  (Intercept) 0.1231   0.3509  
##  Residual             6.6166   2.5723  
## Number of obs: 100, groups:  cluster, 10
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  2.47275    0.47557    5.20
## X_1ij        1.48818    0.04444   33.49
## 
## Correlation of Fixed Effects:
##       (Intr)
## X_1ij -0.808
## True
list(beta_0 = beta_0, beta_1 = beta_1, sigma_u0 = sigma_u0, sigma = sigma)
## $beta_0
## [1] 3
## 
## $beta_1
## [1] 1.5
## 
## $sigma_u0
## [1] 1
## 
## $sigma
## [1] 2