Multi-Level Modeling: Three Levels

Designs

• 3 level models are used when you multiple levels of nesting that you need to account for.
  • Students nested in classrooms, nested in schools
  • Patients nested in doctors, nested in hospitals

3 Levels

3 Level Model

Example

• Lets again examine active learning as it relates to math scores.
  • You measure students math scores (DV) and the proportion of time (IV) they spend using the computer (which you assign)
    • You expect that the more time they spend doing the active learning method, the higher their math test scores will be
  • You collected data in many classrooms in three different schools.
    • The number of classrooms can vary by schools, and the number of students per classroom can vary
    • You think your system might work best for classrooms with a lot of students (you cannot control for that, but you know how many students are in each classroom)
• Let’s go collect (simulate) our study (not shown here because of its length)
• Download Data

MLM.Data.L3<-read.csv("Mixed/Sim3level.csv")
# We need to set Classroom as a factor (since it was #)
MLM.Data.L3$Classroom.F<-as.factor(MLM.Data.L3$Classroom)

Sample Sizes

• Tabulate number of students per classroom per school

names(MLM.Data.L3)

## [1] "Math"        "ActiveTime"  "ClassSize"   "Classroom"   "School"     
## [6] "StudentID"   "Classroom.F"

CrossT <- table(MLM.Data.L3$School,MLM.Data.L3$Classroom.F)

• It seems that the data look somewhat balanced between schools, let’s just ignore the school-level (level 3)

• Let’s just see the overall effect and ignore the classroom-level (level 2)

library(ggplot2)
library(gridExtra)

theme_set(theme_bw(base_size = 7, base_family = "")) 

Active.Plot <-ggplot(data = MLM.Data.L3, aes(x = ActiveTime, y=Math))+
  coord_cartesian(ylim=c(10,80))+
  geom_point()+
  geom_smooth(method = "lm", se = TRUE)+
  xlab("Active Learning Time")+ylab("Math Score")+
  theme(legend.position = "none")

Size.Plot <-ggplot(data = MLM.Data.L3, aes(x = ClassSize, y=Math))+  
  coord_cartesian(ylim=c(10,80))+
  geom_point()+
  geom_smooth(method = "lm", se = TRUE)+
  xlab("Class Size")+ylab("Math Score")+
  theme(legend.position = "top")

 grid.arrange(Size.Plot, Active.Plot, ncol=2)

• Let’s just see: add in the classroom-level (level 2)
• But wait? Can we view Class size as a function of the classroom?
  • No because each classroom only as 1 class size

library(ggplot2)
theme_set(theme_bw(base_size = 10, base_family = "")) 

Active.Class <-ggplot(data = MLM.Data.L3, aes(x = ActiveTime, y=Math,group=Classroom.F))+
  coord_cartesian(ylim=c(10,80))+
  geom_point()+
  geom_smooth(method = "lm", se = TRUE, aes(colour = Classroom.F))+
  xlab("Active Learning Time")+ylab("Math Score")+
  theme(legend.position = "none")
Active.Class

Recap 2-Level MLM

• Let’s ignore the school level and deal only as if this were 2 level MLM
• So students nested in classrooms
• Problem is we have multiple classrooms with the same number: we can fix that in R easily

MLM.Data.L3$Classroom.F2<-paste(MLM.Data.L3$School,MLM.Data.L3$Classroom.F,sep=":")

Random intercept only Model

• Grand mean center

MLM.Data.L3$ActiveTime.GM<-scale(MLM.Data.L3$ActiveTime,scale=F)
MLM.Data.L3$ClassSize.GM<-scale(MLM.Data.L3$ClassSize,scale=F)

• Let test our cluster

library(lme4)
Model.0<-lmer(Math ~ 1
              +(1|Classroom.F2),  
              data=MLM.Data.L3, REML=FALSE)
summary(Model.0)

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Math ~ 1 + (1 | Classroom.F2)
##    Data: MLM.Data.L3
## 
##      AIC      BIC   logLik deviance df.resid 
##   3805.7   3818.7  -1899.8   3799.7      567 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5139 -0.6710 -0.0039  0.6477  2.7589 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  Classroom.F2 (Intercept) 108.41   10.412  
##  Residual                  37.22    6.101  
## Number of obs: 570, groups:  Classroom.F2, 30
## 
## Fixed effects:
##             Estimate Std. Error     df t value Pr(>|t|)    
## (Intercept)   43.990      1.919 29.977   22.93   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• Check to make sure the cluster explains any variance

library(performance)
icc(Model.0)

## # Intraclass Correlation Coefficient
## 
##     Adjusted ICC: 0.744
##   Unadjusted ICC: 0.744

