Intro to MLM: Practical 1

For this practical, we will use data from a sample of the Teaching dataset of The Longitudinal Study of Young People in England, 2004-2006 available here.

We will be fitting a series of school value-added models, which are one of the most prominent examples of multilevel models applied in education research.

A traditional school value-added model is a model that attempts to isolate the “school effects” from the inherent variability/heterogeneity of the pupils. It attempts to ascertain what schools add to the progress of their pupils beyond what is expected of them, given their circumstances.

---

Typical workflow setup and data preparation

1.1. Define a working directory

You can use any directory in your computer. As in the example below:

setwd("C:/myfolder")

Remember to download the data to the folder you will define as working directory, as this makes matters easier.

1.2. Load packages

You can always load packages later on, but it is a good practice to load packages at the beginning of the session on the top of your script or R markdown file.

In this practical, we will use the packages haven, lme4 and ggplot2. Remember that if you haven’t installed them before, you need to do so before you call the library function:

install.packages("haven")
install.packages("lme4")
install.packages("ggplot2")
install.packages("dplyr")

Then you load them as such:

library(haven)
library(lme4)
library(ggplot2)
library(dplyr)

1.3. Read in data

You can download the data from the UKDS website. There are two SPSS datasets and we will be using the one named “lsype_15000_final_2011_05_04.sav”.

You can also download the file from the link below, but this link will stop working by the end of the session, so make sure to register and download the data from the UKDS if you want to continue working with this dataset.

To read this dataset into R, we need to use the package haven:

ype<-read_sav("https://github.com/A-mora/MLM_summer-school/raw/main/data/lsype_15000_final_2011_05_04.sav")

---

Select variables to use

In England and Wales, the Department for Education (DfE) publishes periodically the so-called performance tables, in which schools are assessed (and classified) according to the progress that their pupils make from one stage to another. Secondary schools are judged on the GCSE results of their pupils and the progress they made since the end of primary, when they sat the KS2 tests.

We will select the variables: pupilid, schoolID, ks2stand (KS2 scores), ks4stand (GCSE scores), gender, fsm and indschool.

For this, we need to use the function select of the dplyr package:

valueadded <- select(ype, pupilid, schoolID, 
                     ks2stand, ks4stand, gender, 
                     fsm) 

---

Task 1: Relationship between KS2 and GCSE

Plot the relationship between KS2 and GCSE scores

plot1 <- ggplot(valueadded, aes(x=ks2stand, y=ks4stand)) +
geom_point() + geom_smooth(aes(x=ks2stand, y=ks4stand), method = "lm")

plot1

Questions:

1.1 What can you observe in the plot?

Answer: The relationship between KS2 and GCSE scores is as expected: in general, the better the KS2 score, the better the GCSE score. Nevertheless, there is plenty of variability.

1.2 How correlated are KS2 and GCSE scores?

Answer: Obtain the correlation between scores. This indicates a moderate positive correlation as hinted by the previous plot.

cor(valueadded$ks2stand, valueadded$ks4stand, use="comp")

## [1] 0.6714863

---

Task 2: Variance components (empty) model

Fit an empty multilevel model of pupils within schools

library(lme4)

We will use the lmer functions, which stands for “linear mixed effects regression”. The basic syntax follows the conventions of most R packages running regression. You specify an outcome regressed ~ on variables. Each variable you add needs to be preceded by a + sign. You specify the data.

Note that we will use Maximum Likelihood (ML), not Restricted Maximum Likelihood (REML). REML is the default in lmer, so we need to specify REML = FALSE

Random effects are added within brackets after the fixed effects. 1 indicates that the constant is allowed to vary freely. The random effects are specified like this: (1|level2id). If you want to want random slopes, you specify (1 + variable|level2id)

m0 <- lmer(ks4stand ~ 1 + (1|schoolID), data = valueadded, REML = F)

summary(m0)

## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: ks4stand ~ 1 + (1 | schoolID)
##    Data: valueadded
## 
##      AIC      BIC   logLik deviance df.resid 
## 111792.8 111815.7 -55893.4 111786.8    15401 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7751 -0.6241  0.0163  0.6378  4.0338 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  schoolID (Intercept) 24.33    4.932   
##  Residual             75.85    8.709   
## Number of obs: 15404, groups:  schoolID, 657
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  -0.1125     0.2057  -0.547

This model is called “type 0” value-added model in the literature

Question

2.1. What is the proportion of variation that lies between schools in the empty model?

Answer: You should use the following equation:

\[VPC = \frac{\sigma_u^2}{\sigma_u^2 + \sigma_e^2}\] where \(\sigma_u^2\) is the variance at the school level and \(\sigma_e^2\) is the variance at the pupil level.

Replacing the values:

\[VPC = \frac{24.33}{24.33 + 75.85} = 24.28 \]

2.2. What does that value mean?

Answer: Approximately 24% of the total variance can be attributed to differences between schools.

Task 3: Is the MLM better than a single-level model?

To assess the statistical significance of the MLM, we need to compare it to a single-level model. To fit a single-level model, we type:

single <- lm(ks4stand ~ 1 , data = valueadded)

summary(single)

## 
## Call:
## lm(formula = ks4stand ~ 1, data = valueadded)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -23.026  -7.026  -0.026   6.974  38.974 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  0.02577    0.07999   0.322    0.747
## 
## Residual standard error: 9.928 on 15403 degrees of freedom
##   (366 observations deleted due to missingness)

These results are not very interesting in themselves, so we move on extract the loglikelihood and compare to the MLM.

L1 <- as.numeric(logLik(single)) # store this as numerical to re-use
L2 <- as.numeric(logLik(m0))

(D <- 2*(L2-L1)) # this is the deviance, we put it within brackets to print it immediately

## [1] 2641.369

pchisq(q=D, df=1, lower.tail=F) # To find the p-value

## [1] 0

Question

3.1. Is the addition of the school level statistically significant?

Answer: Yes it is! Even if you hadn’t done the loglikelihood test, the VPC of 24% is strong evidence of the relevance of school effects.

3.2. What does this mean in practice?

Answer: It means that it is extremely important to control for the rather large variability between schools in your models.

---