Miltilevel models

class: center, middle, title-slide, inverse
layout: false

# Multilevel models
## Professor Andy Field

<div>
  <img style="vertical-align:middle; width:30px; height:30px" src="media/twitter_60.png">
  <span style="line-height:40px;">@profandyfield</span>
</div>
<div>
  <img style="vertical-align:middle; width:60px" src="media/youtube.png">
  <span style="line-height:40px;">www.youtube.com/user/ProfAndyField/</span>
</div>
<div>
  <img style="vertical-align:middle; width:30px; height:30px" src="media/ds_com_fav.png">
  <span style="line-height:40px;">www.discoveringstatistics.com</span>
</div>
<div>
  <img style="vertical-align:middle; width:30px; height:30px" src="media/milton_grey_fav.png">
  <span style="line-height:40px;">www.milton-the-cat.rocks</span>
</div>
<div>
  <img style="vertical-align:middle; width:30px; height:30px" src="media/discovr_fav.png">
  <span style="line-height:40px;">www.discovr.rocks</span>
</div>

???
Music: WITTR: Prayer of Transformation

h or ?: Toggle the help window
j: Jump to next slide
k: Jump to previous slide
b: Toggle blackout mode
m: Toggle mirrored mode.
p: Toggle PresenterMode
f: Toggle Fullscreen
t: Reset presentation timer
<number> + <Return>: Jump to slide <number>
c: Create a clone presentation on a new window

---

class: center

![:scale 83%](media/spine_content_map.png)
???

We've seen this map of the process of fitting models before

---
class: center

![:scale 83%](media/spine_lec_04.png)

???

Today we focus back on the model itself to look at the form of the model we're fitting. The faded stuff still applies though - we'll look at bias, robust models, and of course samples and estimation and so on. But they are the same as for other models, what is different is the form of the model we're fitting.

---
#  Learning outcomes

* Understand what hierarchical data are
  + Why we can’t use the ‘normal’ GLM
  
* Understand fixed and random coefficients

* Understand how to build models

* Be able to conduct and interpret models of hierarchical data

---
# Hierarchical data

* Data structures are often hierarchical

+ Children nested within classrooms
  
  + Observations nested within people

+ Employees nested within organisations

+ Patients nested within hospitals
  
  + Patients nested within teams nested within hospitals
  
  + Service users nested within clinicians nested within hospitals nested within NHS trusts!
  
  + Zombies nested within rehabilitation clinics 😉

---
# A two-level hierarchy

.center[
![:scale 100%](media/two_level_hierarchy.png)
]

---
# A three-level hierarchy

.center[
![:scale 80%](media/three_level_hierarchy.png)
]
---
#  Why hierarchies matter

* Data from the same context will be more similar than data from different contexts

+ Children in the same class will perform more similarly than children from different classes
	
	+ People treated in the same clinics should be more similar in response than those treated at different clinics

* Lack of independence

+ Violates the assumption of spherical errors (specifically, independence)
	
	+ Biases SEs, CIs and p -values

---

![:abs 35%, 10%, 5%](media/mr_doubtful.png)

![:abs 35%, 10%, 50%](media/mr_doubtful_class.png)
![:abs 5%, 25%, 45%](media/down_arrow.png)

![:abs 25%, 15%, 43%](media/mr_doubtful_txt.png)

![:abs 35%, 60%, 5%](media/ms_confident.png)

![:abs 35%, 60%, 50%](media/ms_confident_class.png)
![:abs 5%, 75%, 45%](media/down_arrow.png)

![:abs 25%, 65%, 42%](media/ms_confident_txt.png)

---
#  Benefits of multilevel models

* Modelling variability in effects across contexts
	
	+ Model the variability in intercepts
	
	+ Model the variability in slopes
	
* Model violations of the assumption of spherical errors

+ Model differences in the variability of errors
	
	+ Model relationships between errors
		
	  + (Linear model for repeated observations – two weeks time!)
	
* Missing data
	
	+ MLMs cope well with missing data (in general)
	
---
#  A rehabilitative example

* **p_id**: Zombie id

* **clinic_id**: Which of 10 clinics the zombie attended to have the intervention (1 to 10)
	
* **intervention**: Was the participant assigned to waiting list (0) or gene therapy (1)?

* **zombification**: Initial level of zombification.

