Week 4, part 2: Crossed Random Effects

class: center, middle, inverse, title-slide

# <b>Week 4, part 2: Crossed Random Effects</b>
## Multivariate Statistics and Methodology using R (MSMR)<br><br>
### Dan Mirman
### Department of Psychology<br>The University of Edinburgh

---

# Crossed random effects

In most statistical tests, we are evaluating the reliability of some effect in the context of the variability in the data. The goal is to make an inference from the sample to the population from which from that sample was drawn.

Usually, we want to make inferences about reliability/variability across **subjects**: i.e., in general, do participants tend to show this effect or is it just one or two people showing it? That is, would we expect the rest of the participant population to also behave this way?

We might also ask about reliability/variability across **items**. This often comes up in laboratory experimental contexts where we ask participants to solve N problems, or answer N questions, or recognize N words, etc. There will be variability across those N "items" and we should test whether our effects of interest are reliable in the context of that variability. That is, would we expect the same outcome for other problems, questions, words, etc. from the same population?

Historically, this was done by conducting separate by-subjects ("F1") and by-items ("F2") analyses and journals sometimes even required this (a long time ago, I had a paper rejected because a key result was statistically significant by subjects but only marginal by items).

**Multilevel models provide a better solution to this problem.**

---
# Improving problem solving
Made-up data on an intervention to improve problem solving ability

```r
load("./data/problem_solving.Rdata")
summary(problem_solving)
```

```
##       Item      Prob_Type      Subject         Condition          RT       
##  AC     : 115   Hard:1494   S101   :  30   Control  :1397   Min.   :  440  
##  AV     : 109   Easy:1370   S109   :  30   Treatment:1467   1st Qu.: 1161  
##  AE     : 107               S121   :  30                    Median : 1536  
##  AB     : 105               S122   :  30                    Mean   : 1730  
##  AH     : 102               S124   :  30                    3rd Qu.: 2040  
##  AJ     : 102               S134   :  30                    Max.   :12753  
##  (Other):2224               (Other):2684
```

* `Item`: word problem, can be *`Hard`* or *`Easy`*
* `Prob_Type`: difficulty level of word problem (16 hard problems, 14 easy problems)
* `Subject`: Participant ID (N=120)
* `Condition`: whether the participant received the `Treatment` or not
* `RT`: time to solve the problem

Note: there is some missing data because trials where the participants failed to solve the problem are excluded.

---
# A traditional analysis approach

.pull-left[
Calculate subject-by-condition means

```r
ps_subj <- problem_solving %>%
  group_by(Subject, Condition, Prob_Type) %>%
  summarise(RT = mean(RT))
```
]

.pull-right[
![](msmr_lec04-2_RanEf_Crossed_files/figure-html/unnamed-chunk-3-1.png)
]

* Everyone solves `Easy` problems faster than `Hard` ones
* `Treatment` group seems faster at problem solving, esp. for `Hard` problems
---
# Repeated-measures ANOVA

```r
afex::aov_ez(id = "Subject", dv = "RT",
  within = "Prob_Type", between = "Condition",
  data = ps_subj)
```

```
## Anova Table (Type 3 tests)
## 
## Response: RT
##                Effect     df       MSE          F  ges p.value
## 1           Condition 1, 118 255247.21       2.44 .015    .121
## 2           Prob_Type 1, 118  98268.92 223.46 *** .345   <.001
## 3 Condition:Prob_Type 1, 118  98268.92       1.96 .005    .164
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
```

Looks like there is an overall problem difficulty effect, and no effect(s) of the intervention.

But hang on: not all word problems are the same and we're going to make inferences about solving problems **of this type**, not just about solving these particular problems.

---
# By-items analysis
Calculate item-by-condition means

```r
ps_item <- problem_solving %>%
  group_by(Item, Condition, Prob_Type) %>%
  summarise(RT = mean(RT))
# plot the item-level data
ggplot(ps_item, aes(Condition, RT, fill=Prob_Type)) + geom_boxplot()
```

![](msmr_lec04-2_RanEf_Crossed_files/figure-html/unnamed-chunk-5-1.png)

---
# By-items repeated-measures ANOVA

```r
afex::aov_ez(id = "Item", dv = "RT", 
             between = "Prob_Type", within = "Condition", 
             data = ps_item)
```

```
## Anova Table (Type 3 tests)
## 
## Response: RT
##                Effect    df       MSE         F  ges p.value
## 1           Prob_Type 1, 28 342623.14 17.47 *** .374   <.001
## 2           Condition 1, 28  15095.59  12.46 ** .018    .001
## 3 Prob_Type:Condition 1, 28  15095.59    5.26 * .008    .030
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
```

All effects significant. **Hol up, why the difference?**

