Multi-level Models (ELM Section 8.8)

author: Nancy Reid

date: March 4

The junior schools project data

library(faraway)
data(jsp)
head(jsp)
##   school class gender social raven id english math year
## 1      1     1   girl      9    23  1      72   23    0
## 2      1     1   girl      9    23  1      80   24    1
## 3      1     1   girl      9    23  1      39   23    2
## 4      1     1    boy      2    15  2       7   14    0
## 5      1     1    boy      2    15  2      17   11    1
## 6      1     1    boy      2    22  3      88   36    0
jspr = jsp[jsp$year == 2,]

Plotting

========================================================= The data =========================================================

names(jspr)
## [1] "school"  "class"   "gender"  "social"  "raven"   "id"      "english"
## [8] "math"    "year"

“Raven” is test score on entering (measure of ability?)

“Social” is measure of socio-economic status

50 schools; up to 4 classes in each school;

linear model

glin <- lm(math ~ raven * gender * social, data = jspr)
anova(glin)
## Analysis of Variance Table
## 
## Response: math
##                      Df  Sum Sq Mean Sq  F value    Pr(>F)    
## raven                 1 11480.5 11480.5 368.0625 < 2.2e-16 ***
## gender                1    44.1    44.1   1.4142  0.234668    
## social                8   779.4    97.4   3.1233  0.001725 ** 
## raven:gender          1     0.0     0.0   0.0004  0.984718    
## raven:social          8   582.6    72.8   2.3347  0.017460 *  
## gender:social         8   450.1    56.3   1.8038  0.072742 .  
## raven:gender:social   8   234.6    29.3   0.9400  0.482355    
## Residuals           917 28602.8    31.2                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

simpler model

summary(lm(math ~ raven + social, data = jspr))
## 
## Call:
## lm(formula = math ~ raven + social, data = jspr)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -20.8430  -3.2426   0.7726   3.7765  14.0825 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 17.02481    1.37451  12.386   <2e-16 ***
## raven        0.58040    0.03256  17.826   <2e-16 ***
## social2      0.04950    1.12938   0.044   0.9651    
## social3     -0.42893    1.19568  -0.359   0.7199    
## social4     -1.77452    1.05993  -1.674   0.0944 .  
## social5     -0.78228    1.18924  -0.658   0.5108    
## social6     -2.49373    1.26094  -1.978   0.0483 *  
## social7     -3.04851    1.29065  -2.362   0.0184 *  
## social8     -3.11746    1.77494  -1.756   0.0793 .  
## social9     -0.63278    1.12731  -0.561   0.5747    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.632 on 943 degrees of freedom
## Multiple R-squared:  0.2907, Adjusted R-squared:  0.2839 
## F-statistic: 42.93 on 9 and 943 DF,  p-value: < 2.2e-16

random effects model

library(lme4)
## Loading required package: Matrix
## Loading required package: Rcpp
schoolmod <- lmer(math ~ raven*social*gender + 
                    (1 | school) + (1 | school:class),
                  data = jspr)
anova(schoolmod)
## Analysis of Variance Table
##                     Df  Sum Sq Mean Sq  F value
## raven                1 10217.9 10217.9 374.4012
## social               8   615.7    77.0   2.8203
## gender               1    21.6    21.6   0.7917
## raven:social         8   577.5    72.2   2.6449
## raven:gender         1     2.5     2.5   0.0902
## social:gender        8   275.3    34.4   1.2608
## raven:social:gender  8   187.2    23.4   0.8572
jspr$craven <- jspr$raven - mean(jspr$raven)
schoolmod2 <- lmer(math ~ craven*social + (1|school) + (1 | school:class), data = jspr)

output

