---
title: "Chapter 19: Categorical outcomes"
author: "Andy Field"
date: "2025-06-27"
date-modified: today
format: 
  html:
    theme: journal
    embed-resources: true
    toc: true
    code-overflow: wrap
execute: 
    message: false
    warning: false
---

## Setup

```{r}
#| label: setup

library(easystats)
library(tidyverse)


# load data
eel_tib <- here::here("data/eel.csv") |> 
  read_csv() |> 
  mutate(
    cured = as_factor(cured) |>  fct_relevel("Not cured"),
    intervention = as_factor(intervention) |>  fct_relevel("No treatment"),
  )

penalty_tib <- here::here("data/penalty_shootout.csv") |> 
  read_csv() |> 
  mutate(
    scored = as_factor(scored) |>  fct_relevel("Missed penalty")
  )

# or load from the discovr package
#eel_tib <- discovr::eel
#penalty_tib <- discovr::penalty_shootout

```

## An example to make you feel eel


### Peek at the data

```{r}
#| label: eel_peek

head(eel_tib, 10) |> 
  display()
```



### Contingency table

```{r}
#| label: eel_cont

# create contingency table and display
cure_xtbl <- eel_tib |>
  data_tabulate(select = "intervention", by = "cured", remove_na = TRUE)

display(cure_xtbl)

# include proportion of total
eel_tib |>
  data_tabulate(select = "intervention", by = "cured", remove_na = TRUE, proportions = "row") |> 
  display(digits = 2)

# calculate odds ratio

cure_xtbl |> 
  as.table(simplify = TRUE) |> 
  oddsratio() |> 
  display()
```


### Fit the models


```{r}
#| label: eel_mod

intcpt_glm <- glm(cured ~ 1, family = binomial(), data = eel_tib)
intervention_glm <- update(intcpt_glm, .~. + intervention)
duration_glm <- update(intervention_glm, .~. + duration)
full_glm <- update(duration_glm, .~. + intervention:duration)

# compare model fit

compare_performance(intcpt_glm, intervention_glm, duration_glm, full_glm, metrics = "common") |> 
  display()

# test models

test_lrt(intcpt_glm, intervention_glm, duration_glm, full_glm) |> 
  display()

# based on this we'd choose intervention_glm

model_performance(intervention_glm) |> 
  display()
```

### Model parameters

```{r}
#| label: eel_params

# log odds units
model_parameters(intervention_glm) |> 
  display()

# odds units

model_parameters(intervention_glm, exponentiate = TRUE) |> 
  display()

```

### Predictions

```{r}
#| label: eel_pred

# predictions for individual cases
case_predictions <- estimate_prediction(intervention_glm)

slice_sample(.data = case_predictions, n = 10) |> 
  display()

# group-level predictions
estimate_means(model = intervention_glm, by = "intervention") |> 
  display()
```


:::{.callout-note}
## Bonus content!

The following code is for pedagogic purposes only.
:::

To see what Tjur's $R^2$ represents we can use the predicted values. This model has only 1 predictor, the intervention group. So we know that for everyone in the No treatment group, the predicted probability was 0.43 and for everyone in the intervention group it was 0.72. Let's add these predicted values to the data.

```{r}
#| label: Tjur_r2a

eel_tib <- eel_tib |> 
  mutate(pred_prob = ifelse(intervention == "No treatment", 0.43, 0.72))
```

Now, lets' calculated the average predicted probability for those that were cured and those that were not. We can do this by using `group_by()` to tell `R` to apply subsequent calculations separately for each category of the variable **cured**. The 'subsequent calculation' is to compute the mean of the predicted probabilities, which are stored in the variable `pred_prob`. 

```{r}
#| label: Tjur_r2b
eel_tib |> 
  group_by(cured) |> 
  summarize(
    mean_pred = mean(pred_prob, na.rm = T)
  ) |> 
  display(digits = 3)
```

The resulting values tell us that the average predicted probability for those that were cured was

$$
\overline{\hat{\pi}}_{\text{cured}} = 0.613
$$

and the average predicted probability for those that were not cured was

$$
\overline{\hat{\pi}}_{\text{Not cured}} = 0.527
$$

The resulting $R^2_T$ is

$$
\begin{aligned}
R^2_T &= \overline{\hat{\pi}}_{\text{Cured}} - \overline{\hat{\pi}}_{\text{Not cured}} \\
&= 0.613-0.527 \\
&= 0.086
\end{aligned}
$$


### Check residuals

```{r}
#| label: eel_resids

check_model(intervention_glm)

```

## A sporty example

### Peek at the data

```{r}
#| label: pen_peek

# show (the same) random sample of 10 cases. To see different random samples change the seed number
set.seed(356)
slice_sample(.data = penalty_tib, n = 10) |> 
  display()
```

### Fit the models


```{r}
#| label: pen_mod

penalty_glm <- glm(scored ~ ability*position, family = binomial(), data = penalty_tib)

# compare model fit

model_performance(penalty_glm, metrics = "common") |> 
  display()

# test model

test_lrt(penalty_glm) |> 
  display()
```

### Model parameters

```{r}
#| label: pen_params

# log odds units
model_parameters(penalty_glm) |> 
  display()

# odds units

model_parameters(penalty_glm, exponentiate = TRUE) |> 
  display()

```


### Predictions

```{r}
#| label: sporty_pred

# predictions for individual cases
case_predictions <- estimate_prediction(penalty_glm) 

slice_sample(.data = case_predictions, n = 10) |> 
  display()

# group-level predictions
pred_p_ap <- estimate_means(penalty_glm, by = c("ability", "position"))

display(pred_p_ap)

# plot these slopes

pred_p_ap <- estimate_means(penalty_glm, by = c("ability", "position"))

plot(pred_p_ap) +
  scale_colour_viridis_d() +
  scale_fill_viridis_d() +
  labs(x = "Player ability (1-10)", y = "Probability of scoring", colour = "Position", fill = "Position") +
  theme_minimal()

# plot the other way
pred_p_pa <- estimate_means(penalty_glm, by = c("position", "ability"))

plot(pred_p_pa, numeric_as_discrete = 15) +
  scale_colour_viridis_d() +
  scale_fill_viridis_d() +
  labs(x = "Position in shootout", y = "Probability of scoring", colour = "Player ability (1-10)", fill = "Player ability (1-10)") +
  theme_minimal()

```

### Simple slopes

```{r}
#| label: pen_plot

# simple slopes
estimate_slopes(model = penalty_glm,
                term = "ability",
                by = "position",
                predict = "link") |> 
  display()

# simple slopes as ORs
estimate_slopes(model = penalty_glm,
                term = "ability",
                by = "position",
                predict = "link",
                transform = exp) |> 
  display()
```


### Check residuals

```{r}
#| label: pen_resids

check_model(penalty_glm)

```

### Linearity of the logit

```{r}
#| label: lin_logit

logit_lm <- glm(scored ~ position + ability + position:log(position) + ability:log(ability), family = binomial(), data = penalty_tib)

model_parameters(logit_lm) |>  display()
```
