DESCRIPTION OF THE DATASETS

During the exercise sessions, we will use 2 datasets: the hip fracture dataset and the SII dataset. A description of both datasets is given here:

SII Dataset

The data come from the SII (Study of Instructional Improvement), which was carried out to assess the math achievement scores of first- and third-grade pupils in randomly selected classrooms from a national US sample of elementary schools. The dataset includes results for 1190 first-grade pupils sampled from 312 classrooms in 107 schools.
From the sample of elementary schools, all first-grade classrooms were considered (1 to 9). Within a classroom, 1 to 10 children were sampled to take a test on math. The following variables were measured: * schoolid: school’s ID number
* classid: classroom’s ID number
* childid: pupil’s ID number
* sex: an indicator variable for sex (0 for males, 1 for females)
* mathkind: pupil’s math score in the spring of the kindergarten year
* mathgain: pupil’s gain in the math achievement score from the spring of kindergarten to the spring of first grade
* ses: pupil’s socioeconomic status
* housepov: % of households in the neighborhood of the school below the poverty level
* yearstea: years of teacher’s experience in teaching in the first grade

Hip fracture dataset

Hip fractures are a very common and important health problem among elderly people. They are often the result of an unfortunate, seemingly innocent fall in the patient’s home environment. Less fluent movements combined with bad vision and impaired reaction skills make older people more at risk for such falls. The weaker bones and especially the presence of osteoporosis also causes the hip to break more easily. The hip surgery needed is a serious operation and can involve a lot of blood loss. A relatively long recovery period (8 to 12 days) in the hospital is therefore a real necessity.
The combination of these overwhelming events can also cause the mental status of the patients to worsen. When this deterioration is very severe, typically causing a loss of orientation in time and space, this phenomenon is called a delirium. In daily practice, the ‘Mini-Mental State Examination score’ (MMSE) is sometimes being used to measure the cognitive status of patients. It is the sum of 30 binary outcomes which result from correct (1) or incorrect (0) answer to specific questions (e.g. ‘Where are you now?’). Obviously, higher MMSE values indicate better cognitive functioning.

It is of interest to study how this MMSE evolution depends on certain characteristics of the patients. The data set used for this purpose contains the MMSE score for 60 patients, at 1, 3, 5, 8 and 12 days after surgery. For each patient the following variables are measured:

  • IDNR: id of the patient
  • NEURO: pre-operative neuro-status (1: neuro-psychiatric, 0: not neuro-psychiatric)
  • AGE: age of the patient
  • MMSE: Mini-Mental State Examination score
  • TIME: days after surgery: 1, 3, 5, 8 and 12.

1 PART ONE: EXPLORATORY & GRAPHICAL ANALYSIS

Describe in detail why in both datasets, one cannot make use of basic statistical models, like linear regression or ANOVA. What assumptions are violated? What special kind of data is the hip fracture dataset?

Classic ANOVA and regression assume all observations are independent. Here, we have related data - observations are not independent:

  • hip fracture : multiple measurements (of MMSE) over time in the same individual = longitudinal study
  • SII : here there is non-independence at two levels:
    • pupils within one class are not independent.
    • classes within schools are not independent

The purpose of multilevel (or mixed) models is to account for the relatedness between these observations

#Loading the required packages
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(ggplot2)
#> Warning: package 'ggplot2' was built under R version 4.5.3

1.1 SII Dataset

# Read in dataset
siiData<-read.table("SII.txt",header=T)

# Coerce schoolid, classid, childid, sex and ses to a factor
siiData$schoolid<-as.factor(siiData$schoolid)
siiData$classid<-as.factor(siiData$classid)
siiData$childid<-as.factor(siiData$childid)
siiData$sex<-as.factor(siiData$sex)

# Rename the factor levels of 'sex' for easier interpretation
levels(siiData$sex)<-c("male","female")

# Generating the Individual-vs.-Aggregate Dot Plot 
group_means <- siiData %>%
  group_by(schoolid,sex) %>%
  summarize(Mean_mathgain = mean(mathgain), .groups = "drop")

ggplot(group_means, aes(x = schoolid, y = Mean_mathgain, col = sex)) +
  geom_point(alpha = 0.7, size = 2.5) +
  labs(
    title = "Distribution of mathgain for males and females across schools",
    x = "School ID (Cluster)",
    y = "Mathgain",
    color = "Sex of student"
  ) +
  theme_minimal()

1.2 Hip fracture dataset

# Read in dataset
hipData<-read.table("hipfracture2.txt",header = T)

# coercion of categorical variables
hipData$IDNR<-as.factor(hipData$IDNR)
hipData$NEURO<-as.factor(hipData$NEURO)

# Rename the factor levels of 'NEURO' for easier interpretation
levels(hipData$NEURO)<-c("not neuro-psy","neuro-psy")

# one plot with all the patients, using a different color for neuro and non-Neuro patients
ggplot(hipData, aes(x = TIME, y = MMSE, group = IDNR, col=NEURO)) +
  geom_line(alpha = 0.6) +
  geom_point(alpha = 0.6) +
  labs(title = "Individual Trajectory Plot", y = "MMSE")
#> Warning: Removed 30 rows containing missing values or values outside the scale range
#> (`geom_line()`).
#> Warning: Removed 33 rows containing missing values or values outside the scale range
#> (`geom_point()`).


