PSY 504: Advanced Statistics

class: center, middle, inverse, title-slide

.title[
# PSY 504: Advanced Statistics
]
.subtitle[
## Ordinal Regression: Why Dont You Likert Me?
]
.author[
### Jason Geller, Ph.D. (he/him/his)
]
.institute[
### Princeton University
]
.date[
### Updated:2023-02-19
]

---

## Ordinal Response Variables

- In psychology many variables have a natural ordering

- Grades (e.g., A,BC)
  - Education level (e.g., BA, MS, Phd)
  - Competitions (e.g., 1st, 2nd, 3rd)
  - Economic Status (e.g., wealthy, poor)
---
# This is a cat, not a dog?

.pull-left[

<img src="Honey.JPG" width="50%" style="display: block; margin: auto;" />
]

.pull-right[

1. Very likely to be a cat
2. Somewhat likely to be a cat
3. As likely to be cat or dog
4. Somewhat likely to be a dog
5. Very likely to be a dog
]

- This is a Likert scale ("Lick-ert")

---
# Cumulative link model

.pull-left[

- Theory: You have a continuous latent variable that has been categorized into bins:

- Latent "Honey's` catness"
]

.pull-right[
<img src="latent_ord.PNG" width="100%" style="display: block; margin: auto;" />

]
---
## Cumulative Probabilities - Notation

- `$P_k$`: Probability of being category k
- `$C_{pk}$`: Cumulative probability of being category k or lower
- `$1-C_{pk}$`: Probability of being above category k

- Notes:

1. `$C_{pk} = p_1 + p_2 + ... + p_k$`
  2. `$p_1 = C_{p1}, C_{pk} - C_{pk-1}(k=2,...K-1) p_k = 1-C_{pk-1}$`
  
---
## Cumulative Odds and Log-odds

- Odds

`$$\textrm{Odds} = \frac{\# \textrm{successes}}{\# \textrm{failures}}=
\frac{\# \textrm{successes}/n}{\# \textrm{failures}/n}=
\frac{p}{1-p}$$`
- Cumulative Odds

- *Odds* of being in at least in category k to above category k

`$$C_{pk}/1-C_{pk}$$`

- Log odds (cumulative logit)

`$$log(C_{pk}/1-C_{pk})$$`
---
## Ordinal regression model

`$$log (\frac{C_{pk}}{1-C_{pk}}) = \alpha - \beta_{j0}$$`
`$$\begin{array}{rcl} L_1 &=& \alpha_1-\beta_1x_1+\cdots+\beta_p X_p\\ L_2 &=& \alpha_2-\beta_1x_1+\cdots+\beta_p X_p & \\ L_{J-1} &=& \alpha_{J-1}-\beta_1x_1+\cdots+\beta_p X_p \end{array}$$`

- Here we are estimating J-1 equations simultaneously

- Each equation as a different intercept `$\alpha_k$` (thresholds) but a *common slope* `$\beta$`

- Intercepts are always ordered in size `$\alpha_1$` < `$\alpha_2$`

---
## Ordinal regression model

`$$\begin{array}{rcl} L_1 &=& \alpha_1-\beta_1x_1+\cdots+\beta_p X_p\\ L_2 &=& \alpha_2-\beta_1x_1+\cdots+\beta_p X_p& \\ L_{J-1} &=& \alpha_{J-1}-\beta_1x_1+\cdots+\beta_p X_p \end{array}$$`
- Where:

- `$\alpha$` (intercepts/thresholds/cut-offs) = Log-odds of falling into or below category

- `$\beta$` = Slope (constant between categories)

- `$-$` =  Helps with interpretation (positive `$b$` higher chance of being in higher categories)

---
## Proportional odds assumption

- Assumes slope is equal between categories

---
## Data: postgraduate school applications

- Undergraduate students report how likely they were to apply to graduate school: "Unlikely", "Somewhat Likely", "Very likely"

- Got additional information: GPA, parent education (college vs. no college), type of schooling (public vs. private)

```r
library(ordinal) # clm func for regression
library(MASS) # polr func for
```