• Let test our fixed effects

Model.1<-lmer(Math ~ ActiveTime.GM+ClassSize.GM
              +(1|Classroom.F2),  
              data=MLM.Data.L3, REML=FALSE)
summary(Model.1)

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Math ~ ActiveTime.GM + ClassSize.GM + (1 | Classroom.F2)
##    Data: MLM.Data.L3
## 
##      AIC      BIC   logLik deviance df.resid 
##   3390.1   3411.9  -1690.1   3380.1      565 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3148 -0.6516 -0.0386  0.6312  3.1178 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  Classroom.F2 (Intercept) 72.29    8.503   
##  Residual                 17.51    4.184   
## Number of obs: 570, groups:  Classroom.F2, 30
## 
## Fixed effects:
##               Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)    42.9594     1.5847  29.9349  27.109  < 2e-16 ***
## ActiveTime.GM  14.9403     0.6061 540.6478  24.651  < 2e-16 ***
## ClassSize.GM   -1.8671     0.4802  30.0467  -3.888 0.000518 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) AcT.GM
## ActiveTm.GM 0.000        
## ClassSiz.GM 0.167  0.001

• So we have main effects. Lets check for the interaction

Model.2<-lmer(Math ~ ActiveTime.GM*ClassSize.GM
              +(1|Classroom.F2),  
              data=MLM.Data.L3, REML=FALSE)
summary(Model.2)

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Math ~ ActiveTime.GM * ClassSize.GM + (1 | Classroom.F2)
##    Data: MLM.Data.L3
## 
##      AIC      BIC   logLik deviance df.resid 
##   3388.7   3414.8  -1688.3   3376.7      564 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3103 -0.6642 -0.0358  0.6455  3.1392 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  Classroom.F2 (Intercept) 72.41    8.509   
##  Residual                 17.40    4.171   
## Number of obs: 570, groups:  Classroom.F2, 30
## 
## Fixed effects:
##                            Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)                 42.9645     1.5858  29.9350  27.093  < 2e-16 ***
## ActiveTime.GM               14.9057     0.6044 540.6574  24.662  < 2e-16 ***
## ClassSize.GM                -1.8662     0.4805  30.0457  -3.884 0.000524 ***
## ActiveTime.GM:ClassSize.GM   0.3613     0.1939 540.5862   1.863 0.062999 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) AcT.GM ClS.GM
## ActiveTm.GM  0.000              
## ClassSiz.GM  0.167  0.001       
## AT.GM:CS.GM  0.002 -0.031  0.001

• The interaction is weak, model testing will probably show the same result

anova(Model.1,Model.2)

## Data: MLM.Data.L3
## Models:
## Model.1: Math ~ ActiveTime.GM + ClassSize.GM + (1 | Classroom.F2)
## Model.2: Math ~ ActiveTime.GM * ClassSize.GM + (1 | Classroom.F2)
##         npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
## Model.1    5 3390.1 3411.9 -1690.1   3380.1                       
## Model.2    6 3388.7 3414.8 -1688.3   3376.7 3.4595  1    0.06289 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Plot Random slopes

• Lets plot the random and fixed effects (per classroom)
• First, we will predict the model results

MLM.Data.L3$Model.2.FR<-predict(Model.2, newdata=MLM.Data.L3)

• Lets plot the raw data with the fitted of ActiveTime

theme_set(theme_bw(base_size = 7, base_family = "")) 

Active.Class.School <-ggplot(data = MLM.Data.L3, aes(x = ActiveTime, y=Math,group=Classroom.F))+
  coord_cartesian(ylim=c(10,80))+
  geom_point(aes(colour = Classroom.F))+
  geom_smooth(method = "lm", se = TRUE, aes(colour = Classroom.F))+
  xlab("Active Learning Time")+ylab("Math Score")+
  theme(legend.position = "none")

Active.Class.School.fit <-ggplot(data = MLM.Data.L3, aes(x = ActiveTime, y=Model.2.FR,group=Classroom.F))+
  coord_cartesian(ylim=c(10,80))+
  geom_point(aes(colour = Classroom.F))+
  geom_smooth(method = "lm", se = TRUE, aes(colour = Classroom.F))+
  xlab("Active Learning Time")+ylab("Math Score - Model 2")+
  theme(legend.position = "none")

 grid.arrange(Active.Class.School, Active.Class.School.fit, ncol=2)

• Lets add random slopes to this final model and see that happens.
• There are alot of classrooms and we can see in the figure that the slopes do not have a lot of variance

Random Coefficients Model