---
# By-items repeated-measures ANOVA

* Problem type (difficulty) has a large effect, it is significant in both analyses
* Condition (intervention) effect and its interaction with problem type are small

Condition is between-subjects but within-items, so the between-subjects variability is strong in the by-subjects analysis but averaged away in the by-items analysis.

This makes the by-items analysis look overly strong (subject variability is missing) and the by-subjects analysis look overly weak (items consistency is missing).

Also, the idea that effects should be significant in separate by-items and by-subjects analyses (aka F1 and F2) is generally thought to be overly conservative.

---
# A multilevel modeling approach

Multilevel models provide a way to simultaneously model random variability at subject and item levels, as well as the group-level effects that we are interested in.

These data are nested (clustered observations that are not independent), but, unlike prior examples, they are not **hierarchical**. Items are clustered within subjects (subjects solved multiple problems); but also subjects are clustered within items (same problems were presented to multiple subjects). This can be modeled with **crossed random effects**.

```r
mod_ps <- lmer(RT ~ Prob_Type*Condition + 
                 (Prob_Type | Subject) + (Condition | Item), 
                 data=problem_solving, REML=FALSE)
```

```
## boundary (singular) fit: see ?isSingular
```

---
# A multilevel modeling approach

```r
gt(tidy(mod_ps))
```

<div id="kipbpfimvo" 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;
}

#kipbpfimvo .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;
}

#kipbpfimvo .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;
}

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

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

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

#kipbpfimvo .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;
}

#kipbpfimvo .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;
}

#kipbpfimvo .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;
}

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

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

#kipbpfimvo .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%;
}

#kipbpfimvo .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;
}

#kipbpfimvo .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;
}

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

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

#kipbpfimvo .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;
}

#kipbpfimvo .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;
}

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

#kipbpfimvo .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;
}

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

#kipbpfimvo .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;
}

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

#kipbpfimvo .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;
}

#kipbpfimvo .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;
}

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

#kipbpfimvo .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;
}

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

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

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

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

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

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

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

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

#kipbpfimvo .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_left" rowspan="1" colspan="1">effect</th>
      <th class="gt_col_heading gt_columns_bottom_border gt_left" rowspan="1" colspan="1">group</th>
      <th class="gt_col_heading gt_columns_bottom_border gt_left" rowspan="1" colspan="1">term</th>
      <th class="gt_col_heading gt_columns_bottom_border gt_right" rowspan="1" colspan="1">estimate</th>
      <th class="gt_col_heading gt_columns_bottom_border gt_right" rowspan="1" colspan="1">std.error</th>
      <th class="gt_col_heading gt_columns_bottom_border gt_right" rowspan="1" colspan="1">statistic</th>
      <th class="gt_col_heading gt_columns_bottom_border gt_right" rowspan="1" colspan="1">df</th>
      <th class="gt_col_heading gt_columns_bottom_border gt_right" rowspan="1" colspan="1">p.value</th>
    </tr>
  </thead>
  <tbody class="gt_table_body">
    <tr><td class="gt_row gt_left">fixed</td>
<td class="gt_row gt_left">NA</td>
<td class="gt_row gt_left">(Intercept)</td>
<td class="gt_row gt_right">1753.51</td>
<td class="gt_row gt_right">77.16</td>
<td class="gt_row gt_right">22.725</td>
<td class="gt_row gt_right">37.28</td>
<td class="gt_row gt_right">1.933e-23</td></tr>
    <tr><td class="gt_row gt_left">fixed</td>
<td class="gt_row gt_left">NA</td>
<td class="gt_row gt_left">Prob_Type1</td>
<td class="gt_row gt_right">325.24</td>
<td class="gt_row gt_right">73.50</td>
<td class="gt_row gt_right">4.425</td>
<td class="gt_row gt_right">30.88</td>
<td class="gt_row gt_right">1.115e-04</td></tr>
    <tr><td class="gt_row gt_left">fixed</td>
<td class="gt_row gt_left">NA</td>
<td class="gt_row gt_left">Condition1</td>
<td class="gt_row gt_right">58.34</td>
<td class="gt_row gt_right">29.43</td>
<td class="gt_row gt_right">1.982</td>
<td class="gt_row gt_right">111.64</td>
<td class="gt_row gt_right">4.989e-02</td></tr>
    <tr><td class="gt_row gt_left">fixed</td>
<td class="gt_row gt_left">NA</td>
<td class="gt_row gt_left">Prob_Type1:Condition1</td>
<td class="gt_row gt_right">36.14</td>
<td class="gt_row gt_right">17.73</td>
<td class="gt_row gt_right">2.038</td>
<td class="gt_row gt_right">86.73</td>
<td class="gt_row gt_right">4.462e-02</td></tr>
    <tr><td class="gt_row gt_left">ran_pars</td>