+ Were they originally assigned to the low (0) or high (1) level genetic enhancement program at JIG:SAW?

* **months_as_zombie**: Time spent in a zombified state (Months).

+ How long ago did they undergo the original genetic enhancement program (which zombified them)

* **resemblance**: Resemblance of their face to their pre-zombie state (0% to 100%)

---
#  Fixed and random coefficients

* Intercepts and slopes can be fixed or random
	
	+ In OLS regression they are fixed
	
* Fixed coefficients
	
	+ Intercepts/slopes are assumed to be the same across different contexts (in this case clinics)
	
* Random coefficients
  
  + Intercepts/slopes are allowed to vary across different contexts (in this case clinics)

---
#  Intercepts and slopes

![](multilevel_models_files/figure-html/unnamed-chunk-2-1.png)

---
#  The process

.center[
![:scale 75%](media/mlm_process.png)
]
---
.center[
![:scale 55%](media/rehab_data_structure.png)
]

---

.center[
![](multilevel_models_files/figure-html/unnamed-chunk-3-1.png)
]
---
#  The models (using familiar symbols)

* Fixed intercepts and slopes

.center[.eq_med[.ong[
$$
`\begin{aligned}
\text{resemblance}_{i} & = \hat{b}_{0} + \hat{b}_{1}\text{Intervention}_{i}+ e_{i}
\end{aligned}`
$$
]]]

* Random intercepts but fixed slopes

.center[.eq_med[.ong[
$$
`\begin{aligned}
\text{resemblance}_{ij} & = \hat{b}_{0j} + \hat{b}_{1}\text{Intervention}_{ij}+ e_{ij}\\
\hat{b}_{0j} & = \hat{b}_0 + \hat{u}_{0j}
\end{aligned}`
$$
]]]

* Random intercepts and slopes

.center[.eq_med[.ong[
$$
`\begin{aligned}
\text{resemblance}_{ij} & = \hat{b}_{0j} + \hat{b}_{1j}\text{Intervention}_{ij}+ e_{ij}\\
\hat{b}_{0j} & = \hat{b}_0 + \hat{u}_{0j} \\
\hat{b}_{1j} & = \hat{b}_1 + \hat{u}_{1j}
\end{aligned}`
$$
]]]

---
#  The models (using common notation)

* Fixed intercepts and slopes

.center[.eq_med[.ong[
$$
`\begin{aligned}
\text{resemblance}_{i} & = \hat{\gamma}_{0} + \hat{\gamma}_{1}\text{Intervention}_{i}+ e_{i}
\end{aligned}`
$$
]]]

* Random intercepts but fixed slopes

.center[.eq_med[.ong[
$$
`\begin{aligned}
\text{resemblance}_{ij} & = \hat{\gamma}_{0j} + \hat{\gamma}_{1}\text{Intervention}_{ij}+ e_{ij}\\
\hat{\gamma}_{0j} & = \hat{\gamma}_0 + \hat{u}_{0j}
\end{aligned}`
$$
]]]

* Random intercepts and slopes

.center[.eq_med[.ong[
$$
`\begin{aligned}
\text{resemblance}_{ij} & = \hat{\gamma}_{0j} + \hat{\gamma}_{1j}\text{Intervention}_{ij}+ e_{ij}\\
\hat{\gamma}_{0j} & = \hat{\gamma}_0 + \hat{u}_{0j} \\
\hat{\gamma}_{1j} & = \hat{\gamma}_1 + \hat{u}_{1j}
\end{aligned}`
$$
]]]

---
#  Adding predictors

* Add whether original enhancement was low or high intensity (**zombification**)

.center[.eq_med[.ong[
$$
`\begin{aligned}
\text{resemblance}_{ij} & = \hat{\gamma}_{0j} + \hat{\gamma}_{1j}\text{Intervention}_{ij}+ \hat{\gamma}_{2}\text{zombification}_{ij} + e_{ij}\\
\hat{\gamma}_{0j} & = \hat{\gamma}_0 + \hat{u}_{0j} \\
\hat{\gamma}_{1j} & = \hat{\gamma}_1 + \hat{u}_{1j}
\end{aligned}`
$$
]]]
	
* Add time since genetic enhancement (**months_as_zombie**)

