PSY 503: Foundations of Statistical Methods in Psychological Science

class: center, middle, inverse, title-slide

.title[
# PSY 503: Foundations of Statistical Methods in Psychological Science
]
.subtitle[
## Introduction to Multilevel Modeling
]
.author[
### Jason Geller, Ph.D. (he/him/his)
]
.institute[
### Princeton University
]
.date[
### Updated:2022-12-05
]

---

# Today

- Weekly Q&A

- Finish up logistic modeling

- Multilevel/Mixed Models

- Introduction to MLM (connection to GLM)
  
  - How to do it
    - Organizing data for MLM analysis
    - Specifying MLMs
    - Fit and interpret multilevel models
    - Visualizing Data
    - Reporting
---
# R Chunks

- Include this code chunk at top of RMD file

- Gets rid of all the messages and warnings

```r
knitr::opts_chunk$set(
  echo = TRUE,
  message = FALSE, 
  warning = FALSE)
```
---
# Binomial vs. Bernoulli Distributions

- Bernoulli = A binomial distribution with 1 trial (*N*=1)
.pull-left[
<table class="huxtable" style="border-collapse: collapse; border: 0px; margin-bottom: 2em; margin-top: 2em; ; margin-left: auto; margin-right: auto;  " id="tab:unnamed-chunk-2">
<col><col><tr>
<th style="vertical-align: top; text-align: left; white-space: normal; border-style: solid solid solid solid; border-width: 0.4pt 0pt 0.4pt 0.4pt;    padding: 6pt 6pt 6pt 6pt; font-weight: bold;">title</th><th style="vertical-align: top; text-align: left; white-space: normal; border-style: solid solid solid solid; border-width: 0.4pt 0.4pt 0.4pt 0pt;    padding: 6pt 6pt 6pt 6pt; font-weight: bold;">binary</th></tr>
<tr>
<td style="vertical-align: top; text-align: left; white-space: normal; border-style: solid solid solid solid; border-width: 0.4pt 0pt 0pt 0.4pt;    padding: 6pt 6pt 6pt 6pt; background-color: rgb(242, 242, 242); font-weight: normal;">21 & Over</td><td style="vertical-align: top; text-align: left; white-space: normal; border-style: solid solid solid solid; border-width: 0.4pt 0.4pt 0pt 0pt;    padding: 6pt 6pt 6pt 6pt; background-color: rgb(242, 242, 242); font-weight: normal;">FAIL</td></tr>
<tr>
<td style="vertical-align: top; text-align: left; white-space: normal; border-style: solid solid solid solid; border-width: 0pt 0pt 0pt 0.4pt;    padding: 6pt 6pt 6pt 6pt; font-weight: normal;">Dredd 3D</td><td style="vertical-align: top; text-align: left; white-space: normal; border-style: solid solid solid solid; border-width: 0pt 0.4pt 0pt 0pt;    padding: 6pt 6pt 6pt 6pt; font-weight: normal;">PASS</td></tr>
<tr>
<td style="vertical-align: top; text-align: left; white-space: normal; border-style: solid solid solid solid; border-width: 0pt 0pt 0pt 0.4pt;    padding: 6pt 6pt 6pt 6pt; background-color: rgb(242, 242, 242); font-weight: normal;">12 Years a Slave</td><td style="vertical-align: top; text-align: left; white-space: normal; border-style: solid solid solid solid; border-width: 0pt 0.4pt 0pt 0pt;    padding: 6pt 6pt 6pt 6pt; background-color: rgb(242, 242, 242); font-weight: normal;">FAIL</td></tr>
<tr>
<td style="vertical-align: top; text-align: left; white-space: normal; border-style: solid solid solid solid; border-width: 0pt 0pt 0pt 0.4pt;    padding: 6pt 6pt 6pt 6pt; font-weight: normal;">2 Guns</td><td style="vertical-align: top; text-align: left; white-space: normal; border-style: solid solid solid solid; border-width: 0pt 0.4pt 0pt 0pt;    padding: 6pt 6pt 6pt 6pt; font-weight: normal;">FAIL</td></tr>
<tr>
<td style="vertical-align: top; text-align: left; white-space: normal; border-style: solid solid solid solid; border-width: 0pt 0pt 0pt 0.4pt;    padding: 6pt 6pt 6pt 6pt; background-color: rgb(242, 242, 242); font-weight: normal;">42</td><td style="vertical-align: top; text-align: left; white-space: normal; border-style: solid solid solid solid; border-width: 0pt 0.4pt 0pt 0pt;    padding: 6pt 6pt 6pt 6pt; background-color: rgb(242, 242, 242); font-weight: normal;">FAIL</td></tr>
<tr>
<td style="vertical-align: top; text-align: left; white-space: normal; border-style: solid solid solid solid; border-width: 0pt 0pt 0.4pt 0.4pt;    padding: 6pt 6pt 6pt 6pt; font-weight: normal;">47 Ronin</td><td style="vertical-align: top; text-align: left; white-space: normal; border-style: solid solid solid solid; border-width: 0pt 0.4pt 0.4pt 0pt;    padding: 6pt 6pt 6pt 6pt; font-weight: normal;">FAIL</td></tr>
</table>

]
.pull-right[
- Each row here is one trial

- There is only one outcome for each person (*N*=1): pass or fail

```r
glm(pass~bechdel, data=data, family=binomial(link=logit))
```
]
---
# Binomial vs. Bernoulli