summary(schoolmod2)
## Linear mixed model fit by REML ['lmerMod']
## Formula: math ~ craven * social + (1 | school) + (1 | school:class)
##    Data: jspr
## 
## REML criterion at convergence: 5921.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.9872 -0.5250  0.1468  0.6237  2.6341 
## 
## Random effects:
##  Groups       Name        Variance Std.Dev.
##  school:class (Intercept)  1.177   1.085   
##  school       (Intercept)  3.148   1.774   
##  Residual                 27.141   5.210   
## Number of obs: 953, groups:  school:class, 90; school, 48
## 
## Fixed effects:
##                Estimate Std. Error t value
## (Intercept)    31.91123    1.19554  26.692
## craven          0.60584    0.18854   3.213
## social2         0.02361    1.27220   0.019
## social3        -0.63067    1.30888  -0.482
## social4        -1.96700    1.19707  -1.643
## social5        -1.35847    1.30023  -1.045
## social6        -2.26870    1.37374  -1.651
## social7        -2.55185    1.40555  -1.816
## social8        -3.39495    1.80137  -1.885
## social9        -0.83128    1.25347  -0.663
## craven:social2 -0.13206    0.20579  -0.642
## craven:social3 -0.22433    0.21888  -1.025
## craven:social4  0.03582    0.19488   0.184
## craven:social5 -0.15034    0.20889  -0.720
## craven:social6 -0.03860    0.23259  -0.166
## craven:social7  0.39826    0.23177   1.718
## craven:social8  0.25600    0.26154   0.979
## craven:social9 -0.08102    0.20550  -0.394
## 
## Correlation of Fixed Effects:
##             (Intr) craven socil2 socil3 socil4 socil5 socil6 socil7 socil8
## craven      -0.513                                                        
## social2     -0.879  0.479                                                 
## social3     -0.860  0.469  0.805                                          
## social4     -0.942  0.510  0.882  0.862                                   
## social5     -0.863  0.468  0.807  0.784  0.863                            
## social6     -0.826  0.446  0.769  0.754  0.829  0.754                     
## social7     -0.810  0.435  0.750  0.739  0.810  0.743  0.712              
## social8     -0.627  0.348  0.587  0.576  0.627  0.573  0.546  0.539       
## social9     -0.901  0.489  0.839  0.823  0.901  0.823  0.791  0.775  0.596
## craven:scl2  0.476 -0.910 -0.519 -0.435 -0.476 -0.436 -0.417 -0.404 -0.321
## craven:scl3  0.445 -0.860 -0.414 -0.417 -0.442 -0.405 -0.387 -0.375 -0.302
## craven:scl4  0.499 -0.966 -0.466 -0.456 -0.492 -0.456 -0.433 -0.424 -0.337
## craven:scl5  0.463 -0.899 -0.435 -0.424 -0.461 -0.396 -0.403 -0.393 -0.313
## craven:scl6  0.413 -0.810 -0.381 -0.377 -0.407 -0.374 -0.297 -0.344 -0.279
## craven:scl7  0.417 -0.813 -0.389 -0.381 -0.412 -0.377 -0.363 -0.295 -0.283
## craven:scl8  0.374 -0.718 -0.352 -0.346 -0.373 -0.343 -0.325 -0.324 -0.242
## craven:scl9  0.471 -0.916 -0.441 -0.432 -0.470 -0.430 -0.410 -0.399 -0.322
##             socil9 crvn:2 crvn:3 crvn:4 crvn:5 crvn:6 crvn:7 crvn:8
## craven                                                             
## social2                                                            
## social3                                                            
## social4                                                            
## social5                                                            
## social6                                                            
## social7                                                            
## social8                                                            
## social9                                                            
## craven:scl2 -0.455                                                 
## craven:scl3 -0.425  0.786                                          
## craven:scl4 -0.476  0.881  0.832                                   
## craven:scl5 -0.441  0.821  0.775  0.871                            
## craven:scl6 -0.394  0.736  0.700  0.783  0.728                     
## craven:scl7 -0.397  0.740  0.699  0.786  0.733  0.658              
## craven:scl8 -0.359  0.657  0.622  0.695  0.647  0.580  0.587       
## craven:scl9 -0.428  0.834  0.788  0.886  0.826  0.743  0.743  0.656

plots

par(mfrow=c(1,2))
qqnorm(resid(schoolmod2), main = "")
plot(fitted(schoolmod2), resid(schoolmod2),xlab = "fitted", ylab = "residuals")