.center[.ong[
$$
`\begin{aligned}
\text{resemblance}_{ij} & = \hat{\gamma}_{0j} + \hat{\gamma}_{1j}\text{Intervention}_{ij} + \hat{\gamma}_{2}\text{zombification}_{ij} + \hat{\gamma}_{3}\text{months_as_zombie}_{ij} + e_{ij}\\
\hat{\gamma}_{0j} & = \hat{\gamma}_0 + \hat{u}_{0j} \\
\hat{\gamma}_{1j} & = \hat{\gamma}_1 + \hat{u}_{1j}
\end{aligned}`
$$
]]

---
#  Assessing the fit and comparing models

* Assessing fit
	
	+ AIC
	
	+ BIC
	
* Models should be built up gradually
	
	+ Start with fixed coefficients
	
	+ Change 1 aspect of the model and compare to the previous with the change in the -2 log likelihood

---

### Model 1:

.eq_med[.ong[
$$
`\begin{aligned}
\text{resemblance}_{i} = \hat{\gamma}_{0}  + e_{i}
\end{aligned}`
$$
]]

```r
rehab_base <- nlme::gls(resemblance ~ 1,
                        data = rehab_tib,
                        method = "ML"
                        )
```

### Model 2

.eq_med[.ong[
$$
`\begin{aligned}
\text{resemblance}_{ij} & = \hat{\gamma}_{0} + \hat{u}_{0j} + e_{ij}
\end{aligned}`
$$
]]

```r
rehab_ri <- nlme::lme(resemblance ~ 1,
*                     random = ~1|clinic_id,
                      data = rehab_tib,
                      method = "ML"
                      )
```

---
### Model 3
.eq_med[.ong[
$$
`\begin{aligned}
\text{resemblance}_{ij} & = \hat{\gamma}_{0} + \hat{\gamma}_{1}\text{Intervention}_{ij} + \hat{u}_{0j} + e_{ij}
\end{aligned}`
$$
]]

```r
*rehab_int <- nlme::lme(resemblance ~ intervention,
                      random = ~1|clinic_id, 
                      data = rehab_tib,
                      method = "ML"
                      )
```

### Model 4

.eq_med[.ong[
$$
`\begin{aligned}
\text{resemblance}_{ij} & = \hat{\gamma}_{0} + \hat{\gamma}_{1}\text{Intervention}_{ij} + \hat{u}_{0j} + \hat{u}_{1j} + e_{ij}
\end{aligned}`
$$
]]

```r
rehab_rs <- nlme::lme(resemblance ~ intervention, 
*                     random = ~intervention|clinic_id,
                      data = rehab_tib,
                      method = "ML"
                      )
```

---
### Model 5

.eq_med[.ong[
$$
`\begin{aligned}
\text{resemblance}_{ij} & = \hat{\gamma}_{0} + \hat{\gamma}_{1}\text{Intervention}_{ij} + \hat{\gamma}_{2}\text{zombification}_{ij} + e_{ij} + \hat{u}_{0j} + \hat{u}_{1j}
\end{aligned}`
$$
]]

```r
*rehab_zom <- nlme::lme(resemblance ~ intervention + zombification,
                      random = ~intervention|clinic_id, 
                      data = rehab_tib,
                      method = "ML"
                      )
```

### Model 6
.eq_med[.ong[
$$
`\begin{aligned}
\text{resemblance}_{ij} &= \hat{\gamma}_{0} + \hat{\gamma}_{1}\text{Intervention}_{ij} + \hat{\gamma}_{2}\text{zombification}_{ij} + \hat{\gamma}_{3}\text{months_as_zombie}_{ij} \\
&\quad + e_{ij} + \hat{u}_{0j} + \hat{u}_{1j}
\end{aligned}`
$$
]]

```r
*rehab_maz <- nlme::lme(resemblance ~ intervention + zombification + months_as_zombie,
                      random = ~intervention|clinic_id, 
                      data = rehab_tib,
                      method = "ML"
                      )
```

---
#  Comparing models using R

```r
anova(rehab_base, rehab_ri, rehab_int, rehab_rs, rehab_zom, rehab_maz)
```

</br>