# two separate plots for the neuro and non-neuro persons
ggplot(hipData, aes(x = TIME, y = MMSE, group = IDNR, col = NEURO)) +
  geom_line(alpha = 0.6) +
  geom_point(alpha = 0.6) +
  # Split the plot into separate panels based on the NEURO variable
  facet_wrap(~NEURO) + 
  labs(
    title = "Individual Trajectory Plots by Neuro Status", 
    y = "MMSE",
    x = "Time"
  )
#> Warning: Removed 30 rows containing missing values or values outside the scale range
#> (`geom_line()`).
#> Removed 33 rows containing missing values or values outside the scale range
#> (`geom_point()`).

2 PART TWO: EXPLORATORY & GRAPHICAL ANALYSIS

Calculate the ICC and the Design Effect for both datasets. Is there an indication that a multilevel model might be necessary?

2.1 SII Dataset

library(lme4)
#> Loading required package: Matrix
library(performance)
#> Warning: package 'performance' was built under R version 4.5.3

fit_null_sii <- lme4::lmer(mathgain ~ 1 + (1 | schoolid), data = siiData)
summary(fit_null_sii)
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: mathgain ~ 1 + (1 | schoolid)
#>    Data: siiData
#> 
#> REML criterion at convergence: 11776.7
#> 
#> Scaled residuals: 
#>     Min      1Q  Median      3Q     Max 
#> -4.8127 -0.5903 -0.0186  0.5449  5.6632 
#> 
#> Random effects:
#>  Groups   Name        Variance Std.Dev.
#>  schoolid (Intercept)  109.2   10.45   
#>  Residual             1094.0   33.08   
#> Number of obs: 1190, groups:  schoolid, 107
#> 
#> Fixed effects:
#>             Estimate Std. Error t value
#> (Intercept)   57.583      1.442   39.95

# Calculate ICC
ICC_sii <- 109.2/(109.2+1094); ICC_sii
#> [1] 0.09075798
icc(fit_null_sii)
#> # Intraclass Correlation Coefficient
#> 
#>     Adjusted ICC: 0.091
#>   Unadjusted ICC: 0.091

# Calculate Average Cluster Size
avg_n_sii <- mean(table(siiData$schoolid))

# Calculate Design Effect
design_sii <- 1+(avg_n_sii-1)*ICC_sii ; design_sii
#> [1] 1.918606

The ICC is quite low (<0.1) and the design effect is below 2. So there is very few indication that a large proportion of the variance in the data is due to the variance between schools.

2.2 Hip fracture dataset

fit_null_hip <- lme4::lmer(MMSE ~ 1 + (1 | IDNR), data = hipData)
summary(fit_null_hip)
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: MMSE ~ 1 + (1 | IDNR)
#>    Data: hipData
#> 
#> REML criterion at convergence: 1483.8
#> 
#> Scaled residuals: 
#>     Min      1Q  Median      3Q     Max 
#> -4.2710 -0.3834  0.0528  0.4960  3.0246 
#> 
#> Random effects:
#>  Groups   Name        Variance Std.Dev.
#>  IDNR     (Intercept) 56.879   7.542   
#>  Residual              6.782   2.604   
#> Number of obs: 267, groups:  IDNR, 60
#> 
#> Fixed effects:
#>             Estimate Std. Error t value
#> (Intercept)  20.2470     0.9871   20.51
ICC_hip <- 56.879/(56.879+6.782); ICC_hip
#> [1] 0.893467
icc(fit_null_hip)
#> # Intraclass Correlation Coefficient
#> 
#>     Adjusted ICC: 0.893
#>   Unadjusted ICC: 0.893

# Calculate Average Cluster Size
avg_n_hip <- mean(table(hipData$IDNR))

# Calculate Design Effect
design_hip <- 1+(avg_n_hip-1)*ICC_hip ; design_hip
#> [1] 4.573868

The ICC is very high (>0.8) and the design effect is higher than 2. This indicates that a significant proportion of the variance in the data is due to the variance between patients. A multilevel model is necessary.

3 PART THREE: FITTING A MULTILEVEL MODEL + INTERPRETING THE MODEL SUMMARY

3.1 SII Dataset

The research question we want to answer is:
Is there a difference of math achievement score between males and females?

We will now look at a random intercept model. We still only consider the clustering on school-level and ignore the clustering on the class-level (this will be included in a later stage).

Write down the full model with assumptions

\(mathgain_{ij} = \alpha + a_i + \beta \cdot sex_j + \epsilon_{ij}\)

\(a_i \sim N(0, \sigma^2_{school})\), and \(\epsilon_{ij} \sim N(0, \sigma^2_{res})\) independent

Fit the model in R, discuss the within-cluster variability and between-cluster variability

model_sii_1 <- lmer(mathgain ~ sex + (1 | schoolid), data = siiData)