• Lets add the level 1 slope (Activetime) as function of the classroom

Model.2.RC<-lmer(Math ~ ActiveTime.GM*ClassSize.GM
              +(1+ActiveTime.GM|Classroom.F2),  
              data=MLM.Data.L3, REML=FALSE)
summary(Model.2.RC)

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: 
## Math ~ ActiveTime.GM * ClassSize.GM + (1 + ActiveTime.GM | Classroom.F2)
##    Data: MLM.Data.L3
## 
##      AIC      BIC   logLik deviance df.resid 
##   3358.3   3393.1  -1671.1   3342.3      562 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.04260 -0.65433 -0.03363  0.63223  3.07726 
## 
## Random effects:
##  Groups       Name          Variance Std.Dev. Corr 
##  Classroom.F2 (Intercept)   73.65    8.582         
##               ActiveTime.GM 22.19    4.711    -0.31
##  Residual                   15.34    3.916         
## Number of obs: 570, groups:  Classroom.F2, 30
## 
## Fixed effects:
##                            Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)                 42.8612     1.5984 29.9458  26.815  < 2e-16 ***
## ActiveTime.GM               14.8927     1.0450 30.3150  14.251 5.61e-15 ***
## ClassSize.GM                -1.8703     0.4843 30.0394  -3.862 0.000556 ***
## ActiveTime.GM:ClassSize.GM   0.3329     0.3263 33.1054   1.020 0.315013    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) AcT.GM ClS.GM
## ActiveTm.GM -0.258              
## ClassSiz.GM  0.167 -0.043       
## AT.GM:CS.GM -0.042  0.102 -0.250

• The interaction disappeared but the pattern of the other variables held
• However, we are ignoring school-level variable. Lets re-plot our data see if we missed anything

Plotting 3 Levels

• There are many options on how to do this
• Since we have only three schools, lets try facet plotting

theme_set(theme_bw(base_size = 7, base_family = "")) 

Active.Class.School <-ggplot(data = MLM.Data.L3, aes(x = ActiveTime, y=Math,group=Classroom.F))+
  facet_grid(~School)+
  coord_cartesian(ylim=c(10,80))+
  geom_point(aes(colour = Classroom.F))+
  geom_smooth(method = "lm", se = TRUE, aes(colour = Classroom.F))+
  xlab("Active Learning Time")+ylab("Math Score")+
  theme(legend.position = "none")

Size.School <-ggplot(data = MLM.Data.L3, aes(x = ClassSize, y=Math))+
  facet_grid(~School)+
  coord_cartesian(ylim=c(10,80))+
  geom_point()+
  geom_smooth(method = "lm", se = TRUE)+
  xlab("Class Size")+ylab("Math Score")+
  theme(legend.position = "top")

grid.arrange(Active.Class.School, Size.School, ncol=1, nrow=2)

• We can also visualize it this way (by School)

library(dplyr)
Plot.Means<-MLM.Data.L3 %>% group_by(School) %>%  
  dplyr::summarize(MathM=mean(Math, na.rm=TRUE),
           ClassSizeM=mean(ClassSize, na.rm=TRUE))

MathM.school<-ggplot(data = Plot.Means, aes(x = reorder(School, -MathM), y=MathM))+
  geom_point(aes(size = ClassSizeM))+
  xlab("")+ylab("Math Score")+
  theme_bw()+
  theme(legend.position = "top")
MathM.school

• You can see the size the dot changes as a function of the school!
• The school seems to make a huge difference
• The slope for Class Size seems to have vanished (Oh no, Simpson’s paradox again)
• We should control for school as a random intercept!

Confounded Variables

• It is clear from the figure for that Active learning spans the whole range for both classroom and school, but ClassSize is confounded
• Clearly the slopes vary a function of the school, but if we grand mean centering ClassSize it will be confounded with the school
  • In other words, ClassSize seems to vary as a function of the School if we test the fixed effect its going cross-over the schools (level 3 variable)
  • So whatever fixed effect we find could mean that is a function of the school (level 3) or really something about the variable? We cannot know!
  • So the only thing we do is center relative to the cluster
  • This NOT the case for Active learning because we assigned subjects to that condition
• To make this variable meaningful lets re-center ClassSize based on School
• Note: This variable now means something new! It means, “how does class size affect math score, relative the school you are in

library(plyr)
# Level 3 Cluster mean
MLM.Data.L3<-ddply(MLM.Data.L3,.(School), mutate, CS.School.Mean = mean(ClassSize))
MLM.Data.L3$ClassSize.SC<-MLM.Data.L3$ClassSize-MLM.Data.L3$CS.School.Mean