```r
library(tidyverse)
library(emmeans)
library(ggeffects)
library(foreign) # read dta file

# load data
dat <- read.dta("https://stats.idre.ucla.edu/stat/data/ologit.dta")
```

---
## Data

|apply           | pared| public|  gpa|
|:---------------|-----:|------:|----:|
|very likely     |     0|      0| 3.26|
|somewhat likely |     1|      0| 3.21|
|unlikely        |     1|      1| 3.94|
|somewhat likely |     0|      0| 2.81|
|somewhat likely |     0|      0| 2.53|
|unlikely        |     0|      1| 2.59|

```r
dat <- dat %>%
  mutate(pared=as.factor(pared), public=as.factor(public))

# make sure ordered properly

head(ordered(dat$apply))
```

```
## [1] very likely     somewhat likely unlikely        somewhat likely
## [5] somewhat likely unlikely       
## Levels: unlikely < somewhat likely < very likely
```
---
## A simple model

`$$\text{logit}(p(y_i \leq j)) = \theta_j - \beta_2 \text{parent_education}_i$$`

```r
# link = probit would also be acceptable
ols1 = clm(apply ~ pared, data=dat, link = "logit")

summary(ols1)
```

```
## formula: apply ~ pared
## data:    dat
## 
##  link  threshold nobs logLik  AIC    niter max.grad cond.H 
##  logit flexible  400  -361.40 728.79 5(0)  1.25e-10 9.3e+00
## 
## Coefficients:
##        Estimate Std. Error z value Pr(>|z|)    
## pared1   1.1275     0.2634    4.28 1.87e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Threshold coefficients:
##                             Estimate Std. Error z value
## unlikely|somewhat likely      0.3768     0.1103   3.415
## somewhat likely|very likely   2.4519     0.1826  13.430
```
---
## Interpreting output

|term                             |  estimate| std.error| statistic|   p.value|coef.type |
|:--------------------------------|---------:|---------:|---------:|---------:|:---------|
|unlikely&#124;somewhat likely    | 0.3768424| 0.1103421|  3.415217| 0.0006373|intercept |
|somewhat likely&#124;very likely | 2.4518560| 0.1825629| 13.430201| 0.0000000|intercept |
|pared1                           | 1.1274910| 0.2634324|  4.280001| 0.0000187|location  |

- 1 coef and 2 thresholds - What's up with that?
  - Coefficients
  - Can be interpreted similarly to GLM
      - 1-unit increase (no edu -> edu) we expect a change of 1.13 on the log-odds scale
      - Means more likely to apply to college (go to right of scale)
---
## Interpreting output

.pull-left[
- Thresholds (cut-offs)

- Less than or equal to a certain level vs greater than that level

- j = 1: log-odds of rating = 1 vs. 2-3
 - j = 2: log-odds of rating = 1-2 vs. 3
]

.pull-right[
<img src="number_line.png" width="100%" style="display: block; margin: auto;" />
]
---
##  Cumulative odds ratios

- Sometimes odds ratios are more meaningful

- Almost 3x more likely  to apply to college if parent went to college

```r
knitr::kable(model_parameters(ols1, exponentiate=TRUE))
```

|Parameter                        | Coefficient|        SE|   CI|   CI_low|   CI_high|         z| df_error|         p|Component |
|:--------------------------------|-----------:|---------:|----:|--------:|---------:|---------:|--------:|---------:|:---------|
|unlikely&#124;somewhat likely    |    1.457674| 0.1608429| 0.95| 1.174187|  1.809606|  3.415217|      Inf| 0.0006373|intercept |
|somewhat likely&#124;very likely |   11.609874| 2.1195319| 0.95| 8.117639| 16.604481| 13.430201|      Inf| 0.0000000|intercept |
|pared1                           |    3.087899| 0.8134528| 0.95| 1.842591|  5.174843|  4.280001|      Inf| 0.0000187|location  |

---
## Probabilities

`$$p(logit)=\frac{e^{logit}}{1+e^{logit}}\frac{exp(a_k - bx)}{1+exp(a_k - bx)}$$`

```r
## view a summary of the model
ggpredictions_ols1 = ggemmeans(ols1,  terms = c("pared"))

ggpredictions_ols1 <- as.data.frame(ggpredictions_ols1)

knitr::kable(head(ggpredictions_ols1))
```

