---
title: "Chapter 16: Mixed designs"
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
exercise_tib <- here::here("data/exercise.csv") |> 
  read_csv() |> 
  mutate(
    id = as_factor(id),
    intervention = as_factor(intervention) |>  fct_relevel("Wait list"),
    time = as_factor(time) |> fct_relevel("Baseline", "1 month", "6 months", "12 months")
  )

date_tib <- here::here("data/speed_date.csv") |>
  read_csv() |>
   mutate(
    strategy = as_factor(strategy) |> fct_relevel("Normal", "Hard to get"),
    looks = as_factor(looks) |> fct_relevel("Low", "Average", "High"),
    charisma = as_factor(charisma) |> fct_relevel("Low", "Average", "High"),
  )


# or load from the discovr package
# exercise_tib <- discovr::exercise
# date_tib <- discovr::speed_date


```

## Growth models

### Descriptive statistics

```{r}
#| label: growth_describe

exercise_tib |> 
  group_by(time, intervention) |> 
  describe_distribution(select = "wemwbs") |> 
  data_remove("Variable") |> # optional to remove redundant column
  display()
```

### Plot data

```{r}
#| label: plot
ggplot(exercise_tib, aes(time_num, wemwbs, colour = intervention, fill = intervention)) +
  geom_point(size = 1, alpha = 0.6, position = position_jitter(width = 0.2, height = 0.1)) +
  geom_smooth(method = "lm", alpha = 0.3) +
  coord_cartesian(ylim = c(0, 75)) +
  scale_y_continuous(breaks = seq(0, 75, 5)) +
  scale_x_continuous(breaks = c(0, 1, 6, 12), labels = c("0", "1", "6", "12")) +
  scale_colour_viridis_d(begin = 0.3, end = 0.85) +
  scale_fill_viridis_d(begin = 0.3, end = 0.85) +
  labs(x = "Time from baseline (months)", y = "Emotional well-being (WEMWBS)", colour = "Intervention", fill = "Intervention") +
  theme_minimal()
```

Plot a random selection of 10 cases

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

exercise_tib |> 
  filter(id %in% ids) |> 
  ggplot(aes(time_num, wemwbs, colour = id, fill = id)) +
  geom_point(size = 1, alpha = 0.6, position = position_jitter(width = 0.2, height = 0.1)) +
  geom_smooth(method = "lm", alpha = 0.3, se = F) +
  coord_cartesian(ylim = c(0, 75)) +
  scale_y_continuous(breaks = seq(0, 75, 5)) +
  scale_x_continuous(breaks = c(0, 1, 6, 12), labels = c("0", "1", "6", "12")) +
  scale_colour_viridis_d() +
  scale_fill_viridis_d() +
  facet_wrap(~ intervention) +
  labs(x = "Time from baseline (months)", y = "Emotional well-being (WEMWBS)", colour = "Intervention", fill = "Intervention") +
  theme_minimal() +
  theme(legend.position = "none")
```

### Modelling linear change

```{r}
#| label: build_raw_growth

incpt_mlm <-glmmTMB(wemwbs ~ 1 + (1|id), data = exercise_tib)
time_mlm <- glmmTMB(wemwbs ~ time_num + (1|id), data = exercise_tib)
timers_mlm <- glmmTMB(wemwbs ~ time_num + (time_num|id), data = exercise_tib)
ex_mlm <- update(timers_mlm, .~. + intervention)
int_mlm <- update(ex_mlm, .~. + time_num:intervention)

test_lrt(incpt_mlm, time_mlm, timers_mlm, ex_mlm, int_mlm) |> 
  display()

model_performance(int_mlm) |> 
  display()

model_parameters(int_mlm) |> 
  display()

estimate_slopes(model = int_mlm, trend = "time_num", by = "intervention") |> 
  display()

```


## Experimental designs
### Inspect the data

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

```{r}
#| label: date_data

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

### Descriptive statistics

```{r}
#| label: summary_date

date_sum <- date_tib |> 
  group_by(charisma, looks, strategy) |> 
  describe_distribution(select = "date") 

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


