Multilevel data

• We can think of correlated data as a multilevel structure

  • Population elements are aggregated into groups
  • There are observational units and measurements at each level

Two types of effects

• Fixed effects: Effects that are of interest in the study
  • Can think of these as effects whose interpretations would be included in a write up of the study

Practice Questions

Radon is a carcinogen – a naturally occurring radioactive gas whose decay products are also radioactive – known to cause lung cancer in high concentrations. The EPA sampled more than 80,000 homes across the U.S. Each house came from a randomly selected county and measurements were made on each level of each home. Uranium measurements at the county level were included to improve the radon estimates.

1. What is the most basic level of observation (Level One)?
2. What are the group units (Level Two, Level Three, etc…)
3. What is the response variable?
4. Describe the within-group variation.
5. What are the fixed effects? What are the random effects?

Data: Music performance anxiety

The data musicdata.csv come from the Sadler and Miller (2010) study of the emotional state of musicians before performances. The dataset contains information collected from 37 undergraduate music majors who completed the Positive Affect Negative Affect Schedule (PANAS), an instrument produces a measure of anxiety (negative affect) and a measure of happiness (positive affect). This analysis will focus on negative affect as a measure of performance anxiety.

The primary variables we’ll use are

• na: negative affect score on PANAS (the response variable)
• perform_type: type of performance (Solo, Large Ensemble, Small Ensemble)
• instrument: type of intstrument (Voice, Orchestral, Piano)

Look at data

music <- read_csv("data/musicdata.csv")
music %>%
  filter(id %in% c(1, 43)) %>%
  group_by(id) %>%
  slice(1:3) %>%
  select(id, diary, perform_type, na, gender, instrument) %>%
  kable()

• What are the Level One observations? Level Two observations?

• What are the Level One variables? Level Two variables?

Univariate exploratory data analysis

Level One variables

Two ways to approach univariate EDA (visualizations and summary statistics) for Level One variables:

• Use individual observations (i.e., treat observations as independent)

• Use aggregated values for each Level Two observation

Unviariate EDA

p1 <- ggplot(data = music, aes(x = na)) + 
  geom_histogram(fill = "steelblue", color = "black", binwidth = 2) + 
  labs(x = "Individual negative affect", 
       title = "Negative affect scores")

p2 <- music %>%
  group_by(id) %>%
  summarise(mean_na = mean(na)) %>%
  ggplot(aes(x = mean_na)) + 
  geom_histogram(fill = "steelblue", color = "black", binwidth = 2) + 
  labs(x = "Mean negative affect", 
       title = "Mean negative affect scores")

p1 + p2

Bivariate exploratory data analysis

Goals

• Explore general association between the predictor and response variable
• Explore whether subjects at a given level of the predictor tend to have similar mean responses
• Explore whether variation in response differs at different levels of a predictor

Questions we want to answer

The goal is to understand variability in performance anxiety (na) based on performance-level and musician-level characteristics. Specifically:

• What is the association between performance type (large ensemble or not) and performance anxiety? Does the association differ based on instrument type (orchestral or not)?

music <- music %>%
  mutate(orchestra = if_else(instrument == "orchestral instrument", 1, 0), 
         large_ensemble = if_else(perform_type == "Large Ensemble", 1,0))

ols <- lm(na ~ orchestra + large_ensemble + orchestra * large_ensemble, 
          data = music)
tidy(ols) %>% kable(digits = 3)

Other modeling approaches

1️⃣ Condense each musician’s set of responses into a single outcome (e.g., mean max, last observation, etc.) and fit a linear model on these condensed observations

• Leaves few observations (37) to fit the model
• Ignoring a lot of information in the multiple observations for each musician

Two-level model

We’ll start with the Level One model to understand the association between performance type and performance anxiety for the \(i^{th}\) musician.

\[na_{ij} = a_i + b_i ~ LargeEnsemble_{ij} + \epsilon_i, \hspace{5mm} \epsilon_{ij} \sim N(0,\sigma^2)\]

Example Level One model

Below is partial data for observation #22

music %>%
  filter(id == 22) %>%
  select(id, diary, perform_type, instrument, na) %>%
  slice(1:3, 13:15) %>%
  kable()

Level One model

music %>%
  filter(id == 22) %>%
  lm(na ~ large_ensemble, data = .) %>%
  tidy() %>%
  kable(digits = 3)

Level Two Model

The slope and intercept for the \(i^{th}\) musician can be modeled as

\[\begin{aligned}&a_i = \alpha_0 + \alpha_1 ~ Orchestra_i + u_i \\ &b_i = \beta_0 + \beta_1 ~ Orchestra_i + v_i\end{aligned}\]