• See new plot

Size.School.2 <-ggplot(data = MLM.Data.L3, aes(x = ClassSize.SC, y=Math))+
  facet_grid(~School)+
  coord_cartesian(ylim=c(10,80))+
  geom_point()+
  geom_smooth(method = "lm", se = TRUE)+
  xlab("Class Size")+ylab("Math Score")+
  theme_bw()+
  theme(legend.position = "top")
Size.School.2

• Yes thats better, we removed the confound (but we need to be careful how we talk about it if we get any thing, but clearly from the graph we should not get much)

3 Level Random intercepts model

• We simply just need to add more random effects
  • +(1|School) = each school can have its own intercept
  • +(1|School:Classroom.F) = each classroom can have its own intercept relative to which school its nested in
  • Note: I used Classroom.F and not our new Classroom.F2.
    • Classroom.F is not implicitly nested in classroom (there were multiple classroom 1), thus we must tell the function that classroom were nested in schools
    • Classroom.F2is implicitly nested in classroom
    • Thus +(1|School) + (1|School:Classroom.F) = +(1|School) + (1|Classroom.F2)
    • The results will be all the same!

library(lme4)
L3.Model.0<-lmer(Math ~ 1
              +(1|School)
              +(1|School:Classroom.F),  
              data=MLM.Data.L3, REML=FALSE)
summary(L3.Model.0)

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Math ~ 1 + (1 | School) + (1 | School:Classroom.F)
##    Data: MLM.Data.L3
## 
##      AIC      BIC   logLik deviance df.resid 
##   3789.8   3807.2  -1890.9   3781.8      566 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5237 -0.6756 -0.0005  0.6546  2.7491 
## 
## Random effects:
##  Groups             Name        Variance Std.Dev.
##  School:Classroom.F (Intercept) 44.92    6.702   
##  School             (Intercept) 59.88    7.738   
##  Residual                       37.22    6.101   
## Number of obs: 570, groups:  School:Classroom.F, 30; School, 3
## 
## Fixed effects:
##             Estimate Std. Error     df t value Pr(>|t|)   
## (Intercept)   44.411      4.644  3.035   9.563  0.00231 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• Look at the random intercepts: Seems that school should not have been ignored
• Let verify via ICC

icc(L3.Model.0)

## # Intraclass Correlation Coefficient
## 
##     Adjusted ICC: 0.738
##   Unadjusted ICC: 0.738

• Clearly school level was important to explaining the data

• Let’s check our main effects now.

  • Class Size effect should go away based on what we saw in the graph!

L3.Model.1<-lmer(Math ~ ActiveTime.GM+ClassSize.SC
              +(1|School)
              +(1|School:Classroom.F),  
              data=MLM.Data.L3, REML=FALSE)
summary(L3.Model.1)

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: 
## Math ~ ActiveTime.GM + ClassSize.SC + (1 | School) + (1 | School:Classroom.F)
##    Data: MLM.Data.L3
## 
##      AIC      BIC   logLik deviance df.resid 
##   3385.9   3412.0  -1686.9   3373.9      564 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3123 -0.6496 -0.0515  0.6336  3.0742 
## 
## Random effects:
##  Groups             Name        Variance Std.Dev.
##  School:Classroom.F (Intercept) 44.69    6.685   
##  School             (Intercept) 60.42    7.773   
##  Residual                       17.51    4.184   
## Number of obs: 570, groups:  School:Classroom.F, 30; School, 3
## 
## Fixed effects:
##               Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)    44.4402     4.6607   3.0415   9.535  0.00231 ** 
## ActiveTime.GM  14.9488     0.6060 540.9288  24.669  < 2e-16 ***
## ClassSize.SC   -0.0306     0.5688  27.1017  -0.054  0.95749    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) AcT.GM
## ActiveTm.GM 0.001        
## ClassSiz.SC 0.031  0.007

• Interaction could still be real, lets check:

L3.Model.2<-lmer(Math ~ ActiveTime.GM*ClassSize.SC
              +(1|School)
              +(1|School:Classroom.F),  
              data=MLM.Data.L3, REML=FALSE)
