---
title: "Chapter 15: Repeated measures"
author: "Andy Field"
date: "2025-05-22"
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(glmmTMB)
library(tidyverse)


# load data
train_tib <- here::here("data/dog_training.csv") |> 
  read_csv()

sniff_tib <- here::here("data/sniffer_dogs.csv") |>
  read_csv() |>
  mutate(
    entity = as_factor(entity)
  )

scent_tib <- here::here("data/alien_scents.csv") |>
  read_csv() |>
   mutate(
      entity = as_factor(entity) |> fct_relevel("Human", "Shapeshifter", "Alien"),
      scent_mask = as_factor(scent_mask) |> fct_relevel("None", "Human", "Fox")
      )

# or load from the discovr package
# train_tib <- discovr::dog_training
# sniff_tib <- discovr::sniffer_dogs
#scent_tib <- discovr::alien_scents

```

## Growth models

### Descriptive statistics

```{r}
#| label: growth_describe

train_tib |> 
  group_by(block) |> 
  describe_distribution(select = "vocalizations") |> 
  data_remove("Variable") |> # optional to remove redundant column
  display()
```

### Plot data

```{r}
#| label: plot
ggplot(train_tib, aes(x = block, y = vocalizations)) +
  geom_point(size = 1, alpha = 0.6, position = position_jitter(width = 0.1, height = 0.1), colour = "#CC6677") +
  geom_smooth(method = "lm", formula = y ~ x, alpha = 0.3, colour  = "#88CCEE", fill  = "#88CCEE") +
   geom_smooth(method = "lm", formula = y ~ poly(x, 2), alpha = 0.3, colour  = "#DDCC77", fill = "#DDCC77") +
  coord_cartesian(ylim = c(0, 100)) + 
  scale_y_continuous(breaks = seq(0, 100, 10)) +
  scale_x_continuous(breaks = c(1, 2, 3, 5)) +
  labs(x = "Blocks of 100 trials", y = "Vocalizations (% of trials)") +
  theme_minimal()
```

Plot a random selection of 10 cases

```{r}
set.seed(616)
ids <- unique(train_tib$id) |>  sample(size = 10)

train_tib |> 
  filter(id %in% ids) |> 
  ggplot(data = _, aes(x = block, y = vocalizations, colour = id)) +
    geom_smooth(method = "lm", formula = y ~ poly(x, 2), se = F, linewidth = 0.5) +
    geom_point(size = 1, position = position_jitter(width = 0.1, height = 0.1)) +
    coord_cartesian(ylim = c(0, 100)) + 
    scale_y_continuous(breaks = seq(0, 100, 10)) +
    scale_x_continuous(breaks = c(1, 2, 3, 5), labels = c("100", "200", "300", "500")) +
  scale_colour_viridis_d() +
    labs(x = "Block of trials", y = "Vocalizations (% of trials)", colour = "ID") +
    theme_minimal()
```

### Modelling linear change

```{r}
#| label: build_raw_growth

incpt_mlm <-glmmTMB(vocalizations ~ 1 + (1|id), data = train_tib)
block_mlm <- glmmTMB(vocalizations ~ block + (1|id), data = train_tib)
blockrs_mlm <- glmmTMB(vocalizations ~ block + (block|id), data = train_tib)

test_lrt(incpt_mlm, block_mlm, blockrs_mlm) |> 
  display()

model_performance(blockrs_mlm) |> 
  display()

model_parameters(blockrs_mlm) |> 
  display()

```

### Modelling non-linear change

### Build raw model

```{r}
#| label: build_raw_nonlin

quad_mlm <- glmmTMB(vocalizations ~ block + I(block^2) + (block|id), data = train_tib)

test_lrt(incpt_mlm, block_mlm, blockrs_mlm, quad_mlm) |> 
  display()

model_performance(quad_mlm) |> 
  display()

model_parameters(quad_mlm) |> 
  display()
```

### see what `poly()` does

```{r}
#| label: explore_poly

# see what poly() does
poly_tib <- train_tib |> 
  mutate(
    quad = block^2,
    block_linear = poly(block, 2)[, 1],
    block_quadratic = poly(block, 2)[, 2]
  ) |> 
  select(block, quad, block_linear, block_quadratic)

head(poly_tib) |>  display(digits = 3)
cor(poly_tib) |>  round(2) |>  display()
```

### Fit models using `poly()`

```{r}
#| label: build_poly_growth

incpt_mlm <- glmmTMB(vocalizations ~ 1 + (1|id), data = train_tib)
poly_mlm <- glmmTMB(vocalizations ~ poly(block, 2) + (1|id), data = train_tib)
polyrs_mlm <- glmmTMB(vocalizations ~ poly(block, 2) + (poly(block, 1)|id), data = train_tib)
polyrs2_mlm <- glmmTMB(vocalizations ~ poly(block, 2) + (poly(block, 2)|id), data = train_tib)


test_lrt(incpt_mlm, poly_mlm, polyrs_mlm, polyrs2_mlm) |>
  display()

model_performance(quad_mlm) |> 
  display()

model_parameters(polyrs2_mlm) |> display()

```

## One -way experimental designs

### Descriptive statistics (self-test)

```{r}
#| label: summary_by_entity