|x  | predicted| std.error|  conf.low| conf.high|response.level  |group |
|:--|---------:|---------:|---------:|---------:|:---------------|:-----|
|0  | 0.5931113| 0.0266289| 0.5409196| 0.6453030|unlikely        |1     |
|0  | 0.3275858| 0.0239325| 0.2806789| 0.3744926|somewhat likely |1     |
|0  | 0.0793029| 0.0133296| 0.0531773| 0.1054285|very likely     |1     |
|1  | 0.3206800| 0.0532744| 0.2162641| 0.4250959|unlikely        |1     |
|1  | 0.4692270| 0.0333495| 0.4038632| 0.5345908|somewhat likely |1     |
|1  | 0.2100930| 0.0424965| 0.1268014| 0.2933846|very likely     |1     |
---
## Model visualizations

.pull-left[

```r
#Note that ggpredicts doesn't give the original labels for positio - you need to give it the names of the factor labels, which will be in the order of the original model.
ggpredictions_ols1$x = factor(ggpredictions_ols1$x)
levels(ggpredictions_ols1$x) = c("No Edu", "Edu")
colnames(ggpredictions_ols1)[c(1, 6)] = c("Ed", "Apply")
#plot

plot <- ggplot(ggpredictions_ols1, aes(x = as.factor(Apply), y = predicted)) + geom_point(aes(color = Ed), position =position_dodge(width = 0.5)) + geom_errorbar(aes(ymin = conf.low, ymax = conf.high, color = Ed), position = position_dodge(width = 0.5), width = 0.3) + theme_bw() +   scale_y_continuous(labels = scales::percent)  + labs(x = "Apply", y = "predicted probability")
```
]

.pull-right[
<img src="Ordinal_files/figure-html/unnamed-chunk-15-1.png" width="100%" style="display: block; margin: auto;" />

]
---
## Plot ratings

- Predicted Probability

.pull-left[

```r
plot<- emmeans(ols1,
        ~ apply | pared, mode = "prob") %>% # model = prob get predicted probs
  plot() 
```
]
.pull-right[

```r
plot
```

]
---
## Testing proportional odds assumption

.pull-left[
- A few ways:

- Stratified binomial regressions

- Run separate binomial logistic regressions

```r
library(modelsummary) # plot multiple regressions or tables
dat$unlikely <- ifelse(dat$apply == "unlikely", 0, 1)
dat$likely <- ifelse(dat$apply == "very likely", 1, 0)
```
]

.pull-right[

<table class="table" style="width: auto !important; margin-left: auto; margin-right: auto;">
 <thead>
  <tr>
   <th style="text-align:left;">   </th>
   <th style="text-align:center;"> Ordinal 1 </th>
   <th style="text-align:center;">  Ordinal 2 </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> (Intercept) </td>
   <td style="text-align:center;"> −2.441 </td>
   <td style="text-align:center;"> −0.378 </td>
  </tr>
  <tr>
   <td style="text-align:left;">  </td>
   <td style="text-align:center;"> (0.201) </td>
   <td style="text-align:center;"> (0.111) </td>
  </tr>
  <tr>
   <td style="text-align:left;"> pared1 </td>
   <td style="text-align:center;"> 1.094 </td>
   <td style="text-align:center;"> 1.144 </td>
  </tr>
  <tr>
   <td style="text-align:left;box-shadow: 0px 1px">  </td>
   <td style="text-align:center;box-shadow: 0px 1px"> (0.370) </td>
   <td style="text-align:center;box-shadow: 0px 1px"> (0.292) </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Num.Obs. </td>
   <td style="text-align:center;"> 400 </td>
   <td style="text-align:center;"> 400 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> AIC </td>
   <td style="text-align:center;"> 256.2 </td>
   <td style="text-align:center;"> 538.1 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> BIC </td>
   <td style="text-align:center;"> 264.2 </td>
   <td style="text-align:center;"> 546.1 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Log.Lik. </td>
   <td style="text-align:center;"> −126.115 </td>
   <td style="text-align:center;"> −267.038 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> F </td>
   <td style="text-align:center;"> 8.719 </td>
   <td style="text-align:center;"> 15.292 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> RMSE </td>
   <td style="text-align:center;"> 0.30 </td>
   <td style="text-align:center;"> 0.49 </td>
  </tr>
</tbody>
</table>

]