summary(model_sii_1)
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: mathgain ~ sex + (1 | schoolid)
#>    Data: siiData
#> 
#> REML criterion at convergence: 11772.3
#> 
#> Scaled residuals: 
#>     Min      1Q  Median      3Q     Max 
#> -4.8455 -0.6000 -0.0199  0.5417  5.6377 
#> 
#> Random effects:
#>  Groups   Name        Variance Std.Dev.
#>  schoolid (Intercept)  109.3   10.45   
#>  Residual             1093.8   33.07   
#> Number of obs: 1190, groups:  schoolid, 107
#> 
#> Fixed effects:
#>             Estimate Std. Error t value
#> (Intercept)   58.658      1.749  33.539
#> sexfemale     -2.133      1.966  -1.085
#> 
#> Correlation of Fixed Effects:
#>           (Intr)
#> sexfemale -0.566
  • The between-school variability is equal to \(109.3\).
  • The within-school variability is equal to \(1093.8\).
  • The total variability is equal to \(109.3+1093.8=1203.1\).
  • The ICC (variability due to differences between schools) is equal to \(109.3/1203.1=0.091\), meaning 9% of the variance is due to differences between schools, the rest of the variability is due to differences within a school.

What is the average mathgain for female students (across all schools)?

According to the model, the average mathgain for female students across all school is: \[ 58.658 - 2.133 = 56.525 \]

3.2 Hip fracture dataset

We want to study how MMSE develops over time, depending on certain characteristics of a patient. Fit a random intercept model to investigate this research question, taking into account the pre-operative neuro-status of the patient.

Write down the model with assumptions.

\(MMSE_{ij} = \alpha + a_i + \beta_1 \cdot NEURO_j + \beta_2 \cdot TIME_j + \beta_3 \cdot NEURO_jTIME_j + \epsilon_{ij}\)

\(a_i \sim N(0, \sigma^2_{patient})\), and \(\epsilon_{ij} \sim N(0, \sigma^2_{res})\) independent

Discuss the variability between and within patients.

model_hip_1 <- lmer(MMSE ~ TIME + NEURO + TIME:NEURO + (1 | IDNR), data = hipData)

summary(model_hip_1)
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: MMSE ~ TIME + NEURO + TIME:NEURO + (1 | IDNR)
#>    Data: hipData
#> 
#> REML criterion at convergence: 1452.8
#> 
#> Scaled residuals: 
#>     Min      1Q  Median      3Q     Max 
#> -4.1447 -0.4416  0.0153  0.4882  2.8102 
#> 
#> Random effects:
#>  Groups   Name        Variance Std.Dev.
#>  IDNR     (Intercept) 41.970   6.478   
#>  Residual              6.316   2.513   
#> Number of obs: 267, groups:  IDNR, 60
#> 
#> Fixed effects:
#>                     Estimate Std. Error t value
#> (Intercept)         22.01172    1.07702  20.438
#> TIME                 0.19753    0.05240   3.770
#> NEUROneuro-psy      -8.03853    1.86904  -4.301
#> TIME:NEUROneuro-psy -0.05981    0.09642  -0.620
#> 
#> Correlation of Fixed Effects:
#>             (Intr) TIME   NEURO-
#> TIME        -0.254              
#> NEUROnr-psy -0.576  0.146       
#> TIME:NEURO-  0.138 -0.543 -0.258
  • The between-patient variability \(\sigma^2_{patient}\) is equal to \(41.970\).
  • The within-patient variability \(\sigma^2_{res}\) is equal to \(6.316\).
  • The total variability \(\sigma^2_{patient}+\sigma^2_{res}\) is equal to \(41.970+6.316=48.286\).
  • The ICC (variability due to differences between patients) is equal to \(41.970/48.286=0.869\), meaning 86.9% of the variance is due to differences between patients, the rest of the variability is due to differences within a patient.

What is the marginal model (i.e. average model) that describes the MMSE over time for not neuro-psychiatric patients?

\(MMSE_{NotNeuro} = 22.01172 + 0.19753 \cdot TIME\)

4 PART FOUR: PARAMETER INFERENCE

4.1 SII Dataset

Based on the random intercept model that you fit in the previous part, give an answer to the question: Is there a difference of math achievement score between males and females?

library(lmerTest)
#> Warning: package 'lmerTest' was built under R version 4.5.3
#> 
#> Attaching package: 'lmerTest'
#> The following object is masked from 'package:lme4':
#> 
#>     lmer
#> The following object is masked from 'package:stats':
#> 
#>     step
library(pbkrtest)
#> Warning: package 'pbkrtest' was built under R version 4.5.3

model_sii_full <- lme4::lmer(mathgain ~ sex + (1 | schoolid), data = siiData, REML=FALSE)
model_sii_red <- lme4::lmer(mathgain ~  (1 | schoolid), data = siiData, REML=FALSE)

# LRT Test

anova(model_sii_red,model_sii_full)
#> Data: siiData
#> Models:
#> model_sii_red: mathgain ~ (1 | schoolid)
#> model_sii_full: mathgain ~ sex + (1 | schoolid)
#>                npar   AIC   BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)
#> model_sii_red     3 11785 11800 -5889.6     11779                     
#> model_sii_full    4 11786 11806 -5889.0     11778 1.1784  1     0.2777

# Wald test

model_sii_wald <- lmerTest::lmer(mathgain ~ sex + (1 | schoolid), data = siiData)

anova(model_sii_wald, ddf = "Kenward-Roger")
#> Type III Analysis of Variance Table with Kenward-Roger's method
#>     Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)
#> sex 1285.8  1285.8     1 1166.3  1.1755 0.2785