summary(L3.Model.2)

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: 
## Math ~ ActiveTime.GM * ClassSize.SC + (1 | School) + (1 | School:Classroom.F)
##    Data: MLM.Data.L3
## 
##      AIC      BIC   logLik deviance df.resid 
##   3379.9   3410.3  -1682.9   3365.9      563 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.13965 -0.62420 -0.01136  0.62253  2.84136 
## 
## Random effects:
##  Groups             Name        Variance Std.Dev.
##  School:Classroom.F (Intercept) 45.03    6.711   
##  School             (Intercept) 60.36    7.769   
##  Residual                       17.25    4.153   
## Number of obs: 570, groups:  School:Classroom.F, 30; School, 3
## 
## Fixed effects:
##                            Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)                 44.4184     4.6597   3.0422   9.532  0.00231 ** 
## ActiveTime.GM               14.9967     0.6017 540.9271  24.926  < 2e-16 ***
## ClassSize.SC                -0.0368     0.5709  27.1003  -0.064  0.94908    
## ActiveTime.GM:ClassSize.SC  -0.8211     0.2893 540.8396  -2.839  0.00470 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) AcT.GM ClS.SC
## ActiveTm.GM  0.000              
## ClassSiz.SC  0.031  0.007       
## AT.GM:CS.SC  0.002 -0.028  0.004

• Significant interaction!
• But wait…

3 Level Random Coefficients model

Level 1 Random Slopes at Level 3

• Now let think this out for minute, we have level 1 slopes (ActiveTime) as a function of the classroom (level 2), but can they also vary as function of the school (level 3)?
• Lets look and see…

Active.School <-ggplot(data = MLM.Data.L3, aes(x = ActiveTime, y=Math, group=School))+
  geom_point(aes(color=School))+
  geom_smooth(method = "lm", se = TRUE,aes(color=School))+
  xlab("Active Learning Time")+ylab("Math Score")+
  theme_bw()+
  theme(legend.position = "none")
Active.School

• With only our three schools it seems not really (slopes all look mostly the same).
• Also, the random slopes will perfectly correlate with the random intercepts!
• See below:

L3.Model.2a<-lmer(Math ~ ActiveTime.GM*ClassSize.SC
              +(1+ActiveTime.GM|School)
              +(1|School:Classroom.F),  
              data=MLM.Data.L3, REML=FALSE)
summary(L3.Model.2a)

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Math ~ ActiveTime.GM * ClassSize.SC + (1 + ActiveTime.GM | School) +  
##     (1 | School:Classroom.F)
##    Data: MLM.Data.L3
## 
##      AIC      BIC   logLik deviance df.resid 
##   3360.5   3399.6  -1671.2   3342.5      561 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.86601 -0.66736 -0.04249  0.64750  2.62037 
## 
## Random effects:
##  Groups             Name          Variance Std.Dev. Corr 
##  School:Classroom.F (Intercept)   44.874   6.699         
##  School             (Intercept)   60.076   7.751         
##                     ActiveTime.GM  8.174   2.859    -1.00
##  Residual                         16.478   4.059         
## Number of obs: 570, groups:  School:Classroom.F, 30; School, 3
## 
## Fixed effects:
##                             Estimate Std. Error        df t value Pr(>|t|)   
## (Intercept)                 44.56614    4.64705   3.04601   9.590  0.00226 **
## ActiveTime.GM               14.35446    1.75474   3.07930   8.180  0.00347 **
## ClassSize.SC                -0.04435    0.56933  27.83552  -0.078  0.93847   
## ActiveTime.GM:ClassSize.SC  -0.81809    0.28292 539.72042  -2.892  0.00399 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) AcT.GM ClS.SC
## ActiveTm.GM -0.907              
## ClassSiz.SC  0.031  0.001       
## AT.GM:CS.SC  0.001 -0.008  0.004
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

• We could remove the correlation can capture the slight differences in slope

L3.Model.2b<-lmer(Math ~ ActiveTime.GM*ClassSize.SC
              +(1+ActiveTime.GM||School)
              +(1|School:Classroom.F),  
              data=MLM.Data.L3, REML=FALSE)
summary(L3.Model.2b)

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Math ~ ActiveTime.GM * ClassSize.SC + (1 + ActiveTime.GM || School) +  
##     (1 | School:Classroom.F)
##    Data: MLM.Data.L3
## 
##      AIC      BIC   logLik deviance df.resid 
##   3364.2   3399.0  -1674.1   3348.2      562 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8670 -0.6449 -0.0437  0.6589  2.5964 
## 
## Random effects:
##  Groups             Name          Variance Std.Dev.
##  School.Classroom.F (Intercept)   45.305   6.731   
##  School             ActiveTime.GM  7.781   2.789   
##  School.1           (Intercept)   60.824   7.799   
##  Residual                         16.492   4.061   
## Number of obs: 570, groups:  School:Classroom.F, 30; School, 3
## 
## Fixed effects:
##                            Estimate Std. Error       df t value Pr(>|t|)   
## (Intercept)                 44.3918     4.6772   3.0419   9.491  0.00234 **
## ActiveTime.GM               14.3027     1.7210   3.1071   8.311  0.00320 **
## ClassSize.SC                -0.0287     0.5723  27.0995  -0.050  0.96037   
## ActiveTime.GM:ClassSize.SC  -0.8289     0.2831 538.1476  -2.928  0.00356 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) AcT.GM ClS.SC
## ActiveTm.GM  0.000              
## ClassSiz.SC  0.031  0.002       
## AT.GM:CS.SC  0.002 -0.008  0.003

