Week 2: Longitudinal Data Analysis using Multilevel Modeling

class: center, middle, inverse, title-slide

# <b>Week 2: Longitudinal Data Analysis using Multilevel Modeling</b>
## Multivariate Statistics and Methodology using R (MSMR)<br><br>
### Dan Mirman
### Department of Psychology<br>The University of Edinburgh

---

# Longitudinal data are a natural application domain for MLM

* Longitudinal measurements are *nested* within subjects (by definition)
* Longitudinal measurements are related by a continuous variable, spacing can be uneven across participants, and data can be missing
    + These are problems rmANOVA
* Trajectories of longitudinal change can be nonlinear (we'll get to that next week)

## We've already seen some of this...

* Weight maintenance (Week 1 lab)
* Visual search (not longitudinal, but principles are the same)

---
# Another example

County-level percentage of adults who are physically active at recommended levels. Data from PHE for 2012-2015 (`PHE_MentalHealth.Rdata`)

```r
library(tidyverse)
library(lme4)
library(effects)
load("./data/PHE_MentalHealth.Rdata")
str(mh_phe)
```

```
## 'data.frame':	9422 obs. of  7 variables:
##  $ IndicatorID  : int  848 848 848 848 848 848 848 848 848 848 ...
##  $ IndicatorName: chr  "Depression: Recorded prevalence (aged 18+)" "Depression: Recorded prevalence (aged 18+)" "Depression: Recorded prevalence (aged 18+)" "Depression: Recorded prevalence (aged 18+)" ...
##  $ County       : chr  "Hartlepool" "Middlesbrough" "Redcar and Cleveland" "Stockton-on-Tees" ...
##  $ Region       : chr  "North East" "North East" "North East" "North East" ...
##  $ Timeperiod   : chr  "2013/14" "2013/14" "2013/14" "2013/14" ...
##  $ Value        : num  5.32 5.74 7.86 9.19 7.38 ...
##  $ Year         : num  2013 2013 2013 2013 2013 ...
```

---
# Indicators

```r
gt::gt(unique(mh_phe[, 1:2]))
```

<div id="wzrwmntbwk" style="overflow-x:auto;overflow-y:auto;width:auto;height:auto;">
<style>html {
  font-family: -apple-system, BlinkMacSystemFont, 'Segoe UI', Roboto, Oxygen, Ubuntu, Cantarell, 'Helvetica Neue', 'Fira Sans', 'Droid Sans', Arial, sans-serif;
}

#wzrwmntbwk .gt_table {
  display: table;
  border-collapse: collapse;
  margin-left: auto;
  margin-right: auto;
  color: #333333;
  font-size: 16px;
  font-weight: normal;
  font-style: normal;
  background-color: #FFFFFF;
  width: auto;
  border-top-style: solid;
  border-top-width: 2px;
  border-top-color: #A8A8A8;
  border-right-style: none;
  border-right-width: 2px;
  border-right-color: #D3D3D3;
  border-bottom-style: solid;
  border-bottom-width: 2px;
  border-bottom-color: #A8A8A8;
  border-left-style: none;
  border-left-width: 2px;
  border-left-color: #D3D3D3;
}

#wzrwmntbwk .gt_heading {
  background-color: #FFFFFF;
  text-align: center;
  border-bottom-color: #FFFFFF;
  border-left-style: none;
  border-left-width: 1px;
  border-left-color: #D3D3D3;
  border-right-style: none;
  border-right-width: 1px;
  border-right-color: #D3D3D3;
}

#wzrwmntbwk .gt_title {
  color: #333333;
  font-size: 125%;
  font-weight: initial;
  padding-top: 4px;
  padding-bottom: 4px;
  border-bottom-color: #FFFFFF;
  border-bottom-width: 0;
}

#wzrwmntbwk .gt_subtitle {
  color: #333333;
  font-size: 85%;
  font-weight: initial;
  padding-top: 0;
  padding-bottom: 4px;
  border-top-color: #FFFFFF;
  border-top-width: 0;
}