# Bootstrap

bootstrap_test <- PBmodcomp(model_sii_full, model_sii_red, nsim = 500)
bootstrap_test
#> Bootstrap test; time: 11.69 sec; samples: 500; extremes: 147;
#> large : mathgain ~ sex + (1 | schoolid)
#>          stat df p.value
#> LRT    1.1784  1  0.2777
#> PBtest 1.1784     0.2954

All the tests show a p-value larger than 0.05, indicating that there is no significant difference in math achievement scores between males and females.

4.2 Hip fracture dataset

Based on the random intercept model that you fit in the previous part, is there a different evolution of MMSE score in the two neuro-status groups?


model_hip_full <- lme4::lmer(MMSE ~ TIME + NEURO + TIME:NEURO + (1 | IDNR), data = hipData, REML=FALSE)
model_hip_red <- lme4::lmer(MMSE ~ TIME + NEURO + (1 | IDNR), data = hipData, REML=FALSE)

# LRT Test

anova(model_hip_red,model_hip_full)
#> Data: hipData
#> Models:
#> model_hip_red: MMSE ~ TIME + NEURO + (1 | IDNR)
#> model_hip_full: MMSE ~ TIME + NEURO + TIME:NEURO + (1 | IDNR)
#>                npar    AIC    BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)
#> model_hip_red     5 1460.4 1478.4 -725.21    1450.4                     
#> model_hip_full    6 1462.0 1483.5 -725.01    1450.0 0.3915  1     0.5315

# Wald test

model_hip_wald <- lmerTest::lmer(MMSE ~ TIME + NEURO + TIME:NEURO + (1 | IDNR), data = hipData)

anova(model_hip_wald, ddf = "Kenward-Roger")
#> Type III Analysis of Variance Table with Kenward-Roger's method
#>             Sum Sq Mean Sq NumDF   DenDF F value    Pr(>F)    
#> TIME        76.323  76.323     1 206.509 12.0849 0.0006196 ***
#> NEURO      116.823 116.823     1  66.554 18.4975 5.693e-05 ***
#> TIME:NEURO   2.429   2.429     1 206.509  0.3846 0.5358332    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Bootstrap

#bootstrap_test <- PBmodcomp(model_hip_full, model_hip_red, nsim = 500)
#bootstrap_test

Both the LRT test and the Wald test show a p-value above 0.05, indicating that there is no significantly different evolution of MMSE score in the two neuro-status groups. This means that we can remove the interaction term between TIME and NEURO from the model.

Note that if you perform parametric bootstrapping here, you get an error (that’s why the code was not run here). This is due to the fact that there are missing data in the dataset. If you want to perform parametric bootstrapping here, you have to limit your dataset to only include the complete cases.

Based on your conclusion from the previous question, fit a new appropriate model. Is there a significant MMSE evolution over time?


model_hip_full_2 <- lme4::lmer(MMSE ~ TIME + NEURO + (1 | IDNR), data = hipData, REML=FALSE)
model_hip_red_2 <- lme4::lmer(MMSE ~  NEURO + (1 | IDNR), data = hipData, REML=FALSE)

# LRT Test

anova(model_hip_red_2,model_hip_full_2)
#> Data: hipData
#> Models:
#> model_hip_red_2: MMSE ~ NEURO + (1 | IDNR)
#> model_hip_full_2: MMSE ~ TIME + NEURO + (1 | IDNR)
#>                  npar    AIC    BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)    
#> model_hip_red_2     4 1474.6 1488.9 -733.29    1466.6                         
#> model_hip_full_2    5 1460.4 1478.4 -725.21    1450.4 16.161  1  5.819e-05 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Wald test

model_hip_wald_2 <- lmerTest::lmer(MMSE ~ TIME + NEURO + (1 | IDNR), data = hipData)

anova(model_hip_wald_2, ddf = "Kenward-Roger")
#> Type III Analysis of Variance Table with Kenward-Roger's method
#>       Sum Sq Mean Sq NumDF   DenDF F value    Pr(>F)    
#> TIME  105.60  105.60     1 207.139  16.781 6.018e-05 ***
#> NEURO 133.87  133.87     1  58.068  21.274 2.238e-05 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Bootstrap

#bootstrap_test <- PBmodcomp(model_hip_full_2, model_hip_red_2, nsim = 500)
#bootstrap_test

Both the LRT test and the Wald test show a p-value below 0.05, indicating that there is a significant MMSE evolution over time.

Looking at the profiles from the hipfracture dataset (see figure produced in PART ONE of the exercises), it might be of interest to include a quadratic time effect.

hipData$TIME2 <- hipData$TIME* hipData$TIME

Fit the model. Is the quadratic time effect needed here?

model_hip_full_3 <- lme4::lmer(MMSE ~ TIME + NEURO + TIME2 + (1 | IDNR), data = hipData, REML=FALSE)
model_hip_red_3 <- lme4::lmer(MMSE ~  NEURO + TIME + (1 | IDNR), data = hipData, REML=FALSE)

# LRT Test

