Multilevel models

# Multilevel models
## Part 2
### 03.02.22

---

## [Click here for PDF of slides](16-multilevel-models-pt2.pdf)

---

## Announcements

- Tomorrow's lab: article evaluation for mini-project 02

- Looking ahead to the individual final project...
  - Choose a dataset with **correlated observations** and analyze it 
  - More details after spring break
  
---

## Learning goals

- Conduct univariate and bivariate EDA for multilevel models 
 
- Write multilevel model , including assumptions about variance components, in by-level and composite forms

- Interpret the model parameters, fixed effects, and variance components

---

## Data: Music performance anxiety

.midi[The data [`musicdata.csv`](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 instrument (Voice, Orchestral, Piano)

---

## Look at data

.small[
<table>
 <thead>
  <tr>
   <th style="text-align:right;"> id </th>
   <th style="text-align:right;"> diary </th>
   <th style="text-align:left;"> perform_type </th>
   <th style="text-align:right;"> na </th>
   <th style="text-align:left;"> gender </th>
   <th style="text-align:left;"> instrument </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:left;"> Solo </td>
   <td style="text-align:right;"> 11 </td>
   <td style="text-align:left;"> Female </td>
   <td style="text-align:left;"> voice </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 2 </td>
   <td style="text-align:left;"> Large Ensemble </td>
   <td style="text-align:right;"> 19 </td>
   <td style="text-align:left;"> Female </td>
   <td style="text-align:left;"> voice </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 3 </td>
   <td style="text-align:left;"> Large Ensemble </td>
   <td style="text-align:right;"> 14 </td>
   <td style="text-align:left;"> Female </td>
   <td style="text-align:left;"> voice </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 43 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:left;"> Solo </td>
   <td style="text-align:right;"> 19 </td>
   <td style="text-align:left;"> Female </td>
   <td style="text-align:left;"> voice </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 43 </td>
   <td style="text-align:right;"> 2 </td>
   <td style="text-align:left;"> Solo </td>
   <td style="text-align:right;"> 13 </td>
   <td style="text-align:left;"> Female </td>
   <td style="text-align:left;"> voice </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 43 </td>
   <td style="text-align:right;"> 3 </td>
   <td style="text-align:left;"> Small Ensemble </td>
   <td style="text-align:right;"> 19 </td>
   <td style="text-align:left;"> Female </td>
   <td style="text-align:left;"> voice </td>
  </tr>
</tbody>
</table>
]

Draw the data structure, and add the Level One and Level Two observations and 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

**Level Two variables**

- Use a data set that contains one row per Level Two observation

---

## Unviariate EDA

---

## 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

There are two ways to visualize these associations:

- One plot of individual observations (i.e., treat observations as independent)

- Separate plots of responses vs. predictor for each Level Two observation (lattice plots)

---

---

## Fitting the model

---

## 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)?

We will fit the model in two parts:

1️⃣ Fit a separate model for each musician understand the association between performance type and anxiety (Level One models).

2️⃣  Then fit a system of models to predict the fitted coefficients in the Level One models based on instrument type (Level Two models).

---

[Click here](https://forms.gle/evDkJ64T5wrFYvCfA) to submit your response.
]

---

## 1️⃣ Level One models

We'll start with a Level One model to understand the association between performance type and performance anxiety for the `$i^{th}$` musician.

<br>

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

For now, we will estimate `$a_i$` and `$b_i$` using least-squares regression.

---

## Example Level One model

Below is partial data for observation #22

<table>
 <thead>
  <tr>
   <th style="text-align:right;"> id </th>
   <th style="text-align:right;"> diary </th>
   <th style="text-align:left;"> perform_type </th>
   <th style="text-align:left;"> instrument </th>
   <th style="text-align:right;"> na </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:right;"> 22 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:left;"> Solo </td>
   <td style="text-align:left;"> orchestral instrument </td>
   <td style="text-align:right;"> 24 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 22 </td>
   <td style="text-align:right;"> 2 </td>
   <td style="text-align:left;"> Large Ensemble </td>
   <td style="text-align:left;"> orchestral instrument </td>
   <td style="text-align:right;"> 21 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 22 </td>
   <td style="text-align:right;"> 3 </td>
   <td style="text-align:left;"> Large Ensemble </td>
   <td style="text-align:left;"> orchestral instrument </td>
   <td style="text-align:right;"> 14 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 22 </td>
   <td style="text-align:right;"> 13 </td>
   <td style="text-align:left;"> Large Ensemble </td>
   <td style="text-align:left;"> orchestral instrument </td>
   <td style="text-align:right;"> 12 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 22 </td>
   <td style="text-align:right;"> 14 </td>
   <td style="text-align:left;"> Large Ensemble </td>
   <td style="text-align:left;"> orchestral instrument </td>
   <td style="text-align:right;"> 19 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 22 </td>
   <td style="text-align:right;"> 15 </td>
   <td style="text-align:left;"> Solo </td>
   <td style="text-align:left;"> orchestral instrument </td>
   <td style="text-align:right;"> 25 </td>
  </tr>
</tbody>
</table>

---

## Level One model

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

<table>
 <thead>
  <tr>
   <th style="text-align:left;"> term </th>
   <th style="text-align:right;"> estimate </th>
   <th style="text-align:right;"> std.error </th>
   <th style="text-align:right;"> statistic </th>
   <th style="text-align:right;"> p.value </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> (Intercept) </td>
   <td style="text-align:right;"> 24.500 </td>
   <td style="text-align:right;"> 1.96 </td>
   <td style="text-align:right;"> 12.503 </td>
   <td style="text-align:right;"> 0.000 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> large_ensemble1 </td>
   <td style="text-align:right;"> -7.833 </td>
   <td style="text-align:right;"> 2.53 </td>
   <td style="text-align:right;"> -3.097 </td>
   <td style="text-align:right;"> 0.009 </td>
  </tr>
</tbody>
</table>

**Repeat for all 37 musicians.**

---

Suppose a student **only** played in small ensembles/ solos.

- What is the estimated intercept for this student?  
  - What is the estimated slope for this student?

]

---

<div class="figure" style="text-align: center">
<img src="data:image/png;base64,#16-multilevel-models-pt2_files/figure-html/unnamed-chunk-9-1.png" alt="Recreated from BMLR Figure 8.9" width="70%" />
<p class="caption">Recreated from BMLR Figure 8.9</p>
</div>

**Now let's consider if there is an association between the estimated slopes, estimated intercepts, and the type of instrument**

---

## Level Two models

The slope and intercept for the `$i^{th}$` musician can be modeled as

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

Note the response variable in the Level Two models are not observed outcomes but the (fitted) slope and intercept from each musician

**See Part 3: Level Two Models in `lecture-16.Rmd`.**

.question[
How many observations used to estimate a? How many observations are used to estimate b?
[Click here](https://forms.gle/eusrCrTESmTJM1cP7) to submit your response.
]

---

## Estimated coefficients by instrument

---

## Level Two model

**Model for intercepts**

<table>
 <thead>
  <tr>
   <th style="text-align:left;"> term </th>
   <th style="text-align:right;"> estimate </th>
   <th style="text-align:right;"> std.error </th>
   <th style="text-align:right;"> statistic </th>
   <th style="text-align:right;"> p.value </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> (Intercept) </td>
   <td style="text-align:right;"> 16.283 </td>
   <td style="text-align:right;"> 0.671 </td>
   <td style="text-align:right;"> 24.249 </td>
   <td style="text-align:right;"> 0.000 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> orchestra1 </td>
   <td style="text-align:right;"> 1.411 </td>
   <td style="text-align:right;"> 0.991 </td>
   <td style="text-align:right;"> 1.424 </td>
   <td style="text-align:right;"> 0.163 </td>
  </tr>
</tbody>
</table>

**Model for slopes**

<table>
 <thead>
  <tr>
   <th style="text-align:left;"> term </th>
   <th style="text-align:right;"> estimate </th>
   <th style="text-align:right;"> std.error </th>
   <th style="text-align:right;"> statistic </th>
   <th style="text-align:right;"> p.value </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> (Intercept) </td>
   <td style="text-align:right;"> -0.771 </td>
   <td style="text-align:right;"> 0.851 </td>
   <td style="text-align:right;"> -0.906 </td>
   <td style="text-align:right;"> 0.373 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> orchestra1 </td>
   <td style="text-align:right;"> -1.406 </td>
   <td style="text-align:right;"> 1.203 </td>
   <td style="text-align:right;"> -1.168 </td>
   <td style="text-align:right;"> 0.253 </td>
  </tr>
</tbody>
</table>

---

## 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

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

.question[
- What is the predicted average performance anxiety before solos and small ensemble performances for vocalists and keyboardists? For those who place orchestral instruments?

- What is the predicted average performance anxiety before large ensemble performances for those who play orchestral instruments?
]

---

## 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 large ensemble performance)

⚠️  Does not share strength effectively across individuals

---

## Unified approach to two-level modeling

---

## Framework

**Level One**

`$$Y_{ij} = a_i + b_i ~ LargeEnsemble + \epsilon_{ij}$$`

**Level Two**

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

This approach uses likelihood-based methods (instead of least squares) to address the previously mentioned disadvantages

---

## Composite model

Plug in the equations for `$a_i$` and `$b_i$` to get the **composite model**

.eq[
`$$\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}$`

Note that we no longer need to estimate `$a_i$` and `$b_i$` directly as we did earlier. They conceptually connect the Level One and Level Two models.

---

## 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

---

<div class="figure" style="text-align: center">
<img src="data:image/png;base64,#16-multilevel-models-pt2_files/figure-html/unnamed-chunk-13-1.png" alt="Recreated from Figure 8.11" width="70%" />
<p class="caption">Recreated from Figure 8.11</p>
</div>

.question[
Describe what we learn about the association between the slopes and intercepts based on this plot. 
]

---

## Distribution of Level Two errors

Use a **multivariate normal** distribution for the Level Two error terms

.eq[
`$$\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$`?

---

## Visualizing multivariate normal distribution

<div class="figure" style="text-align: center">
<img src="data:image/png;base64,#16-multilevel-models-pt2_files/figure-html/contour-boundary-1.png" alt="Recreated from Figure 8.12" width="60%" />
<p class="caption">Recreated from Figure 8.12</p>
</div>

---

## Fit the model

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

```r
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

<table>
 <thead>
  <tr>
   <th style="text-align:left;"> effect </th>
   <th style="text-align:left;"> group </th>
   <th style="text-align:left;"> term </th>
   <th style="text-align:right;"> estimate </th>
   <th style="text-align:right;"> std.error </th>
   <th style="text-align:right;"> statistic </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> fixed </td>
   <td style="text-align:left;"> NA </td>
   <td style="text-align:left;"> (Intercept) </td>
   <td style="text-align:right;"> 15.930 </td>
   <td style="text-align:right;"> 0.641 </td>
   <td style="text-align:right;"> 24.833 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> fixed </td>
   <td style="text-align:left;"> NA </td>
   <td style="text-align:left;"> orchestra1 </td>
   <td style="text-align:right;"> 1.693 </td>
   <td style="text-align:right;"> 0.945 </td>
   <td style="text-align:right;"> 1.791 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> fixed </td>
   <td style="text-align:left;"> NA </td>
   <td style="text-align:left;"> large_ensemble1 </td>
   <td style="text-align:right;"> -0.911 </td>
   <td style="text-align:right;"> 0.845 </td>
   <td style="text-align:right;"> -1.077 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> fixed </td>
   <td style="text-align:left;"> NA </td>
   <td style="text-align:left;"> orchestra1:large_ensemble1 </td>
   <td style="text-align:right;"> -1.424 </td>
   <td style="text-align:right;"> 1.099 </td>
   <td style="text-align:right;"> -1.295 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> ran_pars </td>
   <td style="text-align:left;"> id </td>
   <td style="text-align:left;"> sd__(Intercept) </td>
   <td style="text-align:right;"> 2.378 </td>
   <td style="text-align:right;"> NA </td>
   <td style="text-align:right;"> NA </td>
  </tr>
  <tr>
   <td style="text-align:left;"> ran_pars </td>
   <td style="text-align:left;"> id </td>
   <td style="text-align:left;"> cor__(Intercept).large_ensemble1 </td>
   <td style="text-align:right;"> -0.635 </td>
   <td style="text-align:right;"> NA </td>
   <td style="text-align:right;"> NA </td>
  </tr>
  <tr>
   <td style="text-align:left;"> ran_pars </td>
   <td style="text-align:left;"> id </td>
   <td style="text-align:left;"> sd__large_ensemble1 </td>
   <td style="text-align:right;"> 0.672 </td>
   <td style="text-align:right;"> NA </td>
   <td style="text-align:right;"> NA </td>
  </tr>
  <tr>
   <td style="text-align:left;"> ran_pars </td>
   <td style="text-align:left;"> Residual </td>
   <td style="text-align:left;"> sd__Observation </td>
   <td style="text-align:right;"> 4.670 </td>
   <td style="text-align:right;"> NA </td>
   <td style="text-align:right;"> NA </td>
  </tr>
</tbody>
</table>

---

## Acknowledgements

The content in the slides is from 
  - [BMLR: Chapter 7 - Correlated data](https://bookdown.org/roback/bookdown-BeyondMLR/ch-corrdata.html)
  - [BMLR: Chapter 8 - Introduction to Multilevel Models](https://bookdown.org/roback/bookdown-BeyondMLR/ch-multilevelintro.html)

- 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.