Level 1 Random Slopes at Level 3 and Level 2

• However, is that variation a function of level 3 (school) is it really a function of level 2 (the classroom)
• Lets say its both and see what happens

L3.Model.2c<-lmer(Math ~ ActiveTime.GM*ClassSize.SC
              +(1+ActiveTime.GM||School)
              +(1+ActiveTime.GM|School:Classroom.F),  
              data=MLM.Data.L3, REML=FALSE)
summary(L3.Model.2c)

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Math ~ ActiveTime.GM * ClassSize.SC + (1 + ActiveTime.GM || School) +  
##     (1 + ActiveTime.GM | School:Classroom.F)
##    Data: MLM.Data.L3
## 
##      AIC      BIC   logLik deviance df.resid 
##   3353.0   3396.5  -1666.5   3333.0      560 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.01413 -0.62593 -0.04133  0.65070  2.84509 
## 
## Random effects:
##  Groups             Name          Variance Std.Dev. Corr
##  School.Classroom.F (Intercept)   46.143   6.793        
##                     ActiveTime.GM 14.102   3.755    0.14
##  School             ActiveTime.GM  6.545   2.558        
##  School.1           (Intercept)   62.121   7.882        
##  Residual                         15.326   3.915        
## Number of obs: 570, groups:  School:Classroom.F, 30; School, 3
## 
## Fixed effects:
##                            Estimate Std. Error       df t value Pr(>|t|)   
## (Intercept)                44.35577    4.72608  3.03267   9.385  0.00245 **
## ActiveTime.GM              14.32683    1.73647  3.17701   8.251  0.00300 **
## ClassSize.SC               -0.02766    0.57719 27.22576  -0.048  0.96213   
## ActiveTime.GM:ClassSize.SC -0.77926    0.42308 29.96709  -1.842  0.07541 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) AcT.GM ClS.SC
## ActiveTm.GM 0.015               
## ClassSiz.SC 0.031  0.008        
## AT.GM:CS.SC 0.004  0.027  0.104

• Seems that maybe there is some variance at both we can capture.
• We lost the interaction: Seems that once we accounted for the random slopes at Level 2 we lost the effect. Why?
  • Remember each classroom had a different class size, adding variance to the slopes of ActiveTime.
  • However, the model was explaining that variance via the fixed effect interaction between ActiveTime and Classsize. When we treated ActiveTime as random that fixed variance went the random slope (in this particular case. Do not think this generalizes to all datasets)
• Lets make sure we need all those random effects (run a LLT)
• Also lets check a simplier fit of the random random structure

L3.Model.2c<-lmer(Math ~ ActiveTime.GM*ClassSize.SC
              +(1+ActiveTime.GM||School)
              +(1+ActiveTime.GM|School:Classroom.F),  
              data=MLM.Data.L3, REML=FALSE)
              
L3.Model.2c1<-lmer(Math ~ ActiveTime.GM*ClassSize.SC
              +(1|School)
              +(1+ActiveTime.GM|School:Classroom.F),  
              data=MLM.Data.L3, REML=FALSE) 

## Data: MLM.Data.L3
## Models:
## L3.Model.2: Math ~ ActiveTime.GM * ClassSize.SC + (1 | School) + (1 | School:Classroom.F)
## L3.Model.2b: Math ~ ActiveTime.GM * ClassSize.SC + (1 + ActiveTime.GM || School) + (1 | School:Classroom.F)
## L3.Model.2c: Math ~ ActiveTime.GM * ClassSize.SC + (1 + ActiveTime.GM || School) + (1 + ActiveTime.GM | School:Classroom.F)
##             npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## L3.Model.2     7 3379.9 3410.3 -1682.9   3365.9                         
## L3.Model.2b    8 3364.2 3399.0 -1674.1   3348.2 17.643  1  2.664e-05 ***
## L3.Model.2c   10 3353.0 3396.5 -1666.5   3333.0 15.205  2  0.0004993 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Data: MLM.Data.L3
## Models:
## L3.Model.2c1: Math ~ ActiveTime.GM * ClassSize.SC + (1 | School) + (1 + ActiveTime.GM | School:Classroom.F)
## L3.Model.2c: Math ~ ActiveTime.GM * ClassSize.SC + (1 + ActiveTime.GM || School) + (1 + ActiveTime.GM | School:Classroom.F)
##              npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
## L3.Model.2c1    9 3354.8 3393.9 -1668.4   3336.8                       
## L3.Model.2c    10 3353.0 3396.5 -1666.5   3333.0 3.7377  1     0.0532 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• Seems that this random structure, (1+ActiveTime.GM||School)+(1+ActiveTime.GM|School:Classroom.F) is marginally better than the simplier one, (1|School)+(1+ActiveTime.GM|School:Classroom.F). Which you use depends on your question. For now, lets use the more complex one for now