#wzrwmntbwk .gt_bottom_border {
  border-bottom-style: solid;
  border-bottom-width: 2px;
  border-bottom-color: #D3D3D3;
}

#wzrwmntbwk .gt_col_headings {
  border-top-style: solid;
  border-top-width: 2px;
  border-top-color: #D3D3D3;
  border-bottom-style: solid;
  border-bottom-width: 2px;
  border-bottom-color: #D3D3D3;
  border-left-style: none;
  border-left-width: 1px;
  border-left-color: #D3D3D3;
  border-right-style: none;
  border-right-width: 1px;
  border-right-color: #D3D3D3;
}

#wzrwmntbwk .gt_col_heading {
  color: #333333;
  background-color: #FFFFFF;
  font-size: 100%;
  font-weight: normal;
  text-transform: inherit;
  border-left-style: none;
  border-left-width: 1px;
  border-left-color: #D3D3D3;
  border-right-style: none;
  border-right-width: 1px;
  border-right-color: #D3D3D3;
  vertical-align: bottom;
  padding-top: 5px;
  padding-bottom: 6px;
  padding-left: 5px;
  padding-right: 5px;
  overflow-x: hidden;
}

#wzrwmntbwk .gt_column_spanner_outer {
  color: #333333;
  background-color: #FFFFFF;
  font-size: 100%;
  font-weight: normal;
  text-transform: inherit;
  padding-top: 0;
  padding-bottom: 0;
  padding-left: 4px;
  padding-right: 4px;
}

#wzrwmntbwk .gt_column_spanner_outer:first-child {
  padding-left: 0;
}

#wzrwmntbwk .gt_column_spanner_outer:last-child {
  padding-right: 0;
}

#wzrwmntbwk .gt_column_spanner {
  border-bottom-style: solid;
  border-bottom-width: 2px;
  border-bottom-color: #D3D3D3;
  vertical-align: bottom;
  padding-top: 5px;
  padding-bottom: 6px;
  overflow-x: hidden;
  display: inline-block;
  width: 100%;
}

#wzrwmntbwk .gt_group_heading {
  padding: 8px;
  color: #333333;
  background-color: #FFFFFF;
  font-size: 100%;
  font-weight: initial;
  text-transform: inherit;
  border-top-style: solid;
  border-top-width: 2px;
  border-top-color: #D3D3D3;
  border-bottom-style: solid;
  border-bottom-width: 2px;
  border-bottom-color: #D3D3D3;
  border-left-style: none;
  border-left-width: 1px;
  border-left-color: #D3D3D3;
  border-right-style: none;
  border-right-width: 1px;
  border-right-color: #D3D3D3;
  vertical-align: middle;
}

#wzrwmntbwk .gt_empty_group_heading {
  padding: 0.5px;
  color: #333333;
  background-color: #FFFFFF;
  font-size: 100%;
  font-weight: initial;
  border-top-style: solid;
  border-top-width: 2px;
  border-top-color: #D3D3D3;
  border-bottom-style: solid;
  border-bottom-width: 2px;
  border-bottom-color: #D3D3D3;
  vertical-align: middle;
}

#wzrwmntbwk .gt_from_md > :first-child {
  margin-top: 0;
}

#wzrwmntbwk .gt_from_md > :last-child {
  margin-bottom: 0;
}

#wzrwmntbwk .gt_row {
  padding-top: 8px;
  padding-bottom: 8px;
  padding-left: 5px;
  padding-right: 5px;
  margin: 10px;
  border-top-style: solid;
  border-top-width: 1px;
  border-top-color: #D3D3D3;
  border-left-style: none;
  border-left-width: 1px;
  border-left-color: #D3D3D3;
  border-right-style: none;
  border-right-width: 1px;
  border-right-color: #D3D3D3;
  vertical-align: middle;
  overflow-x: hidden;
}