anova(model_hip_red_3,model_hip_full_3)
#> Data: hipData
#> Models:
#> model_hip_red_3: MMSE ~ NEURO + TIME + (1 | IDNR)
#> model_hip_full_3: MMSE ~ TIME + NEURO + TIME2 + (1 | IDNR)
#>                  npar    AIC    BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)
#> model_hip_red_3     5 1460.4 1478.4 -725.21    1450.4                     
#> model_hip_full_3    6 1460.3 1481.8 -724.16    1448.3 2.0937  1     0.1479

# Wald test

model_hip_wald_3 <- lmerTest::lmer(MMSE ~ TIME + NEURO + TIME2 + (1 | IDNR), data = hipData)

anova(model_hip_wald_3, ddf = "Kenward-Roger")
#> Type III Analysis of Variance Table with Kenward-Roger's method
#>        Sum Sq Mean Sq NumDF   DenDF F value    Pr(>F)    
#> TIME   38.825  38.825     1 205.200  6.2006   0.01356 *  
#> NEURO 133.636 133.636     1  58.068 21.3422 2.181e-05 ***
#> TIME2  13.036  13.036     1 205.338  2.0819   0.15058    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Both the LRT test and the Wald test show a p-value above 0.05, indicating that there is no need for a quadratic time effect.

Include the age of the patient as an extra covariate in the model. Is there a different effect of age in both neuro groups? 

model_hip_full_4 <- lme4::lmer(MMSE ~ TIME + NEURO + AGE + AGE:NEURO + (1 | IDNR), data = hipData, REML=FALSE)
model_hip_red_4 <- lme4::lmer(MMSE ~  NEURO + TIME + AGE + (1 | IDNR), data = hipData, REML=FALSE)

# LRT Test

anova(model_hip_red_4,model_hip_full_4)
#> Data: hipData
#> Models:
#> model_hip_red_4: MMSE ~ NEURO + TIME + AGE + (1 | IDNR)
#> model_hip_full_4: MMSE ~ TIME + NEURO + AGE + AGE:NEURO + (1 | IDNR)
#>                  npar    AIC    BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)
#> model_hip_red_4     6 1440.6 1462.1 -714.29    1428.6                     
#> model_hip_full_4    7 1442.5 1467.7 -714.27    1428.5 0.0285  1      0.866

# Wald test

model_hip_wald_4 <- lmerTest::lmer(MMSE ~ TIME + NEURO + AGE + AGE:NEURO + (1 | IDNR), data = hipData)

anova(model_hip_wald_4, ddf = "Kenward-Roger")
#> Type III Analysis of Variance Table with Kenward-Roger's method
#>            Sum Sq Mean Sq NumDF   DenDF F value    Pr(>F)    
#> TIME      101.546 101.546     1 207.506 16.1315 8.259e-05 ***
#> NEURO       0.479   0.479     1  55.947  0.0761 0.7837217    
#> AGE       106.330 106.330     1  55.997 16.8915 0.0001305 ***
#> NEURO:AGE   0.166   0.166     1  55.984  0.0263 0.8716731    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Both the LRT test and the Wald test show a p-value above 0.05, indicating that there is not a significantly different effect of age in both neuro groups. This means that we can drop the interaction term from the model.

Is the age effect needed here?

model_hip_full_5 <- lme4::lmer(MMSE ~ TIME + NEURO + AGE + (1 | IDNR), data = hipData, REML=FALSE)
model_hip_red_5 <- lme4::lmer(MMSE ~  NEURO + TIME + (1 | IDNR), data = hipData, REML=FALSE)

# LRT Test

anova(model_hip_red_5,model_hip_full_5)
#> Data: hipData
#> Models:
#> model_hip_red_5: MMSE ~ NEURO + TIME + (1 | IDNR)
#> model_hip_full_5: MMSE ~ TIME + NEURO + AGE + (1 | IDNR)
#>                  npar    AIC    BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)    
#> model_hip_red_5     5 1460.4 1478.4 -725.21    1450.4                         
#> model_hip_full_5    6 1440.6 1462.1 -714.29    1428.6 21.844  1  2.957e-06 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Wald test

model_hip_wald_5 <- lmerTest::lmer(MMSE ~ TIME + NEURO + AGE + (1 | IDNR), data = hipData)

anova(model_hip_wald_5, ddf = "Kenward-Roger")
#> Type III Analysis of Variance Table with Kenward-Roger's method
#>       Sum Sq Mean Sq NumDF   DenDF F value    Pr(>F)    
#> TIME  101.49  101.49     1 207.554  16.123 8.295e-05 ***
#> NEURO 162.88  162.88     1  57.117  25.875 4.223e-06 ***
#> AGE   157.65  157.65     1  57.191  25.043 5.674e-06 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Both the LRT test and the Wald test show a p-value below 0.05, indicating that there is significant effect of AGE on the MMSE value. The variable AGE should therefore be kept in the model.

5 PART FIVE: NESTED EFFECTS

5.1 SII Dataset

Give and fit the model, testing if there is a difference in mathgain between boys and girls, taking into account clustering on both school and class levels.

Write down this model (no other covariates have to be considered yet).

\(mathgain_{ijk}=\alpha +a_{1i} + a_{2i(j)} + \beta \cdot sex_k + \epsilon_{ijk}\)

  • with $a_{1i} = $ school effect
  • with \(a_{2i(j)}=\) class effect (within school)