Estimated coefficients by instrument

musicians <- music %>%
  distinct(id, orchestra) %>%
  bind_cols(model_stats)

p1 <- ggplot(data = musicians, aes(x = intercepts, y = factor(orchestra))) + 
  geom_boxplot(fill = "steelblue", color = "black") + 
  labs(x = "Fitted intercepts", 
       y = "Orchestra")

p2 <- ggplot(data = musicians, aes(x = slopes, y = factor(orchestra))) + 
  geom_boxplot(fill = "steelblue", color = "black") + 
  labs(x = "Fitted slopes", 
       y = "Orchestra")

p1 / p2

Level Two model

a <- lm(intercepts ~ orchestra, data = musicians) 
tidy(a) %>%
  kable(digits = 3)

b <- lm(slopes ~ orchestra, data = musicians) 
tidy(b) %>%
  kable(digits = 3)

Writing out the models

Level One

\[\hat{na}_{ij} = \hat{a}_i + \hat{b}_i ~ LargeEnsemble_{ij}\]

for each musician.

Level Two

\[\begin{aligned}&\hat{a}_i = 16.283 + 1.441 ~ Orchestra_i \\ &\hat{b}_i = -0.771 - 1.406 ~ Orchestra_i\end{aligned}\]

Composite model

\[\begin{aligned}\hat{na}_i &= 16.283 + 1.441 ~ Orchestra_i - 0.771 ~ LargeEnsemble_{ij} \\ &- 1.406 ~ Orchestra:LargeEnsemble_{ij}\end{aligned}\]

Disadvantages to this approach

⚠️ Weighs each musician the same regardless of number of diary entries

⚠️ Drops subjects who have missing values for slope (7 individuals who didn’t play a large ensemble performance)

⚠️ Does not share strength effectively across individuals

Framework

Let \(Y_{ij}\) be the performance anxiety for the \(i^{th}\) musician before performance \(j\).

Level One

\[Y_{ij} = a_i + b_i ~ LargeEnsemble + \epsilon_{ij}\]

lmer Syntax

lmer(response~fixed effect variables + (random effect variables|group variable),data = your_data)

• we determine fixed and random effects based on our composite model
• terms with Greek letters will be fixed effects
• terms with random components, usually \(u\), \(v\), or \(w\) will be random effects

Composite model

Plug in the equations for \(a_i\) and \(b_i\) to get the composite model

\[\begin{aligned}Y_{ij} &= (\alpha_0 + \alpha_1 ~ Orchestra_i + \beta_0 ~ LargeEnsemble_{ij} \\ &+ \beta_1 ~ Orchestra_i:LargeEnsemble_{ij})\\ &+ (u_i + v_i ~ LargeEnsemble_{ij} + \epsilon_{ij})\end{aligned}\]

• The fixed effects to estimate are \(\alpha_0, \alpha_1, \beta_0, \beta_1\)
• The error terms are \(u_i, v_i, \epsilon_{ij}\)

Notation

• Greek letters denote the fixed effect model parameters to be estimated

  • e.g., \(\alpha_0, \alpha_1, \beta_0, \beta_1\)

• Roman letters denote the preliminary fixed effects at lower levels that will not be estimated directly.

  • e.g. \(a_i, b_i\)

• \(\sigma\) and \(\rho\) denote variance components that will be estimated

• \(\epsilon_{ij}, u_i, v_i\) denote error terms

Error terms

• We generally assume that the error terms are normally distributed, e.g. error associated with each performance of a given musician is \(\epsilon_{ij} \sim N(0, \sigma^2)\)

• For the Level Two models, the errors are

  • \(u_i\): deviation of musician \(i\) from the mean performance anxiety before solos and small ensembles after accounting for the instrument
  • \(v_i\): deviance of musician \(i\) from the mean difference in performance anxiety between large ensembles and other performance types after accounting for instrument

• Need to account for fact that \(u_i\) and \(v_i\) are correlated for the \(i^{th}\) musician

Intercept and Slope

musicians %>%
  filter(!is.na(slopes)) %>%
ggplot(aes(x = intercepts, y = slopes)) + 
  geom_point() + 
  geom_smooth(method = "lm") + 
  labs(x = "Fitted intercepts",
       y = "Fitted slopes", 
       title = "Fitted slopes and intercepts", 
       subtitle = paste0("r = ", round(cor(musicians %>% filter(!is.na(slopes)) %>% pull(intercepts), musicians %>% filter(!is.na(slopes)) %>% pull(slopes)),3)))

Distribution of Level Two errors

Use a multivariate normal distribution for the Level Two error terms