- Binomial

- Multiple trials per row

```r
 glm(cbind(Y,N-Y) ~ X, 
       family = binomial,data=DataFrame)
```

- Doesn't really matter which one you use

- I always use Bernoulli
---
# Drop-in-Deviance Test

- In standard logistic regression output, deviance test is analogous to F test in regression

- Can also use it for comparing different models to one another

---
# Interactions

- Use the interaction coefficient to check if there is a significant interaction

<table class="huxtable" style="border-collapse: collapse; border: 0px; margin-bottom: 2em; margin-top: 2em; ; margin-left: auto; margin-right: auto;  " id="tab:unnamed-chunk-6">
<col><col><col><col><col><tr>
<th style="vertical-align: top; text-align: left; white-space: normal; border-style: solid solid solid solid; border-width: 0.4pt 0pt 0.4pt 0.4pt;    padding: 6pt 6pt 6pt 6pt; font-weight: bold;">term</th><th style="vertical-align: top; text-align: right; white-space: normal; border-style: solid solid solid solid; border-width: 0.4pt 0pt 0.4pt 0pt;    padding: 6pt 6pt 6pt 6pt; font-weight: bold;">estimate</th><th style="vertical-align: top; text-align: right; white-space: normal; border-style: solid solid solid solid; border-width: 0.4pt 0pt 0.4pt 0pt;    padding: 6pt 6pt 6pt 6pt; font-weight: bold;">std.error</th><th style="vertical-align: top; text-align: right; white-space: normal; border-style: solid solid solid solid; border-width: 0.4pt 0pt 0.4pt 0pt;    padding: 6pt 6pt 6pt 6pt; font-weight: bold;">statistic</th><th style="vertical-align: top; text-align: right; white-space: normal; border-style: solid solid solid solid; border-width: 0.4pt 0.4pt 0.4pt 0pt;    padding: 6pt 6pt 6pt 6pt; font-weight: bold;">p.value</th></tr>
<tr>
<td style="vertical-align: top; text-align: left; white-space: normal; border-style: solid solid solid solid; border-width: 0.4pt 0pt 0pt 0.4pt;    padding: 6pt 6pt 6pt 6pt; background-color: rgb(242, 242, 242); font-weight: normal;">(Intercept)</td><td style="vertical-align: top; text-align: right; white-space: normal; border-style: solid solid solid solid; border-width: 0.4pt 0pt 0pt 0pt;    padding: 6pt 6pt 6pt 6pt; background-color: rgb(242, 242, 242); font-weight: normal;">4.88  </td><td style="vertical-align: top; text-align: right; white-space: normal; border-style: solid solid solid solid; border-width: 0.4pt 0pt 0pt 0pt;    padding: 6pt 6pt 6pt 6pt; background-color: rgb(242, 242, 242); font-weight: normal;">0.205  </td><td style="vertical-align: top; text-align: right; white-space: normal; border-style: solid solid solid solid; border-width: 0.4pt 0pt 0pt 0pt;    padding: 6pt 6pt 6pt 6pt; background-color: rgb(242, 242, 242); font-weight: normal;">23.8 </td><td style="vertical-align: top; text-align: right; white-space: normal; border-style: solid solid solid solid; border-width: 0.4pt 0.4pt 0pt 0pt;    padding: 6pt 6pt 6pt 6pt; background-color: rgb(242, 242, 242); font-weight: normal;">4.02e-82</td></tr>
<tr>
<td style="vertical-align: top; text-align: left; white-space: normal; border-style: solid solid solid solid; border-width: 0pt 0pt 0pt 0.4pt;    padding: 6pt 6pt 6pt 6pt; font-weight: normal;">age</td><td style="vertical-align: top; text-align: right; white-space: normal; padding: 6pt 6pt 6pt 6pt; font-weight: normal;">-0.0175</td><td style="vertical-align: top; text-align: right; white-space: normal; padding: 6pt 6pt 6pt 6pt; font-weight: normal;">0.00447</td><td style="vertical-align: top; text-align: right; white-space: normal; padding: 6pt 6pt 6pt 6pt; font-weight: normal;">-3.92</td><td style="vertical-align: top; text-align: right; white-space: normal; border-style: solid solid solid solid; border-width: 0pt 0.4pt 0pt 0pt;    padding: 6pt 6pt 6pt 6pt; font-weight: normal;">0.000103</td></tr>
<tr>
<td style="vertical-align: top; text-align: left; white-space: normal; border-style: solid solid solid solid; border-width: 0pt 0pt 0pt 0.4pt;    padding: 6pt 6pt 6pt 6pt; background-color: rgb(242, 242, 242); font-weight: normal;">gendermale</td><td style="vertical-align: top; text-align: right; white-space: normal; padding: 6pt 6pt 6pt 6pt; background-color: rgb(242, 242, 242); font-weight: normal;">-0.446 </td><td style="vertical-align: top; text-align: right; white-space: normal; padding: 6pt 6pt 6pt 6pt; background-color: rgb(242, 242, 242); font-weight: normal;">0.265  </td><td style="vertical-align: top; text-align: right; white-space: normal; padding: 6pt 6pt 6pt 6pt; background-color: rgb(242, 242, 242); font-weight: normal;">-1.68</td><td style="vertical-align: top; text-align: right; white-space: normal; border-style: solid solid solid solid; border-width: 0pt 0.4pt 0pt 0pt;    padding: 6pt 6pt 6pt 6pt; background-color: rgb(242, 242, 242); font-weight: normal;">0.0935  </td></tr>
<tr>
<td style="vertical-align: top; text-align: left; white-space: normal; border-style: solid solid solid solid; border-width: 0pt 0pt 0.4pt 0.4pt;    padding: 6pt 6pt 6pt 6pt; font-weight: normal;">age:gendermale</td><td style="vertical-align: top; text-align: right; white-space: normal; border-style: solid solid solid solid; border-width: 0pt 0pt 0.4pt 0pt;    padding: 6pt 6pt 6pt 6pt; font-weight: normal;">0.0135</td><td style="vertical-align: top; text-align: right; white-space: normal; border-style: solid solid solid solid; border-width: 0pt 0pt 0.4pt 0pt;    padding: 6pt 6pt 6pt 6pt; font-weight: normal;">0.00553</td><td style="vertical-align: top; text-align: right; white-space: normal; border-style: solid solid solid solid; border-width: 0pt 0pt 0.4pt 0pt;    padding: 6pt 6pt 6pt 6pt; font-weight: normal;">2.45</td><td style="vertical-align: top; text-align: right; white-space: normal; border-style: solid solid solid solid; border-width: 0pt 0.4pt 0.4pt 0pt;    padding: 6pt 6pt 6pt 6pt; font-weight: normal;">0.0148  </td></tr>
</table>