### Fit the overall model

```{r}
#| label: scent_mod

date_afx <- afex::aov_4(date ~ strategy*looks*charisma + (looks*charisma|id), data = date_tib)

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


### Plot the data

```{r}
#| label: scent_gg

afex::afex_plot(object = date_afx,
                x = "looks",
                trace = "charisma",
                error = "within",
                panel = "strategy",
                legend_title = "Charisma",
                mapping = c("shape", "color")) +
  scale_colour_viridis_d(begin = 0.3, end = 0.85) +
  labs(x = "Attractiveness", y = "Rating of date (0-100)") +
  theme_minimal()

```

### The main effect of **strategy**

```{r}
#| label: emm_strategy

estimate_means(date_afx, by = "strategy") |> 
  display()

# plot

afex::afex_plot(date_afx,
                x = "strategy",
                mapping = c("shape", "color"),
                data_plot = F,
                point_arg = list(size = 3))  +
  scale_colour_viridis_d(begin = 0.3, end = 0.85) +
  coord_cartesian(ylim = c(40,95)) +
  scale_y_continuous(breaks = seq(40, 95, 5)) +
  labs(x = "Dating strategy", y = "Rating of date (0-100)", colour  = "Dating strategy") +
  theme_minimal() +
  theme(legend.position = "none")

```




### The main effect of **looks**

```{r}
#| label: emm_attract

estimate_means(date_afx, by = "looks") |> 
  display()

# contrasts

estimate_contrasts(date_afx,
                   comparison = "trt.vs.ctrl",
                   contrast = "looks",
                   ref = 2,
                   p_adjust = "bonferroni",
                   backend = "emmeans") |> 
  display()

# plot

afex::afex_plot(date_afx,
                x = "looks",
                error = "within",
                data_plot = F,
                point_arg = list(size = 3))  +
  scale_colour_viridis_d(begin = 0.3, end = 0.85) +
  coord_cartesian(ylim = c(40,95)) +
  scale_y_continuous(breaks = seq(40, 95, 5)) +
  labs(x = "Attractiveness", y = "Rating of date (0-100)", colour  = "Attractiveness") +
  theme_minimal() +
  theme(legend.position = "none")

```


### The main effect of **charisma**

```{r}
#| label: emm_charisma

estimate_means(date_afx, by = "charisma") |> 
  display()

# contrasts

estimate_contrasts(date_afx,
                   comparison = "trt.vs.ctrl",
                   contrast = "charisma",
                   ref = 2,
                   p_adjust = "bonferroni",
                   backend = "emmeans") |> 
  display()

# plot

afex::afex_plot(date_afx,
                x = "charisma",
                mapping = c("shape", "color"),
                error = "within",
                data_plot = F,
                point_arg = list(size = 3))  +
  scale_colour_viridis_d(begin = 0.3, end = 0.85) +
  coord_cartesian(ylim = c(40,95)) +
  scale_y_continuous(breaks = seq(40, 95, 5)) +
  labs(x = "Charisma", y = "Rating of date (0-100)", colour  = "Charisma") +
  theme_minimal() +
  theme(legend.position = "none")

```


### The **strategy** by **looks** interaction

```{r}
#| label: emm_sl

# marginal means
estimate_means(date_afx, by = c("strategy", "looks")) |> 
  display(column_names = c("Strategy", "Looks", "Mean", "SE", "95% CI", "t(440)"))

# contrast

estimate_contrasts(
  date_afx,
  contrast = c("strategy", "looks"),
  interaction = c(strategy = "trt.vs.ctrl", looks = "trt.vs.ctrl"),
  ref = 2,
  p_adjust = "bonferroni",
  backend = "emmeans"
) |> display()

