Multilevel Model:

Random Intercept and Slope Models

Source: LEMMA online multilevel modelling course (University of Bristol)

Adding Level 1 Explanatory Variables

Random Intercept model
Random Slope model
Random Coefficient model

Adding Level 2 Explanatory Variables

Contextual effects
Cross-level interactions

loading libraries and data for this session

rm(list=ls())
setwd("D:/_Lab_Data/R/R_Labs")
library(lme4); library(lattice)

data1<- "HLM/Sample.txt"
mydata <- read.table(file = data1, sep = ",", header = TRUE)

head(mydata)

##   caseid schoolid score cohort90 female sclass schtype schurban schdenom
## 1     18        1     0       -6      1      2       0        1        0
## 2     17        1    10       -6      1      2       0        1        0
## 3     19        1     0       -6      1      4       0        1        0
## 4     20        1    40       -6      1      3       0        1        0
## 5     21        1    42       -6      1      2       0        1        0
## 6     13        1     4       -6      1      2       0        1        0

Random Intercept Models: Adding Level_1 Explanatory Variables

fit <- lmer(score ~ cohort90 + (1 | schoolid), data = mydata, REML = FALSE)
summary(fit)

## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: score ~ cohort90 + (1 | schoolid)
##    Data: mydata
## 
##       AIC       BIC    logLik  deviance  df.resid 
##  280921.6  280955.3 -140456.8  280913.6     33984 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1487 -0.7242  0.0363  0.7339  3.7097 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  schoolid (Intercept)  45.99    6.781  
##  Residual             219.29   14.808  
## Number of obs: 33988, groups:  schoolid, 508
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 30.55915    0.32254   94.74
## cohort90     1.21495    0.01553   78.24
## 
## Correlation of Fixed Effects:
##          (Intr)
## cohort90 -0.002

Plotting: Random Intercept

predscore <- fitted(fit)
data1<-cbind(predscore = predscore, cohort90 = mydata$cohort90, schoolid = mydata$schoolid)
data1<-unique(data1)
class(data1)

## [1] "matrix"

DF1<-data.frame(data1)
xyplot(predscore ~ cohort90, data = DF1, groups = schoolid, type = c("p", "l"), col = "blue")

***

Random Slope Models: Allowing for Different Slopes across Schools

fit <- lmer(score ~ cohort90 + (1 + cohort90 | schoolid), data = mydata, REML = FALSE)
summary(fit)

## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: score ~ cohort90 + (1 + cohort90 | schoolid)
##    Data: mydata
## 
##       AIC       BIC    logLik  deviance  df.resid 
##  280698.2  280748.8 -140343.1  280686.2     33982 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1008 -0.7202  0.0387  0.7264  3.5220 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  schoolid (Intercept)  42.8580  6.5466       
##           cohort90      0.1606  0.4007  -0.39
##  Residual             215.7393 14.6881       
## Number of obs: 33988, groups:  schoolid, 508
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 30.60963    0.31345   97.66
## cohort90     1.23390    0.02531   48.74
## 
## Correlation of Fixed Effects:
##          (Intr)
## cohort90 -0.266

Testing for random slopes

fita <- lmer(score ~ cohort90 + (1 | schoolid), data = mydata, REML = FALSE)
anova(fit, fita)

## Data: mydata
## Models:
## fita: score ~ cohort90 + (1 | schoolid)
## fit: score ~ cohort90 + (1 + cohort90 | schoolid)
##      Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)    
## fita  4 280922 280955 -140457   280914                            
## fit   6 280698 280749 -140343   280686 227.4      2  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Examining intercept and slope residuals for schools

VarCorr(fit)$schoolid

##             (Intercept)  cohort90
## (Intercept)   42.858014 -1.024175
## cohort90      -1.024175  0.160591
## attr(,"stddev")
## (Intercept)    cohort90 
##    6.546603    0.400738 
## attr(,"correlation")
##             (Intercept)   cohort90
## (Intercept)   1.0000000 -0.3903891
## cohort90     -0.3903891  1.0000000

myrandomeff <- ranef(fit, conVar = TRUE)
myrandomeff1<-myrandomeff[[1]]; class(myrandomeff1)

## [1] "data.frame"

head(myrandomeff1)

##   (Intercept)    cohort90
## 1   -5.904642  0.19801952
## 2    2.683568  0.04868368
## 3    1.727237  0.14995853
## 4   -7.802837  0.30472846
## 5    3.295991 -0.45398797
## 6   12.116169 -0.58015477

plot(myrandomeff1, xlab = "Intercept (u0j)", ylab = "Slope of cohort90 (u1j)")
abline(h = 0, col = "red")
abline(v = 0, col = "red")

Plotting: Random Intercept and Slope

