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?

1.2. What sort of pattern is it?

---

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?

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?

---

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

You can also plot individual-level to check for normality

valueadded2$ind_resid <- residuals(m1)

hist(valueadded2$ind_resid)

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?

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

4.3. Is it reasonable to assume homoscedasticity?