See the introduction for an overview. Also see other analyses of this data. See the textbook for a full discussion.
Load the libraries:
library(ggplot2)
library(rstan)
library(reshape2)
library(dplyr)Data
Load in and summarize the data:
data(vision, package="faraway")
vision$npower <- rep(1:4,14)
summary(vision) acuity power eye subject npower
Min. : 94 6/6 :14 left :28 1:8 Min. :1.00
1st Qu.:110 6/18:14 right:28 2:8 1st Qu.:1.75
Median :115 6/36:14 3:8 Median :2.50
Mean :113 6/60:14 4:8 Mean :2.50
3rd Qu.:118 5:8 3rd Qu.:3.25
Max. :124 6:8 Max. :4.00
7:8 Plot the data:
ggplot(vision, aes(y=acuity, x=npower, linetype=eye))+geom_line()+facet_wrap(~ subject)
Fitting
Fit the model. Requires use of STAN command file multilevel.stan. We have used uninformative priors for the treatment effects but slightly informative half-cauchy priors for the three variances. The fixed effects have been collected into a single design matrix. We are using the same STAN command file as for the multilevel.html multilevel example because we have the same two levels of nesting.
lmod <- lm(acuity ~ power, vision)
sdscal <- sd(residuals(lmod))
Xmatrix <- model.matrix(lmod)
vision$subjeye <- factor(paste(vision$subject,vision$eye,sep="."))
visdat <- list(Nobs=nrow(vision),
Npreds=ncol(Xmatrix),
Nlev1=length(unique(vision$subject)),
Nlev2=length(unique(vision$subjeye)),
y=vision$acuity,
x=Xmatrix,
levind1=as.numeric(vision$subject),
levind2=as.numeric(vision$subjeye),
sdscal=sdscal)Break the fitting of the model into three steps.
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
rt <- stanc("multilevel.stan")
sm <- stan_model(stanc_ret = rt, verbose=FALSE)
set.seed(123)
system.time(fit <- sampling(sm, data=visdat)) user system elapsed
0.351 0.088 8.101 Diagnostics
For the error SD:
pname <- "sigmaeps"
muc <- rstan::extract(fit, pars=pname, permuted=FALSE, inc_warmup=FALSE)
mdf <- melt(muc)
ggplot(mdf,aes(x=iterations,y=value,color=chains)) + geom_line() + ylab(mdf$parameters[1])