Plot Random Slopes

• Lets plot the random and fixed effects (per classroom) of our best fit model
• First, we will predict the model results

MLM.Data.L3$L3.Model.2c.FR<-predict(L3.Model.2c, newdata=MLM.Data.L3)

• Lets plot the raw data with the fitted of ActiveTime

theme_set(theme_bw(base_size = 7, base_family = "")) 

Active.Class.School <-ggplot(data = MLM.Data.L3, 
                             aes(x = ActiveTime, y=Math,group=Classroom.F))+
  facet_grid(~School)+
  coord_cartesian(ylim=c(10,80))+
  geom_point(aes(colour = Classroom.F))+
  geom_smooth(method = "lm", se = TRUE, aes(colour = Classroom.F))+
  xlab("Active Learning Time")+ylab("Math Score")+
  theme(legend.position = "none")


Active.Class.School.Fit2c <-ggplot(data = MLM.Data.L3, 
                             aes(x = ActiveTime, y=L3.Model.2c.FR, group=Classroom.F))+
  facet_grid(~School)+
  coord_cartesian(ylim=c(10,80))+
  geom_point(aes(colour = Classroom.F))+
  geom_smooth(method = "lm", se = TRUE, aes(colour = Classroom.F))+
  xlab("Active Learning Time")+ylab("Math Score - L3.Model.2c")+
  theme(legend.position = "none")

 grid.arrange(Active.Class.School, Active.Class.School.Fit2c, ncol=1, nrow=2)

Level 2 Random Slopes at Level 3

• Are all the random issues we need to deal with?
• What do we do with Classsize, we saw earlier it seemed to vary a function of the classroom?
• It purely a level 2 slope and cannot vary as a function of the classroom, but it can vary a function of the school
• Lets remove all that ActiveTime and deal only with class size to understand what is going to happen
• Back to the figure and lets used our grand mean centered variable

Size.School.2

• Now we can test the slopes of ClassSize relative to School, but looking at the figure above, we already know the random intercept will correlate with the slope perfectly.
• Note: Normally when you have many level 3 (lots of schools) this will not always be true (or easy to visualize)

L3.Model.2d<-lmer(Math ~ ActiveTime.GM*ClassSize.SC
              +(1+ClassSize.SC|School)
              +(1|School:Classroom.F),  
              data=MLM.Data.L3, REML=FALSE)
summary(L3.Model.2d, correlation=FALSE)

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Math ~ ActiveTime.GM * ClassSize.SC + (1 + ClassSize.SC | School) +  
##     (1 | School:Classroom.F)
##    Data: MLM.Data.L3
## 
##      AIC      BIC   logLik deviance df.resid 
##   3383.9   3423.0  -1682.9   3365.9      561 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.13981 -0.62403 -0.01151  0.62256  2.84109 
## 
## Random effects:
##  Groups             Name         Variance  Std.Dev. Corr
##  School:Classroom.F (Intercept)  4.503e+01 6.71051      
##  School             (Intercept)  6.043e+01 7.77384      
##                     ClassSize.SC 2.149e-04 0.01466  1.00
##  Residual                        1.725e+01 4.15279      
## Number of obs: 570, groups:  School:Classroom.F, 30; School, 3
## 
## Fixed effects:
##                             Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)                 44.41920    4.66234   3.03562   9.527  0.00234 ** 
## ActiveTime.GM               14.99664    0.60165 540.92818  24.926  < 2e-16 ***
## ClassSize.SC                -0.03996    0.57091  27.07138  -0.070  0.94471    
## ActiveTime.GM:ClassSize.SC  -0.82110    0.28928 540.84064  -2.838  0.00470 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

• Remove random slopes

L3.Model.2e<-lmer(Math ~ ActiveTime.GM*ClassSize.SC
              +(1+ClassSize.SC||School)
              +(1|School:Classroom.F),  
              data=MLM.Data.L3, REML=FALSE)