---
## Test proporitional odds assumption

.pull-left[
- `brant` test

- Likelihood of the full ordinal logistic regression model (which makes the proportional odds assumption) to the likelihood of a reduced model that does not make this assumption
    
    - You want a ns `$\chi^2$` test
]

.pull-right[

```r
library(gofcat)# prop odds assum
#need to fit different model

brant.test(ols1)
```

```
## 
## Brant Test:
##            chi-sq   df   pr(>chi)
## Omnibus     0.017    1        0.9
## pared1      0.017    1        0.9
## 
## H0: Proportional odds assumption holds
```
]
---
## Test proporitional odds assumption

if test is violated, there are a few options:

- Baseline logistic model

- Use lowest level/rank as reference
  
- Adjacent category model

- Multinomial regression

---
## Test proporitional odds assumption

- Partial proportion odds model

```r
ols_nom = clm(apply ~ pared + public + gpa,nominal = ~ public, data=dat, link = "logit")

knitr::kable(tidy(ols_nom))
```

<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>
   <th style="text-align:left;"> coef.type </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> unlikely|somewhat likely.(Intercept) </td>
   <td style="text-align:right;"> 2.1658541 </td>
   <td style="text-align:right;"> 0.7798055 </td>
   <td style="text-align:right;"> 2.7774287 </td>
   <td style="text-align:right;"> 0.0054791 </td>
   <td style="text-align:left;"> intercept </td>
  </tr>
  <tr>
   <td style="text-align:left;"> somewhat likely|very likely.(Intercept) </td>
   <td style="text-align:right;"> 4.4106044 </td>
   <td style="text-align:right;"> 0.8088948 </td>
   <td style="text-align:right;"> 5.4526304 </td>
   <td style="text-align:right;"> 0.0000000 </td>
   <td style="text-align:left;"> intercept </td>
  </tr>
  <tr>
   <td style="text-align:left;"> unlikely|somewhat likely.public1 </td>
   <td style="text-align:right;"> 0.2350037 </td>
   <td style="text-align:right;"> 0.3052548 </td>
   <td style="text-align:right;"> 0.7698608 </td>
   <td style="text-align:right;"> 0.4413825 </td>
   <td style="text-align:left;"> intercept </td>
  </tr>
  <tr>
   <td style="text-align:left;"> somewhat likely|very likely.public1 </td>
   <td style="text-align:right;"> -0.5732669 </td>
   <td style="text-align:right;"> 0.4106292 </td>
   <td style="text-align:right;"> -1.3960695 </td>
   <td style="text-align:right;"> 0.1626936 </td>
   <td style="text-align:left;"> intercept </td>
  </tr>
  <tr>
   <td style="text-align:left;"> pared1 </td>
   <td style="text-align:right;"> 1.0576327 </td>
   <td style="text-align:right;"> 0.2665412 </td>
   <td style="text-align:right;"> 3.9679900 </td>
   <td style="text-align:right;"> 0.0000725 </td>
   <td style="text-align:left;"> location </td>
  </tr>
  <tr>
   <td style="text-align:left;"> public1 </td>
   <td style="text-align:right;"> NA </td>
   <td style="text-align:right;"> NA </td>
   <td style="text-align:right;"> NA </td>
   <td style="text-align:right;"> NA </td>
   <td style="text-align:left;"> location </td>
  </tr>
  <tr>
   <td style="text-align:left;"> gpa </td>
   <td style="text-align:right;"> 0.6105983 </td>
   <td style="text-align:right;"> 0.2607849 </td>
   <td style="text-align:right;"> 2.3413870 </td>
   <td style="text-align:right;"> 0.0192122 </td>
   <td style="text-align:left;"> location </td>
  </tr>
</tbody>
</table>

---
## Test ordinal assumptions