#wzrwmntbwk .gt_stub {
  color: #333333;
  background-color: #FFFFFF;
  font-size: 100%;
  font-weight: initial;
  text-transform: inherit;
  border-right-style: solid;
  border-right-width: 2px;
  border-right-color: #D3D3D3;
  padding-left: 12px;
}

#wzrwmntbwk .gt_summary_row {
  color: #333333;
  background-color: #FFFFFF;
  text-transform: inherit;
  padding-top: 8px;
  padding-bottom: 8px;
  padding-left: 5px;
  padding-right: 5px;
}

#wzrwmntbwk .gt_first_summary_row {
  padding-top: 8px;
  padding-bottom: 8px;
  padding-left: 5px;
  padding-right: 5px;
  border-top-style: solid;
  border-top-width: 2px;
  border-top-color: #D3D3D3;
}

#wzrwmntbwk .gt_grand_summary_row {
  color: #333333;
  background-color: #FFFFFF;
  text-transform: inherit;
  padding-top: 8px;
  padding-bottom: 8px;
  padding-left: 5px;
  padding-right: 5px;
}

#wzrwmntbwk .gt_first_grand_summary_row {
  padding-top: 8px;
  padding-bottom: 8px;
  padding-left: 5px;
  padding-right: 5px;
  border-top-style: double;
  border-top-width: 6px;
  border-top-color: #D3D3D3;
}

#wzrwmntbwk .gt_striped {
  background-color: rgba(128, 128, 128, 0.05);
}

#wzrwmntbwk .gt_table_body {
  border-top-style: solid;
  border-top-width: 2px;
  border-top-color: #D3D3D3;
  border-bottom-style: solid;
  border-bottom-width: 2px;
  border-bottom-color: #D3D3D3;
}

#wzrwmntbwk .gt_footnotes {
  color: #333333;
  background-color: #FFFFFF;
  border-bottom-style: none;
  border-bottom-width: 2px;
  border-bottom-color: #D3D3D3;
  border-left-style: none;
  border-left-width: 2px;
  border-left-color: #D3D3D3;
  border-right-style: none;
  border-right-width: 2px;
  border-right-color: #D3D3D3;
}

#wzrwmntbwk .gt_footnote {
  margin: 0px;
  font-size: 90%;
  padding: 4px;
}

#wzrwmntbwk .gt_sourcenotes {
  color: #333333;
  background-color: #FFFFFF;
  border-bottom-style: none;
  border-bottom-width: 2px;
  border-bottom-color: #D3D3D3;
  border-left-style: none;
  border-left-width: 2px;
  border-left-color: #D3D3D3;
  border-right-style: none;
  border-right-width: 2px;
  border-right-color: #D3D3D3;
}

#wzrwmntbwk .gt_sourcenote {
  font-size: 90%;
  padding: 4px;
}

#wzrwmntbwk .gt_left {
  text-align: left;
}

#wzrwmntbwk .gt_center {
  text-align: center;
}

#wzrwmntbwk .gt_right {
  text-align: right;
  font-variant-numeric: tabular-nums;
}

#wzrwmntbwk .gt_font_normal {
  font-weight: normal;
}

#wzrwmntbwk .gt_font_bold {
  font-weight: bold;
}

#wzrwmntbwk .gt_font_italic {
  font-style: italic;
}

#wzrwmntbwk .gt_super {
  font-size: 65%;
}

#wzrwmntbwk .gt_footnote_marks {
  font-style: italic;
  font-weight: normal;
  font-size: 65%;
}
</style>
<table class="gt_table">
  
  <thead class="gt_col_headings">
    <tr>
      <th class="gt_col_heading gt_columns_bottom_border gt_right" rowspan="1" colspan="1">IndicatorID</th>
      <th class="gt_col_heading gt_columns_bottom_border gt_left" rowspan="1" colspan="1">IndicatorName</th>
    </tr>
  </thead>
  <tbody class="gt_table_body">
    <tr><td class="gt_row gt_right">848</td>