---
# ANOVA Assumptions

- Normality

```r
performance::check_normality()
```

- Homogenity of variance

```r
performance::check_homogeneity()
```

- Independence

---
# Odd's Ratios vs. Probability

- In logistic regression odd's ratios  are used as a measure of effect size

- Can be hard to gauge sometimes
  
    - Couch it within probability

---
# Convert Odd's Ratios to Cohen's D

```r
odds_ratios=c(1.5, 2.0, 3.0)

log_odds <- log(odds_ratios)

cohen_d <- sqrt(3)*log_odds

cohen_d
```

```
## [1] 0.7022862 1.2005661 1.9028523
```
---
# Logistic Activity

```r
data=read_csv("https://raw.githubusercontent.com/jgeller112/psy503-psych_stats/master/static/slides/14-Logistic_Regression/data/gss.csv")
```
Q: Does attitudes about premarital sex predict church attendence?

DV: Attended a religious service in the last week

Independent Variables:

- Age

- Sex

- Socioeconomic Status (SES)

- Negative attitudes toward  premarital sex

---

# Report Logistic Regression

- State your hypothesis, statistical test, its link function, and
justify the use of a logistic regression

> We hypothesized that people who were against premarital
sex would be more likely to have attended a religious event
in the previous week. Religious event attendance was a
binary outcome as the person either attended an event or
they did not, which violated the normality assumption
required for traditional regression. Thus, a logistic
regression with a logit link function was used to model
religious event attendance using R 4.0.4 (R Core Team,
2020).