.pull-left[

- `sure` package: surrogate residuals

```r
library(sure)
library(cowplot)

# for reproducibility
set.seed(1225)

surrogate <- gridExtra::grid.arrange( 
  autoplot.clm(ols1, nsim = 100, what = "qq"),
  autoplot.clm(ols1, nsim = 100, what = "fitted", alpha = 0.5),
  autoplot.clm(ols1, nsim = 100, what = "covariate", x = dat$pared,
           xlab = "Education"),
  ncol = 2
)
```
]

.pull-right[

<img src="Ordinal_files/figure-html/unnamed-chunk-23-1.png" width="100%" style="display: block; margin: auto;" />
]
---
## Model 2: Add Public School + GPA

- Let's run this model:

```r
ols2 = clm(apply ~ pared + public + gpa, data=dat)
```

---
## Interpret the Coefficents

```r
knitr::kable(tidy(ols2))
```

<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>
   <th style="text-align:left;"> coef.type </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> unlikely|somewhat likely </td>
   <td style="text-align:right;"> 2.2033233 </td>
   <td style="text-align:right;"> 0.7795353 </td>
   <td style="text-align:right;"> 2.8264573 </td>
   <td style="text-align:right;"> 0.0047066 </td>
   <td style="text-align:left;"> intercept </td>
  </tr>
  <tr>
   <td style="text-align:left;"> somewhat likely|very likely </td>
   <td style="text-align:right;"> 4.2987674 </td>
   <td style="text-align:right;"> 0.8043147 </td>
   <td style="text-align:right;"> 5.3446338 </td>
   <td style="text-align:right;"> 0.0000001 </td>
   <td style="text-align:left;"> intercept </td>
  </tr>
  <tr>
   <td style="text-align:left;"> pared1 </td>
   <td style="text-align:right;"> 1.0476637 </td>
   <td style="text-align:right;"> 0.2657891 </td>
   <td style="text-align:right;"> 3.9417101 </td>
   <td style="text-align:right;"> 0.0000809 </td>
   <td style="text-align:left;"> location </td>
  </tr>
  <tr>
   <td style="text-align:left;"> public1 </td>
   <td style="text-align:right;"> -0.0586828 </td>
   <td style="text-align:right;"> 0.2978588 </td>
   <td style="text-align:right;"> -0.1970154 </td>
   <td style="text-align:right;"> 0.8438155 </td>
   <td style="text-align:left;"> location </td>
  </tr>
  <tr>
   <td style="text-align:left;"> gpa </td>
   <td style="text-align:right;"> 0.6157458 </td>
   <td style="text-align:right;"> 0.2606311 </td>
   <td style="text-align:right;"> 2.3625186 </td>
   <td style="text-align:right;"> 0.0181512 </td>
   <td style="text-align:left;"> location </td>
  </tr>
</tbody>
</table>

---
## Visualization: stacked area plots (continuous predictors)

.pull-left[

```r
library(effects) # stacked plots

stack <- plot(effect("gpa", ols2), style="stacked")
```
]

.pull-right[
<img src="Ordinal_files/figure-html/unnamed-chunk-27-1.png" width="100%" style="display: block; margin: auto;" />
]
---
## Visualization: stacked area plots (Categorical predictors)

.pull-left[

```r
library(effects) # stacked plots

stack <- plot(effect("public", ols2), style="stacked")
```
]

.pull-right[
<img src="Ordinal_files/figure-html/unnamed-chunk-29-1.png" width="100%" style="display: block; margin: auto;" />
]
---
## Model 3: Add public school + GPA interaction

```r
ols3 = clm(apply ~ public + pared*gpa, data=dat)
```

---
## Model 3: Add public school + GPA interaction

```r
knitr::kable(tidy(ols3))
```