predscore <- fitted(fit)
datapred <- cbind(predscore = predscore, cohort90 = mydata$cohort90, schoolid = mydata$schoolid)
datapred <- data.frame(unique(datapred))
xyplot(predscore ~ cohort90, data = datapred, groups = schoolid, type = c("p", "l"), col = "blue")

Between-school variance as a function of cohort

x <- c(-6:8)
y <- 42.859 - 2.048*x + 0.161*x^2
plot(x, y, type = "l", xlim = c(-6, 10), xlab = "Cohort90", ylab = "Predicted Score")

We can conclude that the mean attainment increased with cohort, and the variation in mean attainment among schools has decreased.

Random Coefficient Model: Adding a avariable for gender (dichotomous x )

(fit2a <- lmer(score ~ cohort90 + female + (1 + cohort90 | schoolid), data = mydata, REML = FALSE))

## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: score ~ cohort90 + female + (1 + cohort90 | schoolid)
##    Data: mydata
##       AIC       BIC    logLik  deviance  df.resid 
##  280558.1  280617.2 -140272.1  280544.1     33981 
## Random effects:
##  Groups   Name        Std.Dev. Corr 
##  schoolid (Intercept)  6.5249       
##           cohort90     0.4016  -0.39
##  Residual             14.6573       
## Number of obs: 33988, groups:  schoolid, 508
## Fixed Effects:
## (Intercept)     cohort90       female  
##      29.585        1.227        1.945

(fit2 <- lmer(score ~ cohort90 + female + (1 + cohort90 + female | schoolid), data = mydata, REML = FALSE))

## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: score ~ cohort90 + female + (1 + cohort90 + female | schoolid)
##    Data: mydata
##       AIC       BIC    logLik  deviance  df.resid 
##  280558.9  280643.2 -140269.4  280538.9     33978 
## Random effects:
##  Groups   Name        Std.Dev. Corr       
##  schoolid (Intercept)  6.3685             
##           cohort90     0.4021  -0.39      
##           female       1.1710   0.21 -0.11
##  Residual             14.6464             
## Number of obs: 33988, groups:  schoolid, 508
## Fixed Effects:
## (Intercept)     cohort90       female  
##      29.589        1.228        1.931

anova(fit2, fit2a)

## Data: mydata
## Models:
## fit2a: score ~ cohort90 + female + (1 + cohort90 | schoolid)
## fit2: score ~ cohort90 + female + (1 + cohort90 + female | schoolid)
##       Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## fit2a  7 280558 280617 -140272   280544                         
## fit2  10 280559 280643 -140269   280539 5.2362      3     0.1553

We cannot reject the null hypothesis and we conclude that the gender effect is the same for each school. We therefore revert to a model with a fixed coefficient for female.

Adding Level 2 Explanatory Variables

mydata$sclass1 <- mydata$sclass == 1
mydata$sclass2 <- mydata$sclass == 2
mydata$sclass4 <- mydata$sclass == 4
fit1a <- lmer(score ~ cohort90 + female + sclass1 + sclass2 + sclass4 + (1 + cohort90 | schoolid), data = mydata, REML = FALSE)
summary(fit1a)

## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: score ~ cohort90 + female + sclass1 + sclass2 + sclass4 + (1 +  
##     cohort90 | schoolid)
##    Data: mydata
## 
##       AIC       BIC    logLik  deviance  df.resid 
##  276712.3  276796.6 -138346.1  276692.3     33978 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.4875 -0.6997  0.0304  0.7071  3.8091 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  schoolid (Intercept)  22.5133  4.7448       
##           cohort90      0.1508  0.3884  -0.32
##  Residual             192.9457 13.8905       
## Number of obs: 33988, groups:  schoolid, 508
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 24.60987    0.27962   88.01
## cohort90     1.18283    0.02431   48.65
## female       1.96134    0.15428   12.71
## sclass1TRUE 11.08568    0.20639   53.71
## sclass2TRUE  5.87520    0.20405   28.79
## sclass4TRUE -3.73774    0.28453  -13.14
## 
## Correlation of Fixed Effects:
##             (Intr) chrt90 female s1TRUE s2TRUE
## cohort90    -0.150                            
## female      -0.296 -0.023                     
## sclass1TRUE -0.395 -0.054  0.008              
## sclass2TRUE -0.386 -0.020  0.009  0.539       
## sclass4TRUE -0.271 -0.036  0.013  0.358  0.357

Contextual effects

(fit2 <- update(fit1a, . ~ . + schtype))

## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: score ~ cohort90 + female + sclass1 + sclass2 + sclass4 + (1 +  
##     cohort90 | schoolid) + schtype
##    Data: mydata
##       AIC       BIC    logLik  deviance  df.resid 
##  276688.9  276781.7 -138333.4  276666.9     33977 
## Random effects:
##  Groups   Name        Std.Dev. Corr 
##  schoolid (Intercept)  4.5354       
##           cohort90     0.3849  -0.26
##  Residual             13.8922       
## Number of obs: 33988, groups:  schoolid, 508
## Fixed Effects:
## (Intercept)     cohort90       female  sclass1TRUE  sclass2TRUE  
##      24.279        1.184        1.964       11.031        5.856  
## sclass4TRUE      schtype  
##      -3.750        4.247

(fit3 <- update(fit2, . ~ . + schurban))

## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: score ~ cohort90 + female + sclass1 + sclass2 + sclass4 + (1 +  
##     cohort90 | schoolid) + schtype + schurban
##    Data: mydata
##       AIC       BIC    logLik  deviance  df.resid 
##  276681.9  276783.1 -138329.0  276657.9     33976 
## Random effects:
##  Groups   Name        Std.Dev. Corr 
##  schoolid (Intercept)  4.467        
##           cohort90     0.385   -0.26
##  Residual             13.893        
## Number of obs: 33988, groups:  schoolid, 508
## Fixed Effects:
## (Intercept)     cohort90       female  sclass1TRUE  sclass2TRUE  
##      25.260        1.182        1.967       11.033        5.847  
## sclass4TRUE      schtype     schurban  
##      -3.740        4.392       -1.437

(fit4 <- update(fit3, . ~ . + schdenom))

## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: score ~ cohort90 + female + sclass1 + sclass2 + sclass4 + (1 +  
##     cohort90 | schoolid) + schtype + schurban + schdenom
##    Data: mydata
##       AIC       BIC    logLik  deviance  df.resid 
##  276683.8  276793.5 -138328.9  276657.8     33975 
## Random effects:
##  Groups   Name        Std.Dev. Corr 
##  schoolid (Intercept)  4.468        
##           cohort90     0.385   -0.27
##  Residual             13.893        
## Number of obs: 33988, groups:  schoolid, 508
## Fixed Effects:
## (Intercept)     cohort90       female  sclass1TRUE  sclass2TRUE  
##     25.2484       1.1820       1.9666      11.0335       5.8473  
## sclass4TRUE      schtype     schurban     schdenom  
##     -3.7405       4.3981      -1.4622       0.1699

Cross-level interactions

mydata$cohort90Xschtype <- mydata$cohort90*mydata$schtype
(fit5 <- update(fit3, . ~ . + cohort90Xschtype))

## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: score ~ cohort90 + female + sclass1 + sclass2 + sclass4 + (1 +  
##     cohort90 | schoolid) + schtype + schurban + cohort90Xschtype
##    Data: mydata
##       AIC       BIC    logLik  deviance  df.resid 
##  276651.0  276760.7 -138312.5  276625.0     33975 
## Random effects:
##  Groups   Name        Std.Dev. Corr 
##  schoolid (Intercept)  4.5182       
##           cohort90     0.3715  -0.23
##  Residual             13.8871       
## Number of obs: 33988, groups:  schoolid, 508
## Fixed Effects:
##      (Intercept)          cohort90            female       sclass1TRUE  
##          25.1871            1.2135            1.9703           11.0194  
##      sclass2TRUE       sclass4TRUE           schtype          schurban  
##           5.8308           -3.7428            5.2908           -1.4038  
## cohort90Xschtype  
##          -0.5994

Overall Conclusion: The mean attainment is higher in independent schools than in state schools, but independent schools experienced a smaller increase in attainment with cohort.

Multilevel Model:

Random Intercept and Slope Models

Source: LEMMA online multilevel modelling course (University of Bristol)

Adding Level 1 Explanatory Variables

Random Intercept model

Random Slope model

Random Coefficient model

Adding Level 2 Explanatory Variables

Contextual effects

Cross-level interactions

loading libraries and data for this session

Random Intercept Models: Adding Level_1 Explanatory Variables

Plotting: Random Intercept

Random Slope Models: Allowing for Different Slopes across Schools

Testing for random slopes

Examining intercept and slope residuals for schools

Plotting: Random Intercept and Slope

Between-school variance as a function of cohort

We can conclude that the mean attainment increased with cohort, and the variation in mean attainment among schools has decreased.

Random Coefficient Model: Adding a avariable for gender (dichotomous x )

We cannot reject the null hypothesis and we conclude that the gender effect is the same for each school. We therefore revert to a model with a fixed coefficient for female.

Adding Level 2 Explanatory Variables

Contextual effects

Cross-level interactions

Overall Conclusion: The mean attainment is higher in independent schools than in state schools, but independent schools experienced a smaller increase in attainment with cohort.