<td class="gt_row gt_left">Depression: Recorded prevalence (aged 18+)</td></tr>
    <tr><td class="gt_row gt_right">41001</td>
<td class="gt_row gt_left">Suicide rate</td></tr>
    <tr><td class="gt_row gt_right">90275</td>
<td class="gt_row gt_left">Percentage of physically active adults - historical method</td></tr>
    <tr><td class="gt_row gt_right">90646</td>
<td class="gt_row gt_left">Depression: QOF incidence (18+) - new diagnosis</td></tr>
  </tbody>
  
  
</table>
</div>

The indicator ID for percentage of physically active adults is 90275

---

```r
filter(mh_phe, IndicatorID == 90275) %>%
  ggplot(aes(Year, Value, color = Region)) +
    stat_summary(fun = mean, geom="point") +
    geom_smooth(method = "lm", se=FALSE)
```

```
## `geom_smooth()` using formula 'y ~ x'
```

![](msmr_lec02-2_LinearLDA_files/figure-html/unnamed-chunk-3-1.png)

Looks like there are some interesting differences. We can ask two kinds of questions:

1.  Did the regions differ in their *baseline (2012)* percentage physically active adults?
2.  Did the regions differ in their *slopes of change* of percentage physically active adults?

---
# What is the baseline?

**Did the regions differ in their *baseline (2012)* percentage physically active adults?**

We can answer this question using the **intercept** coefficients

But those will be estimated at `Year = 0`, and the question is not about physical activity the year Jesus was born.

So we need to adjust the time variable so that 2012 corresponds to time 0 (and we can select just the physical activity values while we're at it):

```r
active_dat <- filter(mh_phe, IndicatorID == 90275) %>%
  mutate(Time = Year - 2012)
```

Now `Time` is a variable just like `Year`, but going from 0 to 3 instead of 2012 to 2015

---
# Fit the models

```r
# base model: just change over time
m <- lmer(Value ~ Time + (Time | County), data = active_dat, REML = F)
# add baseline differences between regions
m.0 <- lmer(Value ~ Time + Region + (Time | County), data = active_dat, REML = F)
# add slope differences between regions
m.1 <- lmer(Value ~ Time * Region + (Time | County), data = active_dat, REML = F)
```

---
# Compare models

```r
anova(m, m.0, m.1)
```

```
## Data: active_dat
## Models:
## m: Value ~ Time + (Time | County)
## m.0: Value ~ Time + Region + (Time | County)
## m.1: Value ~ Time * Region + (Time | County)
##     npar  AIC  BIC logLik deviance Chisq Df Pr(>Chisq)    
## m      6 3299 3325  -1643     3287                        
## m.0   14 3275 3337  -1624     3247 39.34  8    4.3e-06 ***
## m.1   22 3281 3378  -1619     3237  9.86  8       0.28    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

Regions differed in baseline percentage of physically active adults ( `$\chi^2(8) = 39.3, p < 0.0001$` ), but there not in terms of the slope of change ( `$\chi^2(8) = 9.86, p = 0.28$` ).

---
# Visualize effects, Option 1: model-estimated percentages

```r
ef <- as.data.frame(effect("Time:Region", m.1))
ggplot(ef, aes(Time, fit, color=Region)) + 
  geom_line() +
  theme_bw() + scale_color_brewer(palette = "Set1")
```

![](msmr_lec02-2_LinearLDA_files/figure-html/unnamed-chunk-7-1.png)

---
# Visualize effects, Option 2: estimated parameters

[Side note: it took me a while to work out the data wrangling required to get the coefficients into a convenient form for plotting]

```r
ggplot(est, aes(Region, estimate, ymin=(estimate-std.error), ymax=(estimate+std.error))) +
  facet_wrap(~ Type, scales="free_x") +
  geom_pointrange() + coord_flip()
```

![](msmr_lec02-2_LinearLDA_files/figure-html/unnamed-chunk-9-1.png)