<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>
   <th style="text-align:left;"> coef.type </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> unlikely|somewhat likely </td>
   <td style="text-align:right;"> 2.1310806 </td>
   <td style="text-align:right;"> 0.8463483 </td>
   <td style="text-align:right;"> 2.5179711 </td>
   <td style="text-align:right;"> 0.0118033 </td>
   <td style="text-align:left;"> intercept </td>
  </tr>
  <tr>
   <td style="text-align:left;"> somewhat likely|very likely </td>
   <td style="text-align:right;"> 4.2278488 </td>
   <td style="text-align:right;"> 0.8669972 </td>
   <td style="text-align:right;"> 4.8764271 </td>
   <td style="text-align:right;"> 0.0000011 </td>
   <td style="text-align:left;"> intercept </td>
  </tr>
  <tr>
   <td style="text-align:left;"> public1 </td>
   <td style="text-align:right;"> -0.0621122 </td>
   <td style="text-align:right;"> 0.2984302 </td>
   <td style="text-align:right;"> -0.2081296 </td>
   <td style="text-align:right;"> 0.8351278 </td>
   <td style="text-align:left;"> location </td>
  </tr>
  <tr>
   <td style="text-align:left;"> pared1 </td>
   <td style="text-align:right;"> 0.5875071 </td>
   <td style="text-align:right;"> 2.1304962 </td>
   <td style="text-align:right;"> 0.2757607 </td>
   <td style="text-align:right;"> 0.7827319 </td>
   <td style="text-align:left;"> location </td>
  </tr>
  <tr>
   <td style="text-align:left;"> gpa </td>
   <td style="text-align:right;"> 0.5918130 </td>
   <td style="text-align:right;"> 0.2826551 </td>
   <td style="text-align:right;"> 2.0937637 </td>
   <td style="text-align:right;"> 0.0362810 </td>
   <td style="text-align:left;"> location </td>
  </tr>
  <tr>
   <td style="text-align:left;"> pared1:gpa </td>
   <td style="text-align:right;"> 0.1484331 </td>
   <td style="text-align:right;"> 0.6819411 </td>
   <td style="text-align:right;"> 0.2176627 </td>
   <td style="text-align:right;"> 0.8276920 </td>
   <td style="text-align:left;"> location </td>
  </tr>
</tbody>
</table>

---
## Visualization: Interactions

.pull-left[

```r
interact <- ggemmeans(ols2, terms= c("gpa", "pared"))
```

]

.pull-right[

<img src="Ordinal_files/figure-html/unnamed-chunk-33-1.png" width="100%" style="display: block; margin: auto;" />
]
---
## Testing Significnce

- Likelihood ratio tests (LRT)

- Model comparisons
  
    - Main Effects vs. Interaction

```r
#main effects model vs. interaction
ols_test <- anova(ols2, ols3)

knitr::kable(ols_test)
```

<table>
 <thead>
  <tr>
   <th style="text-align:left;">   </th>
   <th style="text-align:right;"> no.par </th>
   <th style="text-align:right;"> AIC </th>
   <th style="text-align:right;"> logLik </th>
   <th style="text-align:right;"> LR.stat </th>
   <th style="text-align:right;"> df </th>
   <th style="text-align:right;"> Pr(&gt;Chisq) </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> ols2 </td>
   <td style="text-align:right;"> 5 </td>
   <td style="text-align:right;"> 727.0249 </td>
   <td style="text-align:right;"> -358.5124 </td>
   <td style="text-align:right;"> NA </td>
   <td style="text-align:right;"> NA </td>
   <td style="text-align:right;"> NA </td>
  </tr>
  <tr>
   <td style="text-align:left;"> ols3 </td>
   <td style="text-align:right;"> 6 </td>
   <td style="text-align:right;"> 728.9774 </td>
   <td style="text-align:right;"> -358.4887 </td>
   <td style="text-align:right;"> 0.0474444 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 0.8275715 </td>
  </tr>
</tbody>
</table>

---

```r
#main effects model vs. interaction
# USE TYPE III IF INTERACTIONS ARE IMPORTANT
ols_test <- car::Anova(ols3, type="III")

knitr::kable(ols_test)
```

<table>
 <thead>
  <tr>
   <th style="text-align:left;">   </th>
   <th style="text-align:right;"> Df </th>
   <th style="text-align:right;"> Chisq </th>
   <th style="text-align:right;"> Pr(&gt;Chisq) </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> public </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 6.3401786 </td>
   <td style="text-align:right;"> 0.0118033 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> pared </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 23.7795415 </td>
   <td style="text-align:right;"> 0.0000011 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> gpa </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 0.0433179 </td>
   <td style="text-align:right;"> 0.8351278 </td>
  </tr>
  <tr>
   <td style="text-align:left;"> pared:gpa </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 0.0760440 </td>
   <td style="text-align:right;"> 0.7827319 </td>
  </tr>