---
# Report Logistic Regression

- State the full model and your results, including the Odds Ratio

> Religious event attendance was modelled as a function of attitudes about premarital sex, controlling for the age, sex, and
socioeconomic status of the respondent. As shown in Figure 1, this analysis revealed that attitudes about premarital sex
significantly predicted religious event attendance, b = 0.72, SE
= 0.05, z(1303) = 13.24, p < 0.01, Odds Ratio = 2.05:1. This implies that every 1-unit increase in negative attitudes toward premarital sex predicted a 2-fold increase in the likelihood that the person had attended religious services in the previous week. Table 1 provides the model estimates for all predictors.

---
class: middle

# A Paradox

- Simpson's Paradox

<img src="mlm_files/figure-html/unnamed-chunk-12-1.png" width="80%" />
---

```r
p <- prep %>% ggplot(aes(study,mcat)) 
ggplot(prep, aes(x = study, y = mcat))+
  geom_point2(size=10) +
  theme_lucid() + 
  geom_smooth(method="lm")
```

<img src="mlm_files/figure-html/unnamed-chunk-13-1.png" width="80%" />
---

```r
ggplot(prep, aes(x = study, y = mcat, color=classroom))+
  geom_point2(size=10) +
  theme_lucid() + 
   geom_smooth(method = "lm", se = FALSE)
```

<img src="mlm_files/figure-html/unnamed-chunk-14-1.png" width="80%" />
---
<img src="images/important.png" width="100%" style="display: block; margin: auto;" />
---
# One Analysis, Many Names

<img src="images/many_names.JPG" width="40%" style="display: block; margin: auto;" />
---
# What are Multilevel Models?
<br>
<br>
<br>
- A broad class of analyses that deal with hierarchy in your data

---
# What is a “Hierarchy?”

*Clustering = Nesting = Grouping = Hierarchies*

- Key idea: More than one dimension sampling simultaneously

- Repeated-measures and longitudinal designs

- “Nested” designs

- Any complex mixed design
  
*NO NEED FOR REPEATED ANOVA, MIXED ANOVA, ANCOVA*

---
# Repeated Measures Designs