summary(L3.Model.2e, correlation=FALSE)

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Math ~ ActiveTime.GM * ClassSize.SC + (1 + ClassSize.SC || School) +  
##     (1 | School:Classroom.F)
##    Data: MLM.Data.L3
## 
##      AIC      BIC   logLik deviance df.resid 
##   3381.9   3416.7  -1682.9   3365.9      562 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.13965 -0.62420 -0.01136  0.62253  2.84136 
## 
## Random effects:
##  Groups             Name         Variance Std.Dev.
##  School.Classroom.F (Intercept)  45.03    6.711   
##  School             ClassSize.SC  0.00    0.000   
##  School.1           (Intercept)  60.36    7.769   
##  Residual                        17.25    4.153   
## Number of obs: 570, groups:  School:Classroom.F, 30; School, 3
## 
## Fixed effects:
##                            Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)                 44.4184     4.6597   3.0422   9.532  0.00231 ** 
## ActiveTime.GM               14.9967     0.6017 540.9271  24.926  < 2e-16 ***
## ClassSize.SC                -0.0368     0.5709  27.1003  -0.064  0.94908    
## ActiveTime.GM:ClassSize.SC  -0.8211     0.2893 540.8396  -2.839  0.00470 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

• Wait how did we get a zero random effect? Does this mean it has zero variance?
• We have few schools and only 3 slopes, and we are trying to test its fixed and random effects
• If we removed the fixed effects, the random slopes might appear

L3.Model.2f<-lmer(Math ~ ActiveTime.GM
              +(1+ClassSize.SC||School)
              +(1|School:Classroom.F),  
              data=MLM.Data.L3, REML=FALSE)
summary(L3.Model.2f)

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: 
## Math ~ ActiveTime.GM + (1 + ClassSize.SC || School) + (1 | School:Classroom.F)
##    Data: MLM.Data.L3
## 
##      AIC      BIC   logLik deviance df.resid 
##   3385.9   3412.0  -1686.9   3373.9      564 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3121 -0.6498 -0.0517  0.6337  3.0737 
## 
## Random effects:
##  Groups             Name         Variance Std.Dev.
##  School.Classroom.F (Intercept)  44.69    6.685   
##  School             ClassSize.SC  0.00    0.000   
##  School.1           (Intercept)  60.47    7.776   
##  Residual                        17.51    4.184   
## Number of obs: 570, groups:  School:Classroom.F, 30; School, 3
## 
## Fixed effects:
##               Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)     44.448      4.660   3.037   9.537  0.00233 ** 
## ActiveTime.GM   14.949      0.606 540.983  24.670  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## ActiveTm.GM 0.000 
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

• We got a number, but it is really small. Seems we should leave this level 2 random slope as fixed effect

Plot Final Results (Fixed effects only)

• Our best fit was L3.Model.2c
• Lets plot the interaction (but remember its tiny and not really significant)

library(effects)
# Set levels of support
SLevels<-c(mean(MLM.Data.L3$ClassSize.SC)-sd(MLM.Data.L3$ClassSize.SC),
                               mean(MLM.Data.L3$ClassSize.SC),
           mean(MLM.Data.L3$ClassSize.SC)+sd(MLM.Data.L3$ClassSize.SC))
#extract fixed effects
Results.BestFit<-Effect(c("ActiveTime.GM","ClassSize.SC"),L3.Model.2c,
     xlevels=list(ActiveTime.GM=seq(-.5,.5,.2), 
                  ClassSize.SC=SLevels))
#Convert to data frame for ggplot
Results.BestFit<-as.data.frame(Results.BestFit)
#Label Support for graphing
Results.BestFit$ClassSize.F<-factor(Results.BestFit$ClassSize.SC,
                        levels=SLevels,
                         labels=c("Small Class", "Mean Class", "Large Class"))

#Plot fixed effect
Final.Fixed.Plot.1 <-ggplot(data = Results.BestFit, 
                            aes(x = ActiveTime.GM, y =fit, group=ClassSize.F))+
  geom_line(size=2, aes(color=ClassSize.F,linetype=ClassSize.F))+
  coord_cartesian(xlim = c(-.5, .5),ylim = c(30, 60))+ 
  geom_ribbon(aes(ymin=fit-se, ymax=fit+se, group=ClassSize.F, fill=ClassSize.F),alpha=.2)+
  xlab("Proportion of Time Engaged in Active Learning \nCentered")+
  ylab("Math Score")+ 
  theme_bw()+
  theme(legend.position = "top", 
        legend.title=element_blank())
Final.Fixed.Plot.1