\[\left[ \begin{array}{c} u_{i} \\ v_{i} \end{array} \right] \sim N \left( \left[ \begin{array}{c} 0 \\ 0 \end{array} \right], \left[ \begin{array}{cc} \sigma_{u}^{2} & \rho_{uv}\sigma_{u}\sigma_v \\ \rho_{uv}\sigma_{u}\sigma_v & \sigma_{v}^{2} \end{array} \right] \right)\]

where \(\sigma^2_u\) and \(\sigma^2_v\) are the variance of \(u_i\)’s and \(v_i\)’s respectively, and \(\sigma_{uv} = \rho_{uv}\sigma_u\sigma_v\) is covariance between \(u_i\) and \(v_i\)

• What does it mean for \(\rho_{uv} > 0\)?

• What does it mean for \(\rho_{uv} < 0\)?

• Being able to write out these mammoth variance-covariance matrices is less important than recognizing the number of variance components that must be estimated by our intended model.

Visualizing multivariate normal distribution

Recreated from Figure 8.12

Fit the model (1/2)

Fit multilevel model using the lmer function from the lme4 package. Display results using hte tidy() function from the broom.mixed package.

library(lme4)
library(broom.mixed)
music_model <- lmer(na ~ orchestra + large_ensemble + 
                      orchestra:large_ensemble + (large_ensemble|id), 
                    REML = TRUE, data = music)

tidy(music_model) %>% kable(digits = 3)

Fit the model in R (2/2)

na ~ orchestra + large_ensemble + orchestra:large_ensemble: Represents the fixed effects + intercept

Tidy output

Display results using the tidy function from the broom.mixed package.

library(broom.mixed)
tidy(music_model) 

• Get fixed effects only

Model output: Fixed effects

Model output: Random effects

Interpret the effects

• Split into 3 groups.

• Each group will write the interpretation for one main effect and one variance component. The terms on each slide do not necessarily have a direct correspondence to each other.

• One person write the group’s interpretations on the board.

Comparing fixed effects estimates

Below are the estimated coefficients and standard errors for four modeling approaches.

Pseudo R-squared Values

\[\textrm{Pseudo }R^2_{L1} = \frac{\hat{\sigma}^{2}(\textrm{Model A})-\hat{\sigma}^{2}(\textrm{Model B})}{\hat{\sigma}^{2}(\textrm{Model A})}\] - Pseudo R-squared values are not universally reliable as measures of model performance. - Because of the complexity of estimating fixed effects and variance components at various levels of a multilevel model, it is not unusual to encounter situations in which covariates in a Level Two equation for, say, the intercept remain constant (while other aspects of the model change), yet the associated pseudo R-squared values differ or are negative. For this reason, pseudo R-squared values in multilevel models should be interpreted cautiously.

ML and REML

Maximum Likelihood (ML) and Restricted (Residual) Maximum Likelihood (REML) are the two most common methods for estimating the fixed effects and variance components

ML and REML

Restricted Maximum Likelihood (REML)

• Estimate the variance components using the sample residuals not the sample data
• It is conditional on the fixed effects, so it accounts for uncertainty in fixed effects estimates. This results in unbiased estimates of variance components.

ML or REML?

• Research has not determined one method absolutely superior to the other

• REML (REML = TRUE; default in lmer) is preferable when

  • the number of parameters is large or primary, or
  • primary objective is to obtain estimates of the model parameters

• ML (REML = FALSE) must be used if you want to compare nested fixed effects models using a likelihood ratio test (e.g., a drop-in-deviance test).

  • For REML, the goodness-of-fit and likelihood ratio tests can only be used to draw conclusions about variance components

Comparing ML and REML

Strategy for building multilevel models

• Conduct exploratory data analysis for Level One and Level Two variables

• Fit model with no covariates to assess variability at each level

• Create Level One models. Start with a single term, then add terms as needed.

• Create Level Two models. Start with a single term, then add terms as needed. Start with equation for intercept term.

• Begin with the full set of variance components, then remove variance terms as needed.

Acknowledgements

The content in the slides is from

• BMLR: Chapter 8 - Introduction to Multilevel Models
  • Initial versions of the slides are by Dr. Maria Tackett, Duke University
• Singer, J. D. & Willett, J. B. (2003). Applied longitudinal data analysis: Modeling change and event occurrence. Oxford university press.
  
• Sadler, Michael E., and Christopher J. Miller. 2010. “Performance Anxiety: A Longitudinal Study of the Roles of Personality and Experience in Musicians.” Social Psychological and Personality Science 1 (3): 280–87. http://dx.doi.org/10.1177/1948550610370492.
• Initial versions of the slides are by Dr. Maria Tackett, Duke University