<table>
 <thead>
  <tr>
   <th style="text-align:right;"> Model </th>
   <th style="text-align:right;"> df </th>
   <th style="text-align:right;"> AIC </th>
   <th style="text-align:right;"> BIC </th>
   <th style="text-align:right;"> logLik </th>
   <th style="text-align:left;"> Test </th>
   <th style="text-align:right;"> L.Ratio </th>
   <th style="text-align:right;"> p-value </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 2 </td>
   <td style="text-align:right;"> 1724.75 </td>
   <td style="text-align:right;"> 1731.25 </td>
   <td style="text-align:right;"> -860.38 </td>
   <td style="text-align:left;">  </td>
   <td style="text-align:right;"> NA </td>
   <td style="text-align:right;"> NA </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 2 </td>
   <td style="text-align:right;"> 3 </td>
   <td style="text-align:right;"> 1715.50 </td>
   <td style="text-align:right;"> 1725.24 </td>
   <td style="text-align:right;"> -854.75 </td>
   <td style="text-align:left;"> 1 vs 2 </td>
   <td style="text-align:right;"> 11.25 </td>
   <td style="text-align:right;"> 0.00 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 3 </td>
   <td style="text-align:right;"> 4 </td>
   <td style="text-align:right;"> 1710.38 </td>
   <td style="text-align:right;"> 1723.36 </td>
   <td style="text-align:right;"> -851.19 </td>
   <td style="text-align:left;"> 2 vs 3 </td>
   <td style="text-align:right;"> 7.12 </td>
   <td style="text-align:right;"> 0.01 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 4 </td>
   <td style="text-align:right;"> 6 </td>
   <td style="text-align:right;"> 1702.11 </td>
   <td style="text-align:right;"> 1721.59 </td>
   <td style="text-align:right;"> -845.05 </td>
   <td style="text-align:left;"> 3 vs 4 </td>
   <td style="text-align:right;"> 12.27 </td>
   <td style="text-align:right;"> 0.00 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 5 </td>
   <td style="text-align:right;"> 7 </td>
   <td style="text-align:right;"> 1703.08 </td>
   <td style="text-align:right;"> 1725.81 </td>
   <td style="text-align:right;"> -844.54 </td>
   <td style="text-align:left;"> 4 vs 5 </td>
   <td style="text-align:right;"> 1.03 </td>
   <td style="text-align:right;"> 0.31 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 6 </td>
   <td style="text-align:right;"> 8 </td>
   <td style="text-align:right;"> 1103.53 </td>
   <td style="text-align:right;"> 1129.50 </td>
   <td style="text-align:right;"> -543.76 </td>
   <td style="text-align:left;"> 5 vs 6 </td>
   <td style="text-align:right;"> 601.55 </td>
   <td style="text-align:right;"> 0.00 </td>
  </tr>
</tbody>
</table>

---
#  Model parameters (Fixed effects)

```r
broom.mixed::tidy(rehab_maz, effects = "fixed")
```

</br>

<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;"> df </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;"> 80.933 </td>
   <td style="text-align:right;"> 1.599 </td>
   <td style="text-align:right;"> 177 </td>
   <td style="text-align:right;"> 50.622 </td>
   <td style="text-align:right;"> 0.000 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> interventionGene therapy </td>
   <td style="text-align:right;"> 6.793 </td>
   <td style="text-align:right;"> 2.706 </td>
   <td style="text-align:right;"> 177 </td>
   <td style="text-align:right;"> 2.511 </td>
   <td style="text-align:right;"> 0.013 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> zombificationHigh intensity </td>
   <td style="text-align:right;"> -4.297 </td>
   <td style="text-align:right;"> 0.548 </td>
   <td style="text-align:right;"> 177 </td>
   <td style="text-align:right;"> -7.843 </td>
   <td style="text-align:right;"> 0.000 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> months_as_zombie </td>
   <td style="text-align:right;"> -1.977 </td>
   <td style="text-align:right;"> 0.028 </td>
   <td style="text-align:right;"> 177 </td>
   <td style="text-align:right;"> -69.631 </td>
   <td style="text-align:right;"> 0.000 </td>
  </tr>
</tbody>
</table>

---
#  Model parameters (random effects)

```r
broom.mixed::tidy(rehab_maz, effects = "ran_pars")
```

</br>