Is there a difference in mathgain between boys and girls?


model_nested <- lme4::lmer(mathgain~ sex + (1|schoolid/classid), data=siiData, REML=F)
model_nested_red <- lme4::lmer(mathgain~ 1+ (1|schoolid/classid), data=siiData, REML=F)

anova(model_nested_red,model_nested )
#> Data: siiData
#> Models:
#> model_nested_red: mathgain ~ 1 + (1 | schoolid/classid)
#> model_nested: mathgain ~ sex + (1 | schoolid/classid)
#>                  npar   AIC   BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)
#> model_nested_red    4 11779 11800 -5885.7     11771                     
#> model_nested        5 11780 11805 -5884.9     11770 1.4574  1     0.2273

The p-value for the LRT is greater than 0.05, so there is no significant difference in mathgain between boys and girls.

Discuss the different variance components from the model in the previous exercise. * What fraction of the total variance is due to the variability between classes? * What fraction of the total variance is due to the variability between schools?


model_nested <- lme4::lmer(mathgain~ sex + (1|schoolid/classid), data=siiData)

summary(model_nested)
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: mathgain ~ sex + (1 | schoolid/classid)
#>    Data: siiData
#> 
#> REML criterion at convergence: 11764.1
#> 
#> Scaled residuals: 
#>     Min      1Q  Median      3Q     Max 
#> -4.6789 -0.5889 -0.0324  0.5396  5.6126 
#> 
#> Random effects:
#>  Groups           Name        Variance Std.Dev.
#>  classid:schoolid (Intercept)  101.04  10.052  
#>  schoolid         (Intercept)   76.96   8.773  
#>  Residual                     1026.86  32.045  
#> Number of obs: 1190, groups:  classid:schoolid, 312; schoolid, 107
#> 
#> Fixed effects:
#>             Estimate Std. Error t value
#> (Intercept)   58.615      1.747  33.547
#> sexfemale     -2.353      1.948  -1.208
#> 
#> Correlation of Fixed Effects:
#>           (Intr)
#> sexfemale -0.564
  • The total variance is equal to \(101.04+76.96+1026.86=1204.86\)
  • \(101.04/1204.86=0.084 \rightarrow 8.4 \%\) of the total variance is due to the variability between classes.
  • \(76.96/1204.86=0.064 \rightarrow 6.4\%\) of the total variance is due to the variability between schools.

Fit the same model, but consider the random effects for school and class levels to be crossed, rather than nested. Look at the variance components now. Do you notice a difference? Why is that?


model_crossed <- lme4::lmer(mathgain~ sex + (1|schoolid) + (1|classid), data=siiData)
summary(model_crossed)
#> Linear mixed model fit by REML ['lmerMod']
#> Formula: mathgain ~ sex + (1 | schoolid) + (1 | classid)
#>    Data: siiData
#> 
#> REML criterion at convergence: 11764.1
#> 
#> Scaled residuals: 
#>     Min      1Q  Median      3Q     Max 
#> -4.6789 -0.5889 -0.0324  0.5396  5.6126 
#> 
#> Random effects:
#>  Groups   Name        Variance Std.Dev.
#>  classid  (Intercept)  101.04  10.052  
#>  schoolid (Intercept)   76.96   8.773  
#>  Residual             1026.86  32.045  
#> Number of obs: 1190, groups:  classid, 312; schoolid, 107
#> 
#> Fixed effects:
#>             Estimate Std. Error t value
#> (Intercept)   58.615      1.747  33.547
#> sexfemale     -2.353      1.948  -1.208
#> 
#> Correlation of Fixed Effects:
#>           (Intr)
#> sexfemale -0.564

There is no difference in the estimated variance components compared to the previous model specification. This is because the variable classid is coded with globally unique identifiers for each class across the entire dataset, rather than with identifiers that repeat within schools (e.g., class 1, 2, 3 in school 1 and again class 1, 2, 3 in school 2). As a result, each class is uniquely associated with a single school, making the specification with nested random effects mathematically equivalent to the specification with crossed random effects for this dataset.

6 PART SIX: MODEL SELECTION

6.1 SII Dataset

Given the model from the previous part, with sex as a fixed effect and a random effect for school and class. Is it necessary to take into account the nested structure of classes within schools? Test whether the nested effects model performs better than a model that only includes a random effect for school, and a model that only includes a random effect for class.

model_nested <- lme4::lmer(mathgain~ sex + (1|schoolid/classid), data=siiData)
model_school <- lme4::lmer(mathgain~ sex + (1|schoolid), data=siiData)
model_class <- lme4::lmer(mathgain~ sex + (1|classid), data=siiData)

anova(model_school,model_nested)
#> refitting model(s) with ML (instead of REML)
#> Data: siiData
#> Models:
#> model_school: mathgain ~ sex + (1 | schoolid)
#> model_nested: mathgain ~ sex + (1 | schoolid/classid)
#>              npar   AIC   BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)   
#> model_school    4 11786 11806 -5889.0     11778                        
#> model_nested    5 11780 11805 -5884.9     11770 8.1813  1   0.004232 **
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(model_class,model_nested)
#> refitting model(s) with ML (instead of REML)
#> Data: siiData
#> Models:
#> model_class: mathgain ~ sex + (1 | classid)
#> model_nested: mathgain ~ sex + (1 | schoolid/classid)
#>              npar   AIC   BIC  logLik -2*log(L)  Chisq Df Pr(>Chisq)   
#> model_class     4 11785 11806 -5888.7     11777                        
#> model_nested    5 11780 11805 -5884.9     11770 7.4614  1   0.006303 **
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

