---
title: "Chapter 18: Categorical outcomes"
author: "Andy Field"
date: "2025-06-18"
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(ggmosaic)
library(tidyverse)


# load data

dance_tib <- here::here("data/cat_dance.csv") |> 
  read_csv() |> 
  mutate(
    dance = as_factor(dance) |> fct_relevel("No"),
    training = as_factor(training) |> fct_relevel("Food")
  )

catreg_tib <- here::here("data/cat_reg.csv") |> 
  read_csv()

animal_tib <- here::here("data/animal_dance.csv") |> 
  read_csv() |> 
  mutate(
    dance = as_factor(dance) |> fct_relevel("No"),
    training = as_factor(training) |> fct_relevel("Food"),
    animal = as_factor(animal) |> fct_relevel("Cat")
  )

# or load from the discovr package
#dance_tib <- discovr::cat_dance
#catreg_tib <- discovr::cat_reg
#animal_tib <- discovr::animal_dance
```

Inspect the data

```{r}
#| label: see_data

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

## Load data as frequencies

```{r}
#| label: cat_freq

# create the table
cat_freq <- tribble(
  ~training, ~dance, ~n,
  "Food", "No", 10,
  "Food", "Yes", 28,
  "Affection", "No", 114,
  "Affection", "Yes", 48,
)
# display the table
cat_freq |> 
  display()
```

Create these frequencies from raw data

```{r}
#| label: raw_to_freq
# create the table
cat_freq <- dance_tib |> 
  group_by(training, dance) |> 
  tally() |> 
  ungroup()

# display the table
cat_freq |> 
  display()
```

## Creating a contingency table

From the raw data

```{r}
#| label: cat_tbly_raw

# create the table
cat_xtbl <- dance_tib |>
  data_tabulate(select = "training", by = "dance", remove_na = TRUE)
 
# display the table
cat_xtbl |>
  display()
```

From  frequency data

```{r}
#| label: cat_tbly_freq

# create the table
cat_xtbl <- cat_freq |>
  data_tabulate(select = "training", by = "dance", weights = "n", remove_na = TRUE)

# display the table
cat_xtbl |> 
  display()
```

## Tables of proportions

### Total proportions

```{r}
#| label: prop_tot

# create the table
prop_total <- dance_tib |>
  data_tabulate(select = "training", by = "dance", proportions = "full", remove_na = TRUE)

# display the table
prop_total |> 
  display()
```

### Proportions within `dance`

```{r}
#| label: prop_dance

# create the table
prop_dance <- dance_tib |>
  data_tabulate(select = "training", by = "dance", proportions = "column", remove_na = TRUE)

prop_dance |> 
  display()
```

### Proportions within `training`


```{r}
#| label: prop_train

# create the table
prop_training <- dance_tib |>
  data_tabulate(select = "training", by = "dance", proportions = "row", remove_na = TRUE)

# display the table
prop_training |> 
  display()
```

## Chi-square test

```{r}
cat_chi <- cat_xtbl |> 
  as.table(simplify = TRUE) |> 
  chisq.test()


model_parameters(cat_chi) |> 
  display()


# check residual
get_residuals(cat_chi)  |> 
  display()

# check expected values
get_predicted(cat_chi)  |> 
  display()
```


## Odds ratio

```{r}
#| label: or

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


## Fisher's exact test

```{r}
#| label: catfish

cat_fish <- cat_xtbl |> 
  as.table(simplify = TRUE) |> 
  fisher.test()

#to view the results
model_parameters(cat_fish) |> 
  display()
```

## Bayesian measure of association

```{r}
#| label: bayes_factor

cat_bf <- cat_xtbl |> 
  as.table(simplify = TRUE) |> 
  as.matrix() |> 
  BayesFactor::contingencyTableBF(sampleType = "indepMulti", fixedMargin = "rows")

model_parameters(cat_bf) |> 
  display()
```


## Loglinear analysis

```{r}
#| label: see_animal_data

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

### Load data as frequencies