afex::afex_plot(date_afx,
                x = "looks",
                trace = "strategy",
                mapping = c("shape", "color"),
                data_plot = F,
                point_arg = list(size = 3),
                legend_title = "Dating strategy")  +
  scale_colour_viridis_d(begin = 0.3, end = 0.85) +
  coord_cartesian(ylim = c(40,95)) +
  scale_y_continuous(breaks = seq(40, 95, 5)) +
  labs(x = "Attractiveness", y = "Rating of date (0-100)") +
  theme_minimal() +
  theme(legend.position = "top")
```

### The **strategy** by **charisma** interaction

```{r}
#| label: emm_sp

# marginal means
estimate_means(date_afx, by = c("strategy", "charisma")) |> 
  display(column_names = c("Strategy", "Charisma", "Mean", "SE", "95% CI", "t(440)"))

# contrast

estimate_contrasts(
  date_afx,
  contrast = c("strategy", "charisma"),
  interaction = c(strategy = "trt.vs.ctrl", charisma = "trt.vs.ctrl"),
  ref = 2,
  p_adjust = "bonferroni",
  backend = "emmeans"
) |> display()

afex::afex_plot(date_afx,
                x = "charisma",
                trace = "strategy",
                mapping = c("shape", "color"),
                data_plot = F,
                point_arg = list(size = 3),
                legend_title = "Dating strategy")  +
  scale_colour_viridis_d(begin = 0.3, end = 0.85) +
  coord_cartesian(ylim = c(40,95)) +
  scale_y_continuous(breaks = seq(40, 95, 5)) +
  labs(x = "Charisma", y = "Rating of date (0-100)") +
  theme_minimal() +
  theme(legend.position = "top")
```

### The **looks** by **charisma** interaction

```{r}
#| label: emm_lp

# marginal means
estimate_means(date_afx, by = c("looks", "charisma")) |> 
  display(column_names = c("Looks", "Charisma", "Mean", "SE", "95% CI", "t(440)"))

# contrast

estimate_contrasts(
  date_afx,
  contrast = c("looks", "charisma"),
  interaction = c(looks = "trt.vs.ctrl", charisma = "trt.vs.ctrl"),
  ref = 2,
  p_adjust = "bonferroni",
  backend = "emmeans"
) |> display()

afex::afex_plot(date_afx,
                x = "looks",
                trace = "charisma",
                mapping = c("shape", "color"),
                data_plot = F,
                point_arg = list(size = 3),
                legend_title = "Charisma")  +
  scale_colour_viridis_d(begin = 0.3, end = 0.85) +
  coord_cartesian(ylim = c(40,95)) +
  scale_y_continuous(breaks = seq(40, 95, 5)) +
  labs(x = "Attractiveness", y = "Rating of date (0-100)") +
  theme_minimal() +
  theme(legend.position = "top")
```


### The **looks** by **charisma** by **strategy** interaction

```{r}

# marginal means
estimate_means(date_afx, by = c("looks", "charisma", "strategy")) |> 
  display(column_names = c("Looks", "Charisma", "Strategy", "Mean", "SE", "95% CI", "t(440)"))

# contrasts
estimate_contrasts(
  date_afx,
  contrast = c("strategy", "looks", "charisma"),
  interaction = c(strategy = "trt.vs.ctrl", looks = "trt.vs.ctrl", charisma = "trt.vs.ctrl"),
  ref = 2,
  p_adjust = "bonferroni",
  backend = "emmeans"
) |> display()


# plot

afex::afex_plot(date_afx,
                x = "looks",
                trace = "charisma",
                panel = "strategy",
                mapping = c("shape", "color"),
                error = "within",
                point_arg = list(size = 3), 
                legend_title = "Charisma")  +
  scale_colour_viridis_d(begin = 0.3, end = 0.85) +
  coord_cartesian(ylim = c(40,95)) +
  scale_y_continuous(breaks = seq(40, 95, 5)) +
  labs(x = "Attractiveness", y = "Rating of date (0-100)") +
  theme_minimal() +
  theme(legend.position = "bottom")

# simple effects

estimate_contrasts(date_afx,
                   contrast = "looks",
                   by = c("charisma", "strategy"),
                   comparison = "joint",
                   p_adjust = "bonferroni") |> 
  display()

```