AIC(model_school,model_class,model_nested)
#>              df      AIC
#> model_school  4 11780.30
#> model_class   4 11779.89
#> model_nested  5 11774.14
BIC(model_school,model_class,model_nested)
#>              df      BIC
#> model_school  4 11800.63
#> model_class   4 11800.22
#> model_nested  5 11799.55

The LRT, AIC and BIC all indicate that the fully nested model, with random effects for school and class performs the best.

Fit one model including all covariates as main effects and as an interaction with sex, and taking into account all levels of clustering. Discuss significance of the covariates. Perform stepwise backward model selection on the fixed effects to select the best model

model_fullynested <- lmer(mathgain~ sex*(ses+housepov+yearstea+mathkind) + (1|schoolid/classid), data=siiData, REML=F)
#> Warning: Some predictor variables are on very different scales: consider
#> rescaling
#> Warning: Some predictor variables are on very different scales: consider
#> rescaling
anova(model_fullynested)
#> Type III Analysis of Variance Table with Satterthwaite's method
#>              Sum Sq Mean Sq NumDF   DenDF  F value    Pr(>F)    
#> sex             231     231     1 1143.39   0.3141   0.57532    
#> ses           15624   15624     1 1176.79  21.1982 4.591e-06 ***
#> housepov       2305    2305     1  118.84   3.1268   0.07958 .  
#> yearstea        202     202     1  247.62   0.2744   0.60088    
#> mathkind     317401  317401     1 1173.82 430.6319 < 2.2e-16 ***
#> sex:ses           0       0     1 1136.29   0.0006   0.98007    
#> sex:housepov   2641    2641     1 1125.28   3.5837   0.05861 .  
#> sex:yearstea   2853    2853     1 1122.22   3.8704   0.04939 *  
#> sex:mathkind    665     665     1 1143.29   0.9028   0.34224    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

model_step2 <- lmer(mathgain~ sex*(housepov+yearstea+mathkind) + ses + (1|schoolid/classid), data=siiData, REML=F)
#> Warning: Some predictor variables are on very different scales: consider
#> rescaling
#> Warning: Some predictor variables are on very different scales: consider
#> rescaling
anova(model_step2)
#> Type III Analysis of Variance Table with Satterthwaite's method
#>              Sum Sq Mean Sq NumDF   DenDF  F value    Pr(>F)    
#> sex             238     238     1 1140.69   0.3231   0.56985    
#> housepov       2304    2304     1  118.71   3.1265   0.07960 .  
#> yearstea        203     203     1  247.36   0.2757   0.60003    
#> mathkind     317401  317401     1 1173.78 430.6343 < 2.2e-16 ***
#> ses           15626   15626     1 1176.76  21.2000 4.587e-06 ***
#> sex:housepov   2708    2708     1 1128.21   3.6747   0.05550 .  
#> sex:yearstea   2855    2855     1 1122.36   3.8739   0.04929 *  
#> sex:mathkind    697     697     1 1140.32   0.9457   0.33103    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

model_step3 <- lmer(mathgain~ sex*(housepov+yearstea)+ ses + mathkind + (1|schoolid/classid), data=siiData, REML=F)
anova(model_step3)
#> Type III Analysis of Variance Table with Satterthwaite's method
#>              Sum Sq Mean Sq NumDF   DenDF  F value    Pr(>F)    
#> sex            3038    3038     1 1143.60   4.1153   0.04273 *  
#> housepov       2321    2321     1  118.71   3.1435   0.07879 .  
#> yearstea        216     216     1  246.71   0.2931   0.58875    
#> ses           15491   15491     1 1176.63  20.9834 5.126e-06 ***
#> mathkind     318745  318745     1 1173.73 431.7661 < 2.2e-16 ***
#> sex:housepov   3177    3177     1 1127.79   4.3041   0.03825 *  
#> sex:yearstea   2842    2842     1 1122.78   3.8502   0.04999 *  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

model_step4 <- lmer(mathgain~ sex*(housepov)+ ses + mathkind + yearstea + (1|schoolid/classid), data=siiData, REML=F)
anova(model_step4)
#> Type III Analysis of Variance Table with Satterthwaite's method
#>              Sum Sq Mean Sq NumDF   DenDF  F value    Pr(>F)    
#> sex             667     667     1 1114.14   0.9015   0.34259    
#> housepov       2184    2184     1  118.73   2.9514   0.08841 .  
#> ses           15953   15953     1 1175.51  21.5611 3.813e-06 ***
#> mathkind     320784  320784     1 1172.86 433.5494 < 2.2e-16 ***
#> yearstea        231     231     1  246.76   0.3123   0.57678    
#> sex:housepov   2340    2340     1 1118.40   3.1620   0.07564 .  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

