Intro to MLM: Practical 4

Random slopes model

Following up from practical 2, we will continue to build up our multilevel model for school effects by adding a random slope for prior attainment.

---

Task 1: Fit a random slopes model

We have theoretical reasons to believe that the effect of KS2 varies across schools, since some schools do seem to make their pupils progress more than others. To test this hypothesis, we fit a random slopes model by running the following code:

m1_rs <- lmer(ks4stand ~ as.vector(scale(ks2stand))+
              (1 + as.vector(scale(ks2stand))|schoolID), data = valueadded, REML = F)

# Note that KS2 have been added like this "as.vector(scale(ks2stand))". This is done to prevent a convergence error. The function "scale" rescales KS2 scores to units of standard deviation, i.e. it "standardises".

summary(m1_rs)

## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: 
## ks4stand ~ as.vector(scale(ks2stand)) + (1 + as.vector(scale(ks2stand)) |  
##     schoolID)
##    Data: valueadded
## 
##      AIC      BIC   logLik deviance df.resid 
##  94083.5  94128.8 -47035.7  94071.5    14046 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.9474 -0.5315  0.0255  0.5764  4.5486 
## 
## Random effects:
##  Groups   Name                       Variance Std.Dev. Corr
##  schoolID (Intercept)                 8.964   2.994        
##           as.vector(scale(ks2stand))  1.295   1.138    0.50
##  Residual                            43.160   6.570        
## Number of obs: 14052, groups:  schoolID, 628
## 
## Fixed effects:
##                            Estimate Std. Error t value
## (Intercept)                -0.11737    0.13375  -0.877
## as.vector(scale(ks2stand))  6.42098    0.07797  82.353
## 
## Correlation of Fixed Effects:
##             (Intr)
## as.vct((2)) 0.282

Question

1.1. What is the correlation between the intercept and the slope?

Answer: The correlation is the second random effect reported in the model summary. This is 0.5

1.2. What sort of pattern is it?

Answer: A positive correlation between the intercept and the slope means there is a “fanning out” pattern. We will see this pattern very clearly in task 3.

---

Task 2: Compare the fit of the random slopes model to the random intercepts model

To compare the random intercepts models with the random slopes model, we first run the random intercepts model again

m1 <- lmer(ks4stand ~ ks2stand + (1|schoolID), data = valueadded, REML = F)
anova(m1, m1_rs)

## Data: valueadded
## Models:
## m1: ks4stand ~ ks2stand + (1 | schoolID)
## m1_rs: ks4stand ~ as.vector(scale(ks2stand)) + (1 + as.vector(scale(ks2stand)) | schoolID)
##       npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)    
## m1       4 94184 94215 -47088    94176                         
## m1_rs    6 94083 94129 -47036    94071 104.87  2  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

And then we compare using the anova function in R:

anova(m1, m1_rs)

## Data: valueadded
## Models:
## m1: ks4stand ~ ks2stand + (1 | schoolID)
## m1_rs: ks4stand ~ as.vector(scale(ks2stand)) + (1 + as.vector(scale(ks2stand)) | schoolID)
##       npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)    
## m1       4 94184 94215 -47088    94176                         
## m1_rs    6 94083 94129 -47036    94071 104.87  2  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Question:

2.1. Is the random slope specification a significant addition to the model?

Answer: The results indicate that, even though more complex (2 extra parameters), the random slopes model has a significantly better fit than the random intercepts model.

---

Task 3: Visualising results

To plot the school predicted lines, we can retrieve the fitted values from the model m1_rs, as such:

valueadded2 <- filter(valueadded, !is.na(ks4stand) & !is.na(ks2stand)) # this filter is necessary to avoid issues with missing values

valueadded2$pred_rs<-fitted(m1_rs)

Then we can plot the predicted school results by typing:

school_plot_rs<-ggplot(valueadded2, aes(x=ks2stand, y=pred_rs, group=factor(schoolID))) + 
  geom_smooth(method="lm", colour="black") +
  xlab("Standardised KS2 score") +
  ylab("Predicted KS4 score") +
  theme_bw()

school_plot_rs

Voilà! School predicted scores have varying slopes for the relationship between KS2 and GCSE scores. You can see that pupils in some schools make more progress than others on average and some make less.

Question

3.1. What sort of pattern is this?

Answer: Confirming what the model reported in task 1, we clearly see a “fanning out” pattern.

---

Task 4: Checking assumptions

We will check the assumptions for the random intercepts model (m1).

You can also plot the higher-level (school) residuals to check for normality

u0 <- ranef(m1, condVar = TRUE) # These are the residuals from model "m1"

hist(u0$schoolID$`(Intercept)`) # Histogram for school residuals

qqnorm(u0$schoolID$`(Intercept)`, pch = 1, frame = FALSE) # Q-Q plot for school residuals
qqline(u0$schoolID$`(Intercept)`, col = "steelblue", lwd = 2) # Q-Q plot for school residuals

You can also plot individual-level to check for normality

valueadded2$ind_resid <- residuals(m1)

hist(valueadded2$ind_resid)

qqnorm(valueadded2$ind_resid, pch = 1, frame = FALSE)
qqline(valueadded2$ind_resid, col = "steelblue", lwd = 2)

And finally, you can plot individual-level residuals against the predicted values:

valueadded2$pred <- fitted(m1) # retrieve predicted values from model m1

homoscedasticity <- ggplot(valueadded2, aes(y = ind_resid, x = pred)) + geom_point()

homoscedasticity

Questions

4.1. Is the normality assumption reasonable at the school level?

Answer: Residuals rarely ever look like a neat normal distribution. The right tail of the distribution seems longer than the left tail, which means that there may be variables that we still need to control for.

4.2. Is the normality assumption reasonable at the pupil level?

Answer: This time, the left tail of the distribution seems longer than the right tail, so again, that there may be variables that we still need to control for.

4.3. Is it reasonable to assume homoscedasticity?

Answer: The variability looks uneven. The right-hand side is narrower than the middle section. This is evidence of heteroscedasticity, but remember that this plot is generated from the residuals of the random intercepts model. The random slopes will suffer less from this issue because it allows for correlation between intercept and slope.