</tbody>
</table>
---
## Pairwise comparisons

- Main effects

```r
# pairwise contrasts
emmeans(ols3, list(pairwise ~ pared, pairwise ~ public))
```

```
## $`emmeans of pared`
##  pared emmean    SE  df asymp.LCL asymp.UCL
##  0     -1.436 0.169 Inf    -1.767    -1.104
##  1     -0.403 0.271 Inf    -0.934     0.128
## 
## Results are averaged over the levels of: public 
## Confidence level used: 0.95 
## 
## $`pairwise differences of pared`
##  1               estimate    SE  df z.ratio p.value
##  pared0 - pared1    -1.03 0.275 Inf  -3.759  0.0002
## 
## Results are averaged over the levels of: public 
## 
## $`emmeans of public`
##  public emmean    SE  df asymp.LCL asymp.UCL
##  0      -0.888 0.148 Inf     -1.18    -0.599
##  1      -0.950 0.295 Inf     -1.53    -0.373
## 
## Results are averaged over the levels of: pared 
## Confidence level used: 0.95 
## 
## $`pairwise differences of public`
##  1                 estimate    SE  df z.ratio p.value
##  public0 - public1   0.0621 0.298 Inf   0.208  0.8351
## 
## Results are averaged over the levels of: pared
```
---
## Simple effects

```r
# pairwise contrasts
emmeans(ols3,pairwise ~ pared*gpa)
```

```
## $emmeans
##  pared gpa emmean    SE  df asymp.LCL asymp.UCL
##  0       3 -1.436 0.169 Inf    -1.767    -1.104
##  1       3 -0.403 0.271 Inf    -0.934     0.128
## 
## Results are averaged over the levels of: public 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast                                                estimate    SE  df
##  pared0 gpa2.99892500132322 - pared1 gpa2.99892500132322    -1.03 0.275 Inf
##  z.ratio p.value
##   -3.759  0.0002
## 
## Results are averaged over the levels of: public
```

- What is wrong with this?
---
# Simple slopes

```r
emtrends(ols3,pairwise ~ pared, var="gpa")
```

```
## $emtrends
##  pared gpa.trend    SE  df asymp.LCL asymp.UCL
##  0         0.592 0.283 Inf    0.0378      1.15
##  1         0.740 0.629 Inf   -0.4922      1.97
## 
## Results are averaged over the levels of: public 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast        estimate    SE  df z.ratio p.value
##  pared0 - pared1   -0.148 0.682 Inf  -0.218  0.8277
## 
## Results are averaged over the levels of: public
```
---
## Model fit

- Pseudo-$R^2$

1 - `$LL_{mod} /LL_{null}$`

```r
library(easystats)

#model goodness
r2_mcfadden(ols2)
```

```
## # R2 for Generalized Linear Regression
##        R2: 0.033
##   adj. R2: 0.030
```
---
## Sample write-up

>  An ordered logit model was estimated to investigate factors (parent education, GPA, and public schooling) that influence whether undergraduates apply to  graduate school (“unlikely,” “somewhat likely,” “very likely”). Parent education predicted the likelihood of applyng to graduate school, *b* = 1.04, *z* = 3.942, *p* < .001, OR = 2.82.Students with parents that went to college were 4% more likely to apply to graduate school.  GPA was also a significant predictor, *b* = 0.615, *z* = 2.363, *p* < .001, OR = 1.84. Each point increase on GPA was associated with a 84% increase in the likelihood of applying to college. The overall  McFadden’s pseudo-R2 = .042. 
---
## Extensions

- Ordinal Regression in `brms`

- Bayesian implementation (tomororw)

- Familiar output
  
- Figure one liner

```r
ols2_brm = brm(ordered(apply) ~  gpa, data=dat, family = cumulative, cores = 4,chains = 4, backend = "cmdstanr")
```

---
## Multilevel Ordinal Regressions

- Repeated measures designs

- Clustered/nested designs

```r
ols2_clmm = clmm()
```