```{r}
#| label: animal_freq

# create the table
animal_freq <- tribble(
  ~animal, ~training, ~dance, ~n,
  "Cat", "Food", "No", 10,
  "Cat", "Food", "Yes", 28,
  "Cat", "Affection", "No", 114,
  "Cat", "Affection", "Yes", 48,
  "Dog", "Food", "No", 14,
  "Dog", "Food", "Yes", 20,
  "Dog", "Affection", "No", 7,
  "Dog", "Affection", "Yes", 29
)
# display the table
animal_freq |> 
  display()
```

### Create frequencies from the raw data

```{r}
#| label: animal_tbly_raw

# create the frequency table from raw scores
animal_freq <- animal_tib |>
  group_by(animal, training, dance) |> 
  tally() |> 
  ungroup()

```

### Create a contingency tables


```{r}
#| label: animal_tbly_freq

# create the table
animal_xtbl <- animal_freq |>
  pivot_wider(
    values_from = "n",
    names_from = "dance"
  )

# display the table
animal_xtbl |> 
  display()
```

### Fit the model

```{r}
#| label: models_animal

intcpt_glm <- glm(n ~ 1, family = "poisson", data = animal_freq)
dance_glm <- update(intcpt_glm, .~. + dance)
train_glm <- update(dance_glm, .~. + training)
all_main_glm <- update(train_glm, .~. + animal)
train_dance_glm <- update(all_main_glm, .~. + training:dance)
animal_dance_glm <- update(train_dance_glm, .~. + animal:dance)
all_two_way_glm <- update(animal_dance_glm, .~. + animal:training)
saturated_glm <- update(all_two_way_glm, .~. + dance:training:animal)

# compare blocks

compare_performance(saturated_glm, all_two_way_glm, all_main_glm, intcpt_glm, 
                    metrics = "common") |> 
  display()

# test blocks

test_lrt(saturated_glm, all_two_way_glm, all_main_glm, intcpt_glm) |> 
  display()

# In a parallel universe where the three-way interaction was not significant, test individual interactions

test_lrt(all_two_way_glm, animal_dance_glm, train_dance_glm, all_main_glm) |> 
  display()


model_parameters(saturated_glm) |> 
  display()

# Expected frequencies

estimate_expectation(saturated_glm) |> 
  display()

```


### Mosaic plot

```{r}
ggplot(data = animal_tib) +
  geom_mosaic(aes(x = product(training, dance), fill = animal)) + 
  scale_fill_viridis_d(begin = 0.3, end = 0.85) +
  labs(fill = "Animal", x = "Animal:Dance", y = "Type of training") + 
  theme_minimal()
```



### Follow-up on cats

This is identical to your earlier analysis.

### Follow-up on dogs

```{r}
#| label: dog_chi

dog_xtbl <- animal_freq |>
  filter(animal == "Dog") |>
  data_tabulate(select = "training", by = "dance", weights = "n", remove_na = TRUE)

# display the table
dog_xtbl |> 
  display()

# proportions within training

# create the table
prop_dog_train <- animal_freq |>
  filter(animal == "Dog") |>
  data_tabulate(select = "training", by = "dance", weights = "n", proportions = "row", remove_na = TRUE)
  

# display the table
prop_dog_train |> 
  display()


# chi-square test

dog_chi <- dog_xtbl |> 
  as.table(simplify = TRUE) |> 
  chisq.test()

model_parameters(dog_chi) |> 
  display()


# check residual
get_residuals(dog_chi) |> 
  display()

# check expected values
display(dog_chi$expected)

# odds ratios

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

```





## Self-test

```{r}
#| label: chi_lm

# saturated model
caturated_lm <- lm(ln_observed  ~ training*dance, data = catreg_tib)
model_parameters(caturated_lm) |> 
  display(digits = 3)

# un-saturated model

cat_lm <- lm(ln_expected  ~ training + dance, data = catreg_tib)
model_parameters(cat_lm) |> 
  display()

```