sniff_sum <- sniff_tib |> 
  group_by(entity) |> 
  describe_distribution(select = "vocalizations") 

sniff_sum |> 
  data_remove("Variable") |> # optional to remove redundant column
  display()
```

### Plot the data (self-test)

```{r}
#| label: sniff_plot

ggplot(sniff_tib, aes(x = entity, y = vocalizations)) +
  geom_violin(colour = "#999933", fill = "#DDCC77", alpha = 0.3) +
  stat_summary(fun.data = "mean_cl_normal", colour = "#999933") +
  scale_y_continuous(breaks = 0:12) +
  labs(y = "Vocalizations (in 1 minute)", x = "Entity sniffed") +
  theme_minimal()
```

### Fit the overall model

```{r}
#| label: sniff_mod

sniff_afx <- afex::aov_4(vocalizations ~ entity + (entity|dog_name), data = sniff_tib)

model_parameters(sniff_afx, es_type = "omega") |> 
  display(use_symbols = TRUE)
```

For pedagogic purposes only, refit using the Huynh–Feldt correction

```{r}
#| label: hf_epsilon

hf_afx <- afex::aov_4(vocalizations ~ entity + (entity|dog_name),
                      data = sniff_tib,
                      anova_table = list(correction = "HF"))

model_parameters(hf_afx, es_type = "omega") |> 
  display(use_symbols = TRUE)
```

### Contrasts

::: callout-caution
## Danger zone!

Given the main effect was not significant, this analysis is for pedagogic purposes only.
:::



```{r}
#| label: sniff_con

#check the order of entities in the data:

levels(sniff_tib$entity)

sniff_cons <- cbind(
  aliens_vs_non = c(1/2, -1/2, -1/2, 1/2),
  alien_vs_shape = c(1/2, 0, 0, -1/2),
  human_vs_manquin = c(0, 1/2, -1/2, 0)
  )

# inspect the contrast weights to check they're c orrect
sniff_cons


estimate_contrasts(sniff_afx, contrast = "entity", comparison = sniff_cons) |> 
  display()
```

### *Post hoc* tests

::: callout-caution
## Danger zone!

Given the main effect was not significant, this analysis is for pedagogic purposes only.
:::

```{r}
#| label: sniff_ph

estimate_contrasts(sniff_afx, contrast = "entity", p_adjust = "bonferroni") |> 
  display()
```

### Robust tests

```{r}
#| label: sniff_rob

WRS2::rmanova(
  y = sniff_tib$vocalizations,
  groups = sniff_tib$entity,
  blocks = sniff_tib$dog_name
  )
WRS2::rmmcp(
  y = sniff_tib$vocalizations,
  groups = sniff_tib$entity,
  blocks = sniff_tib$dog_name
  )
```

## Factorial repeated measures designs
### Inspect the data

Let's have a quick look at the structure of the data.

```{r}
#| label: scent_data

head(scent_tib, 15) |> 
  display()
```

### Descriptive statistics

```{r}
#| label: summary_scent

scent_sum <- scent_tib |> 
  group_by(entity, scent_mask) |> 
  describe_distribution(select = "vocalizations") 

scent_sum |> 
  data_remove("Variable") |> # optional to remove redundant column
  display()
```


### Fit the overall model

```{r}
#| label: scent_mod

scent_afx <- afex::aov_4(vocalizations ~ entity*scent_mask + (entity*scent_mask|dog_id), data = scent_tib)

model_parameters(scent_afx, es_type = "omega") |> 
  display(use_symbols = TRUE)
```


### Plot the data

```{r}
#| label: scent_gg

afex::afex_plot(object = scent_afx,
                x = "scent_mask",
                trace = "entity",
                error = "within",
                legend_title = "Entity") +
  labs(x = "Scent used", y = "Number of vocalizations") +
  theme_minimal()

```

### The main effect of **entity**

```{r}
#| label: emm_entity

estimate_means(scent_afx, by = "entity") |> 
  display()

# post hoc tests

estimate_contrasts(scent_afx, contrast = "entity", p_adjust = "bonferroni") |> 
  display()
```

### The main effect of **scent_mask**

```{r}
#| label: emm_scent_mask

estimate_means(scent_afx, by = "scent_mask") |> 
  display()

# post hoc tests

estimate_contrasts(scent_afx, contrast = "scent_mask", p_adjust = "bonferroni") |> 
  display()

```

### The interaction effect

```{r}
#| label: emm_interaction

# marginal means
estimate_means(scent_afx, by = c("scent_mask", "entity")) |> 
  display(column_names = c("Scent", "Entity", "Mean", "SE", "95% CI", "t(440)"))

# simple effect of entity within scent
estimate_contrasts(model = scent_afx, contrast = "scent_mask", by = "entity", comparison = "joint", p_adjust = "bonferroni") |> 
  display()

# simple effect of scent within entity
estimate_contrasts(model = scent_afx, contrast = "entity", by = "scent_mask", comparison = "joint", p_adjust = "bonferroni") |> 
  display()

# post hoc tests of entity within each scent

estimate_contrasts(model = scent_afx, contrast = "entity", by = "scent_mask", p_adjust = "bonferroni") |> 
  display()

```