model_step5 <- lmer(mathgain~ sex*housepov + ses + mathkind  + (1|schoolid/classid), data=siiData, REML=F)
anova(model_step5)
#> Type III Analysis of Variance Table with Satterthwaite's method
#>              Sum Sq Mean Sq NumDF   DenDF  F value    Pr(>F)    
#> sex             669     669     1 1113.92   0.9048   0.34170    
#> housepov       2354    2354     1  116.91   3.1834   0.07698 .  
#> ses           16013   16013     1 1175.67  21.6527 3.638e-06 ***
#> mathkind     320733  320733     1 1173.15 433.6978 < 2.2e-16 ***
#> sex:housepov   2368    2368     1 1118.53   3.2020   0.07382 .  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

model_step6 <- lmer(mathgain~ sex + housepov + ses + mathkind  + (1|schoolid/classid), data=siiData, REML=F)
anova(model_step6)
#> Type III Analysis of Variance Table with Satterthwaite's method
#>          Sum Sq Mean Sq NumDF  DenDF  F value    Pr(>F)    
#> sex         459     459     1 1132.7   0.6192   0.43149    
#> housepov   2364    2364     1  116.7   3.1917   0.07661 .  
#> ses       16501   16501     1 1175.0  22.2751 2.646e-06 ***
#> mathkind 321639  321639     1 1172.0 434.1773 < 2.2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

model_step7 <- lmer(mathgain~ housepov + ses + mathkind  + (1|schoolid/classid), data=siiData, REML=F)
anova(model_step7)
#> Type III Analysis of Variance Table with Satterthwaite's method
#>          Sum Sq Mean Sq NumDF   DenDF  F value    Pr(>F)    
#> housepov   2394    2394     1  116.53   3.2282   0.07497 .  
#> ses       16638   16638     1 1175.56  22.4353 2.438e-06 ***
#> mathkind 322799  322799     1 1171.62 435.2651 < 2.2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

model_step8 <- lmer(mathgain~ ses + mathkind  + (1|schoolid/classid), data=siiData)
anova(model_step8)
#> Type III Analysis of Variance Table with Satterthwaite's method
#>          Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
#> ses       18190   18190     1 1164.0  24.466 8.675e-07 ***
#> mathkind 319907  319907     1 1164.9 430.269 < 2.2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The least significant variables that are excluded in each step are:
1) sex:ses
2) sex:mathkind
3) sex:yearsta
4) yearsta
5) sex:housepov
6) sex
7) housepov

The final model includes fixed effects for ses and mathkind and random effects for school and class.

7 PART SEVEN: MODEL ASSUMPTIONS

7.1 SII Dataset

Check the model assumptions of the final model you have chosen in the previous exercise.

model_final <- lmer(mathgain~ ses + mathkind  + (1|schoolid/classid), data=siiData)

# Residuals vs. Fitted Values Plot (Tests Equal Variance)
plot(model_final, 
     main = "Residuals vs. Fitted Values",
     xlab = "Fitted Values", ylab = "Residuals")

# Q-Q Plot of Residuals (Tests Normality of Residuals)
qqnorm(resid(model_final), main = "Normal Q-Q Plot: Residuals")
qqline(resid(model_final), col = "red")

# Extract Empirical Bayes estimates (random effects)
eb_estimates <- ranef(model_final)

# Histogram to check for Class-level outliers
ggplot(mapping = aes(x = eb_estimates$classid[[1]])) +
  geom_histogram(bins = 30, fill = "lightblue", color = "black") +
  labs(title = "Distribution of Class Random Effects", x = "Class Intercept Deviation")

# Histogram to check for School-level outliers
ggplot(mapping = aes(x = eb_estimates$schoolid[[1]])) +
  geom_histogram(bins = 15, fill = "lightgreen", color = "black") +
  labs(title = "Distribution of School Random Effects", x = "School Intercept Deviation")

# Scatterplot
mapping_df <- unique(siiData[, c("classid", "schoolid")])
mapping_df <- na.omit(mapping_df)

mapping_df$class_eb <- eb_estimates$classid[as.character(mapping_df$classid), 1]
mapping_df$school_eb      <- eb_estimates$schoolid[as.character(mapping_df$schoolid), 1]

# Plot the relationship between the levels
ggplot(mapping_df, aes(x = school_eb, y = class_eb)) +
  geom_point(alpha = 0.5, color = "darkblue") +
  geom_smooth(method = "lm", color = "red", se = FALSE) +
  labs(title = "Class Effects vs. School Effects", x = "School Random Effect", y = "Class Random Effect")
#> `geom_smooth()` using formula = 'y ~ x'
#> Warning: Removed 31 rows containing non-finite outside the scale range
#> (`stat_smooth()`).
#> Warning: Removed 31 rows containing missing values or values outside the scale range
#> (`geom_point()`).

  • Linearity - Residuals vs. Fitted values - There is no clear trend in the points, so linearity is fulfilled.
  • Homoscedasticity - Residuals vs. Fitted values - The point are evenly distributed around 0, so constant variance is fulfilled.
  • Normality - QQplot - The points follow the diagonal line nicely, so normality of the residuals is fulfilled.
  • Outliers - Histogram of EB estimates - There might be some outlying values to be investigated.
  • Relationship between random effects - Scatterplot of EB values - There is a clear trend between the random effects, indicating that there is strong correlation between the random effects.