library(tidyverse)
library(knitr)
library(lme4)
library(broom.mixed)Lecture 13 AE: Multilevel models
Interpretation + Estimation
music <- read_csv("data/music-data.csv") |>
mutate(orchestra = factor(if_else(instrument == "orchestral instrument",
1, 0)),
large_ensemble = factor(if_else(perform_type == "Large Ensemble",
1, 0)))
Important
Start modeling building with exploratory data analysis of Level One and Level Two variables.
Model 1: Random intercepts model
This is also known as the unconditional means model.
#Model1 (Unconditional means model)
model1 <- lmer(na ~ 1 + (1 | id), data = music)
tidy(model1) |>
kable(digits = 3)| effect | group | term | estimate | std.error | statistic |
|---|---|---|---|---|---|
| fixed | NA | (Intercept) | 16.237 | 0.428 | 37.943 |
| ran_pars | id | sd__(Intercept) | 2.225 | NA | NA |
| ran_pars | Residual | sd__Observation | 4.739 | NA | NA |
Model 2: Add Level One covariate
Add large_ensemble as a Level One covariate. This is called a random slopes and intercepts model.
model2 <- lmer(na ~ large_ensemble + (large_ensemble | id),
data = music)
tidy(model2) |>
kable(digits = 3)| effect | group | term | estimate | std.error | statistic |
|---|---|---|---|---|---|
| fixed | NA | (Intercept) | 16.730 | 0.491 | 34.09 |
| fixed | NA | large_ensemble1 | -1.676 | 0.542 | -3.09 |
| ran_pars | id | sd__(Intercept) | 2.517 | NA | NA |
| ran_pars | id | cor__(Intercept).large_ensemble1 | -0.759 | NA | NA |
| ran_pars | id | sd__large_ensemble1 | 0.862 | NA | NA |
| ran_pars | Residual | sd__Observation | 4.666 | NA | NA |
Model 3: Add Level Two covariates
Add orchestra as a Level Two covariate.
model3 <- lmer(na ~ large_ensemble + orchestra + orchestra:large_ensemble +
(large_ensemble | id), data = music)
tidy(model3) |>
kable(digits = 3)| effect | group | term | estimate | std.error | statistic |
|---|---|---|---|---|---|
| fixed | NA | (Intercept) | 15.930 | 0.641 | 24.833 |
| fixed | NA | large_ensemble1 | -0.911 | 0.845 | -1.077 |
| fixed | NA | orchestra1 | 1.693 | 0.945 | 1.791 |
| fixed | NA | large_ensemble1:orchestra1 | -1.424 | 1.099 | -1.295 |
| ran_pars | id | sd__(Intercept) | 2.378 | NA | NA |
| ran_pars | id | cor__(Intercept).large_ensemble1 | -0.635 | NA | NA |
| ran_pars | id | sd__large_ensemble1 | 0.672 | NA | NA |
| ran_pars | Residual | sd__Observation | 4.670 | NA | NA |
Model 4: Add mpqnem as an additional Level Two covariate
Ex 1
Write the Level One and Level Two models.
Ex 2
Write the composite model.
model4 <- lmer(na ~ large_ensemble + orchestra + mpqnem +
orchestra:large_ensemble + mpqnem:large_ensemble +
(large_ensemble | id), data = music)
tidy(model4) |>
kable(digits = 3)| effect | group | term | estimate | std.error | statistic |
|---|---|---|---|---|---|
| fixed | NA | (Intercept) | 11.568 | 1.221 | 9.478 |
| fixed | NA | large_ensemble1 | -0.280 | 1.834 | -0.153 |
| fixed | NA | orchestra1 | 1.001 | 0.817 | 1.225 |
| fixed | NA | mpqnem | 0.148 | 0.038 | 3.893 |
| fixed | NA | large_ensemble1:orchestra1 | -0.949 | 1.106 | -0.858 |
| fixed | NA | large_ensemble1:mpqnem | -0.030 | 0.052 | -0.575 |
| ran_pars | id | sd__(Intercept) | 1.813 | NA | NA |
| ran_pars | id | cor__(Intercept).large_ensemble1 | -0.379 | NA | NA |
| ran_pars | id | sd__large_ensemble1 | 0.746 | NA | NA |
| ran_pars | Residual | sd__Observation | 4.670 | NA | NA |
Ex 3
Interpret the effect of mpqnem for solo and small ensembles.
Ex 4
Interpret the effect of mpqnem for large ensembles.
Compare models
You can get model summary statistics using the glance function from the broom.mixed package.
Ex 5
Calculate AIC and BIC for Models 1 - 4. Use these statistics to choose a final model.
# add code here