<img src="images/repeat_trials.png" width="100%" style="display: block; margin: auto;" />
---
# Repeated Designs
<img src="images/repeat1.png" width="100%" style="display: block; margin: auto;" />
---
# Repeated Designs
<img src="images/repeat2.png" width="100%" style="display: block; margin: auto;" />
---
# Repeated Designs
<img src="images/repeat3.png" width="100%" style="display: block; margin: auto;" />
---
# Repeated Designs
<img src="images/repeat4.png" width="100%" style="display: block; margin: auto;" />
---
# Repeated Designs
<img src="images/repeat5.png" width="100%" style="display: block; margin: auto;" />
---
# Why MLM is Awesome #1

- Classic Analysis:

- Aggregate to level of the group  (e.g., with means)

- Drawback of classic analysis:

- Loss of resolution!
  - Loss of power!

- MLM Approach:

- Can analyze any complex data structure
---
# Longitudinal Designs

<img src="images/long_fixed.png" width="100%" style="display: block; margin: auto;" />
---
# Why MLM is Awesome #2

- Classic Analysis:

- Repeated-measures ANOVA

- Drawback of this approach:

- Missing Data

- Must exclude entire people OR, you can interpolate missing  data

- MLM Approach:

- Can analyze all the observations you have!

---
# Nested Designs

---
# Two-level Hierarchy

- Nested designs

<img src="images/nested.png" width="100%" style="display: block; margin: auto;" />
---
# Three-level Hierarchy

- For now we will focus on data with two levels:

- Level one: Most basic level of observation
  - Level two: Groups formed from aggregated level-one observation
---
# Crossed vs. Nested Designs

- Crossed Designs

---
# Why MLM is Awesome #3

- Classic Analysis:

- Repeated-measures ANOVA  
  
- Drawback of ANOVA:

- Only use categorical predictors

- MLM Approach:

- Can use any combo of categorical and continuous predictors!
---
# Broad Intro

- Multilevel Models are linear models with extra error terms

- *Understand regression/GLM, and you can do any MLM*
---
# LM

`$$y = b_0 + b_1*x + e$$`

---
# Explaining Variance is the Grand Prize

- Almost any classical statistic compares:

`$$\frac{Variance {Explained} }{Unexplained}$$`
---
# One Source of Variance

`$$y = b_0 + b_1*x + e$$`
<img src="images/glm.png" width="50%" style="display: block; margin: auto;" />

---
# Assumptions of LM

- Assumptions:

- Normality
  - Homosecdacisity
  - Independence
  - No Multicollinearity
  - No missing data
---
# Violations of Assumptions

- Transform them, if you can

- Use robust SEs

- Accept a decrease in statistical power, if you can

- Find a test that doesn’t require that assumption

---
# Key Assumptions Relative to MLM

- **All observations must be independent**

-  Residuals are expected to be independent of each  other and normally-distributed

- **No missing data**

---
# But The World is Interdependent!

<img src="images/interdep.png" width="100%" style="display: block; margin: auto;" />
---
# Multilevel Models

- When to use them:

- Within-subjects data
  
  - Mix of between- and within-subjects data
  
- Why use them:

- Captures variance occurring between groups and within groups
  
- What they are:

- Linear model with extra residuals
---
# Jumping Right In

- Words you hear constantly in  MLM Land:

- *Fixed effects*
  - *Random effects*
  - *Random intercepts*  
  - *Random slopes*
  
- What do they mean?

---
# Fixed and Random Effects

- **Fixed effect**:

- Population-level (i.e., average) effects that should persist across experiments
  
  - Usually experimental manipulations
  
  - Can be continuous or categorical
---
# Fixed and Random Effects

- **Random effects**:

- Randomly sampled observations over which you plan to generalize

- Participants
      - Words
      - Pictures
      
  - Can help account for individual variation
  
---
# Fixed and Random Effects

<img src="images/fixedxvsran.png" width="100%" style="display: block; margin: auto;" />
---
<img src="images/age.png" width="100%" style="display: block; margin: auto;" />
  
