rstanarm: Mixed Model

Let’s look at a mixed model for another demonstration

• The average reaction time per day for subjects in a sleep deprivation study
• On day 0 the subjects had their normal amount of sleep
• Subsequently restricted to 3 hours of sleep per night
• The observations represent the average reaction time on a series of tests

We’ll have a random intercept and random coefficient for Days

library(lme4)
sleepstudy_lmer <- lmer(Reaction ~ Days + (1 + Days|Subject), 
                        data = sleepstudy)
summary(sleepstudy_lmer)

Linear mixed model fit by REML ['lmerMod']
Formula: Reaction ~ Days + (1 + Days | Subject)
   Data: sleepstudy

REML criterion at convergence: 1743.6

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.9536 -0.4634  0.0231  0.4634  5.1793 

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 Subject  (Intercept) 612.09   24.740       
          Days         35.07    5.922   0.07
 Residual             654.94   25.592       
Number of obs: 180, groups:  Subject, 18

Fixed effects:
            Estimate Std. Error t value
(Intercept)  251.405      6.825  36.838
Days          10.467      1.546   6.771

Correlation of Fixed Effects:
     (Intr)
Days -0.138

Again, rstanarm sticks with the same style

sleepstudy_blmer <- stan_lmer(Reaction ~ Days + (1 + Days|Subject), 
                              data = sleepstudy)
summary(sleepstudy_blmer, digits = 3)

stan_lmer
 family:       gaussian [identity]
 formula:      Reaction ~ Days + (1 + Days | Subject)
 observations: 180
------
            Median  MAD_SD 
(Intercept) 251.616   6.503
Days         10.451   1.629

Auxiliary parameter(s):
      Median MAD_SD
sigma 25.853  1.541

Error terms:
 Groups   Name        Std.Dev. Corr
 Subject  (Intercept) 24.258       
          Days         6.901   0.08
 Residual             25.959       
Num. levels: Subject 18 

Sample avg. posterior predictive distribution of y:
         Median  MAD_SD 
mean_PPD 298.572   2.716

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg

In the Bayesian model, the random effects are not BLUPS, but are parameters estimates in the model

In this case, we see a little more shrinkage relative to the standard approach

The following are obtained from the same ranef function used in lme4

lme4	bayesian
2.3	2.6
-40.4	-36.6
-39.0	-35.3
23.7	20.3
22.3	19.4
9.0	7.6
16.8	14.5
-7.2	-6.5
-0.3	-1.2
34.9	31.5
-25.2	-22.6
-13.1	-11.7
4.6	3.3
20.9	18.4
3.3	2.6
-25.6	-23.0
0.8	0.6
12.3	10.8