<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>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> ran_pars </td>
   <td style="text-align:left;"> clinic_id </td>
   <td style="text-align:left;"> sd_(Intercept) </td>
   <td style="text-align:right;"> 4.47 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> ran_pars </td>
   <td style="text-align:left;"> clinic_id </td>
   <td style="text-align:left;"> cor_interventionGene therapy.(Intercept) </td>
   <td style="text-align:right;"> 0.66 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> ran_pars </td>
   <td style="text-align:left;"> clinic_id </td>
   <td style="text-align:left;"> sd_interventionGene therapy </td>
   <td style="text-align:right;"> 8.28 </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;"> 3.60 </td>
  </tr>
</tbody>
</table>

---
#  Random Intercepts (*SD* = 4.47)

.center[
![](multilevel_models_files/figure-html/unnamed-chunk-17-1.png)
]

---
#  Random Intercepts (*SD* = 4.47)

.center[
![](multilevel_models_files/figure-html/unnamed-chunk-18-1.png)
]

---

#  Random Intercepts (*SD* = 4.47)

.center[
![](multilevel_models_files/figure-html/unnamed-chunk-19-1.png)
]

---
#  Random slopes (*SD* = 8.28)

.center[
![](multilevel_models_files/figure-html/unnamed-chunk-20-1.png)
]

---
#  Random slopes (*SD* = 8.28)

.center[
![](multilevel_models_files/figure-html/unnamed-chunk-21-1.png)
]

---

#  Random slopes (*SD* = 8.28)

.center[
![](multilevel_models_files/figure-html/unnamed-chunk-22-1.png)
]

---
.pull-left[

```r
plot(rehab_maz)
```

![](multilevel_models_files/figure-html/unnamed-chunk-23-1.png)
]

.pull-right[

```r
qqnorm(rehab_maz)
```

![](multilevel_models_files/figure-html/unnamed-chunk-24-1.png)
]

---
#  Robust model parameters

```r
rehab_maz_rob <- robustlmm::rlmer(
  resemblance ~ intervention + zombification + months_as_zombie + (intervention|clinic_id),
  data = rehab_tib)
broom.mixed::tidy(rehab_maz_rob, effects = "fixed")
```

<table>
 <thead>
  <tr>
   <th style="text-align:left;"> effect </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;"> (Intercept) </td>
   <td style="text-align:right;"> 81.69 </td>
   <td style="text-align:right;"> 1.60 </td>
   <td style="text-align:right;"> 51.13 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> fixed </td>
   <td style="text-align:left;"> interventionGene therapy </td>
   <td style="text-align:right;"> 7.01 </td>
   <td style="text-align:right;"> 3.25 </td>
   <td style="text-align:right;"> 2.15 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> fixed </td>
   <td style="text-align:left;"> zombificationHigh intensity </td>
   <td style="text-align:right;"> -4.58 </td>
   <td style="text-align:right;"> 0.49 </td>
   <td style="text-align:right;"> -9.28 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> fixed </td>
   <td style="text-align:left;"> months_as_zombie </td>
   <td style="text-align:right;"> -1.99 </td>
   <td style="text-align:right;"> 0.03 </td>
   <td style="text-align:right;"> -77.85 </td>
  </tr>
</tbody>
</table>

---
#  Robust random effects

```r
broom.mixed::tidy(rehab_maz_rob, effects = "ran_pars")
```

<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>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> ran_pars </td>
   <td style="text-align:left;"> clinic_id </td>
   <td style="text-align:left;"> sd__(Intercept) </td>
   <td style="text-align:right;"> 4.47 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> ran_pars </td>
   <td style="text-align:left;"> clinic_id </td>
   <td style="text-align:left;"> cor__(Intercept).interventionGene therapy </td>
   <td style="text-align:right;"> 0.55 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> ran_pars </td>
   <td style="text-align:left;"> clinic_id </td>
   <td style="text-align:left;"> sd__interventionGene therapy </td>
   <td style="text-align:right;"> 9.83 </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;"> 3.18 </td>
  </tr>
</tbody>
</table>
---
#  To sum up …

* Data can be hierarchical and this hierarchical structure can be important.

+ The OLS linear model simply ignores the hierarchy.

* Hierarchical models are just a fancy linear model in which you estimate the variability in the slopes and intercepts within contexts

+ i.e. slopes and intercepts can be random variables (allowed to vary) rather than fixed (assumed to be equal in different situations).

* Start with a model that ignores the hierarchy and then add in random intercepts and slopes to see if they improve the fit of the model.