---

---
<img src="images/trial.png" width="100%" style="display: block; margin: auto;" />

---
## Practice

.question[
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...)
2. What is the response variable?
3. What are the fixed effects? What are the random effects?
]

<div class="countdown" id="timer_6a700b85" data-update-every="1" tabindex="0" style="right:0;bottom:0;margin:1.25%;">
<div class="countdown-controls"><button class="countdown-bump-down">−</button><button class="countdown-bump-up">+</button></div>
<code class="countdown-time"><span class="countdown-digits minutes">03</span><span class="countdown-digits colon">:</span><span class="countdown-digits seconds">00</span></code>
</div>
---
# OLS Approach

- Fixed effects model

---
# Random Intercepts

.pull-left[

```r
library(lme4)

lmer(DV ~ IV + IV2 (1|participant))
```

]

.pull-right[
<img src="mlm_files/figure-html/unnamed-chunk-38-1.png" width="100%" style="display: block; margin: auto;" />
]
---
# Random Intercepts and Slopes

*You should only use a random slope for a fixed effect if the fixed effect varies within the grouping factor*

.pull-left[

```r
library(lme4)

lmer(DV ~ IV + IV2 (1+IV1|participant))
```
]

.pull-right[
<img src="mlm_files/figure-html/unnamed-chunk-40-1.png" width="100%" style="display: block; margin: auto;" />
]
---
# Random Coefficents

- Estimated per group

- Level 1 units used to estimate each group’s random slopes & intercept

---
# How Do Your Groups Differ?

- Different averages

- **Random intercept**
      - Each group gets its own intercept
    
- Different relationships between x and y

- **Random Slope**
      - Each group gets its own slope
      - Only if var differs within subject
    
- Per-group slopes and intercepts are represented with new residuals

---

<img src="images/ranslopes.png" width="100%" style="display: block; margin: auto;" />
---
# Implications

- Multiple sources of variance?

- Just add more residuals!

- Each source of variance gets its own residual term

- Residuals capture variance

- Residuals are added rendering them conditionally independent

- Otherwise, MLM is the same as GLM
---
# How Are They “The Same"

- The fixed effects (usually) hold your hypothesis tests

- Fixed effects output:  Looks like GLM output

- Can (essentially) be interpreted like GLM output  For most people, this is all that matters

---
# How Are LM and MLM Different?

- MLM has random effects output

- Variance explained by the random effect

- This may be interesting, in and  amongst itself

- Fixed effects of random terms are the average estimates across groups

- Fixed effects and random effects are conditionally independent

---
class: middle

Note: 
# Golden Rule: Should have > 5 levels per random effect
---
# Number of Random Effects

The size of the groups affects how many random effects can be estimated for each group

Max random coefficients = ((# level 1 observations)/  group) - 1

---
# All Together

- In regression you just estimate the outcome, `$y$`

- In MLM, you estimate parameters on the right side of the  equation, too:

- Intercept: `$b_0$`  
  
  - Slopes: `$b_1$`, `$b_2$`
---

---
- Level 1: 
`$$y_{ij} = \beta_{0j} + \beta_{1j}x_{ij} + e_{ij}$$`
- Level 2:

- Multilevel Equation

- Predicted Intercept:

`$$\beta_{0j} = \beta_{0} + u_{0j}$$`

- Predicted Slope:

`$$\beta_{1j} = \beta_{1} + u_{1j}$$`
- Mixed Model Equation

`$$y_{ij} = (\beta_{0} + u_{0j}) + (\beta_{1} + u_{1j})x_{ij} + e_{ij}$$`
---
# Error Terms

- Captures all unexplained variance: `$y_{ij} - y$`

- We want to explain more of it by considering groups, `$y{ij}$` - `$y_j$`

- Since each group j has its own intercept and/or slope, you are more accurate at predicting $y_{ij} for any individual in the group

- Moreover, you are now accounting for the shared variance among group members

- Each new residual makes additional sources of variance conditionally independent