<td class="gt_row gt_left">Subject</td>
<td class="gt_row gt_left">sd__(Intercept)</td>
<td class="gt_row gt_right">272.06</td>
<td class="gt_row gt_right">NA</td>
<td class="gt_row gt_right">NA</td>
<td class="gt_row gt_right">NA</td>
<td class="gt_row gt_right">NA</td></tr>
    <tr><td class="gt_row gt_left">ran_pars</td>
<td class="gt_row gt_left">Subject</td>
<td class="gt_row gt_left">cor__(Intercept).Prob_Type1</td>
<td class="gt_row gt_right">1.00</td>
<td class="gt_row gt_right">NA</td>
<td class="gt_row gt_right">NA</td>
<td class="gt_row gt_right">NA</td>
<td class="gt_row gt_right">NA</td></tr>
    <tr><td class="gt_row gt_left">ran_pars</td>
<td class="gt_row gt_left">Subject</td>
<td class="gt_row gt_left">sd__Prob_Type1</td>
<td class="gt_row gt_right">94.41</td>
<td class="gt_row gt_right">NA</td>
<td class="gt_row gt_right">NA</td>
<td class="gt_row gt_right">NA</td>
<td class="gt_row gt_right">NA</td></tr>
    <tr><td class="gt_row gt_left">ran_pars</td>
<td class="gt_row gt_left">Item</td>
<td class="gt_row gt_left">sd__(Intercept)</td>
<td class="gt_row gt_right">391.54</td>
<td class="gt_row gt_right">NA</td>
<td class="gt_row gt_right">NA</td>
<td class="gt_row gt_right">NA</td>
<td class="gt_row gt_right">NA</td></tr>
    <tr><td class="gt_row gt_left">ran_pars</td>
<td class="gt_row gt_left">Item</td>
<td class="gt_row gt_left">cor__(Intercept).Condition1</td>
<td class="gt_row gt_right">1.00</td>
<td class="gt_row gt_right">NA</td>
<td class="gt_row gt_right">NA</td>
<td class="gt_row gt_right">NA</td>
<td class="gt_row gt_right">NA</td></tr>
    <tr><td class="gt_row gt_left">ran_pars</td>
<td class="gt_row gt_left">Item</td>
<td class="gt_row gt_left">sd__Condition1</td>
<td class="gt_row gt_right">37.31</td>
<td class="gt_row gt_right">NA</td>
<td class="gt_row gt_right">NA</td>
<td class="gt_row gt_right">NA</td>
<td class="gt_row gt_right">NA</td></tr>
    <tr><td class="gt_row gt_left">ran_pars</td>
<td class="gt_row gt_left">Residual</td>
<td class="gt_row gt_left">sd__Observation</td>
<td class="gt_row gt_right">736.17</td>
<td class="gt_row gt_right">NA</td>
<td class="gt_row gt_right">NA</td>
<td class="gt_row gt_right">NA</td>
<td class="gt_row gt_right">NA</td></tr>
  </tbody>
  
  
</table>
</div>

---
# Plotting results

Computing means and SE is a little tricky for trial-level data

Conveniently, the `effects` package will do that based on the fitted model

```r
library(effects)
efx <- effect("Prob_Type:Condition", mod_ps) %>% as.data.frame(.)
efx
```

```
##   Prob_Type Condition  fit    se lower upper
## 1      Hard   Control 2173 120.5  1937  2410
## 2      Easy   Control 1450 120.4  1214  1687
## 3      Hard Treatment 1984 104.2  1780  2188
## 4      Easy Treatment 1406 101.4  1207  1605
```

---
# Plotting results

```r
ggplot(efx, aes(Condition, fit, color=Prob_Type)) + 
  geom_point(size=3) + 
  geom_linerange(aes(ymin=fit-se, ymax=fit+se), size=2) +
  geom_linerange(aes(ymin=lower, ymax=upper), size=1) +
  theme_bw() + labs(y="Response Time", color="Problem\nType")
```

![](msmr_lec04-2_RanEf_Crossed_files/figure-html/unnamed-chunk-11-1.png)

---
# Not so different after all
This week we covered some more complex random effect structures

* 3-level nesting
* Crossed random effects

These are just kinds of MLM and the principles from other lectures also apply to these analyses

* p-value (df) estimation using Satterthwaite method; model comparisons would've been a good alternative
* use full random effect structure, where "full" is defined by study design; can simplify if model doesn't converge by removing correlations and random "slopes"
* be aware of how your categorical variables as coded; can conduct pairwise comparisons using a single model
* use logistic regression for binary outcomes