---
title: "Chapter 12: Categorical and continuous predictors"
author: "Andy Field"
date: "2025-04-29"
date-modified: today
format: 
  html:
    theme: journal
    embed-resources: true
    toc: true
execute: 
    message: false
    warning: false
---

```{r}
#| label: setup

library(easystats)
library(rstanarm)
library(tidyverse)
```

## Self-tests

### Load data

```{r}
#| label: st_data

# load data
puptreat_tib <- here::here("data/puppy_rct.csv") |>
  read_csv() |>
  mutate(
    dose = as_factor(dose) |> fct_relevel("No puppies", "15 mins", "30 mins")
  )


# alternatively
# puptreat_tib <- discovr::puppy_rct
```

### Means per group

```{r}
#| label: st1_solution
#| eval: false

puptreat_tib |> 
  group_by(dose) |> 
  describe_distribution() |> 
  display()
```

```{r}
#| label: st1_fancy_solution

puptreat_tib |> 
  group_by(dose) |> 
  describe_distribution() |> 
  select(dose:SD) |> 
  pivot_wider(
    names_from = "Variable",
    values_from = c("Mean", "SD")
  ) |> 
  display(column_names = c("Treatment", "*M* (Pre)", "*M* (Post)", "*SD* (Pre)", "*SD* (Post)"))
```

### Overall means


```{r}
#| label: grand_mean

puptreat_tib |> 
  describe_distribution(iqr = FALSE, range = FALSE) |>
  display()
```


### Violin plots

```{r}
#| label: violin_post

ggplot(puptreat_tib, aes(x = dose, y = post_happy, colour = dose)) +
  geom_point(position = position_jitter(width = 0.1), alpha = 0.6) +
  geom_violin(alpha = 0.2) + 
  stat_summary(fun.data = "mean_cl_normal", geom = "pointrange", position = position_dodge(width = 0.9)) +
  coord_cartesian(ylim = c(0, 10)) +
  scale_y_continuous(breaks = 0:10) +
  scale_colour_viridis_d(begin = 0.3, end = 0.8) +
  labs(x = "Puppy therapy group", y = "Post-therapy post_happy (0-10)", colour = "Puppy therapy group") +
  theme_minimal() +
  theme(legend.position = "none")
```


### Scatterplot of covariate

```{r}
#| label: scatter

ggplot(puptreat_tib, aes(x = pre_happy, y = post_happy)) +
  geom_smooth(method = "lm", colour = "#CC6677", fill = "#CC6677", alpha = 0.2) +
  geom_point(colour = "#882255") +
  coord_cartesian(ylim = c(0, 10), xlim = c(0, 10)) +
  scale_x_continuous(breaks = 0:10) +
  scale_y_continuous(breaks = 0:10) +
  labs(x = "Pre-treatment post_happy (0-10)", y = "Post-treatment post_happy (0-10)", colour = "Treatment", fill = "Treatment") +
 facet_wrap(~dose) + 
  theme_minimal()
```


### Sums of squares

```{r}
#| label: st3_model

# fit model with pre-treatment post_happy first
cov_first_lm <- lm(post_happy ~ pre_happy + dose, data = puptreat_tib)
anova(cov_first_lm) |> 
  model_parameters() |> 
  display()

# view the model parameters
model_parameters(cov_first_lm) |> 
  display(digits = 3)


# fit model with dose first
dose_first_lm <- lm(post_happy ~ dose + pre_happy, data = puptreat_tib)
anova(dose_first_lm) |> 
  model_parameters() |> 
  display()
```

### Self test

```{r}
#| label: self_test_1

pup_lm <- lm(post_happy ~ pre_happy + dose, data = puptreat_tib)
model_parameters(pup_lm) |> 
  display()
```


### Set contrasts

```{r}
#| label: contrasts

# set contrasts
puppy_vs_none <- c(-2/3, 1/3, 1/3)
long_vs_short <- c(0, -1/2, 1/2)

contrasts(puptreat_tib$dose) <- cbind(puppy_vs_none, long_vs_short)
```


### Fit the model taking a modelling approach: non-parallel slopes


```{r}
#| label: modelling_approach

intcpt_lm <- lm(post_happy ~ 1, data = puptreat_tib)
pre_lm <- lm(post_happy ~ pre_happy, data = puptreat_tib)
dose_lm <- lm(post_happy ~ pre_happy + dose, data = puptreat_tib)
interact_lm <- lm(post_happy ~ pre_happy + dose + dose:pre_happy, data = puptreat_tib)

test_wald(intcpt_lm, pre_lm, dose_lm, interact_lm) |> 
  display()

model_parameters(interact_lm) |> 
  display()

# The interaction doesn't signbificantly improve fit, so we might ultimately look at this model.

model_parameters(dose_lm) |> 
  display()
```

Alternatively use `update()`

```{r}
#| label: modelling_update
#| eval: false

# alternatively
intcpt_lm <- lm(post_happy ~ 1, data = puptreat_tib)
pre_lm <- update(intcpt_lm, .~. + pre_happy)
dose_lm <- update(pre_lm, .~. + dose)
interact_lm <- update(interact_lm, .~. + dose:pre_happy)

test_wald(intcpt_lm, pre_lm, dose_lm, interact_lm) |> 
  display()

model_parameters(interact_lm) |> 
  display()
```

### Testing homogeneity of regression slopes

```{r}
#| label: hors

main_effects_lm <- lm(post_happy ~ dose + pre_happy, data = puptreat_tib)
full_lm <- update(main_effects_lm, .~. + dose:pre_happy)

test_wald(main_effects_lm, full_lm) |> 
  display()
```

### Testing the independence of the treatment variable and covariate

```{r}
#| label: ind_treat

pre_lm <- lm(pre_happy ~ dose, data = puptreat_tib) 

test_wald(pre_lm) |> 
  display()
```

### Fit the model

```{r}
#| label: anova_style

puptreat_lm <- lm(post_happy ~ pre_happy + dose, data = puptreat_tib)
car::Anova(puptreat_lm, type = 3) |> 
  model_parameters(es_type = "omega", ci = 0.95) |> 
  display(use_symbols = TRUE)

# parameter estimates
model_parameters(puptreat_lm) |> 
  display()
```

### Evaluate the model


```{r}
#| label: modelling_check

check_model(puptreat_lm)
```

### Fit a robust model

```{r}
#| label: robust_lm

# parameter estimates CI and tests using heteroscedasticity consistent SEs
model_parameters(puptreat_lm, vcov = "HC4") |> 
  display()
```

### Adjusted means

```{r}
#| label: adj_m

estimate_means(puptreat_lm, by = c("dose")) |> 
  display()
```

### *Post hoc* tests

```{r}
#| label: post_hoc

estimate_contrasts(puptreat_lm, p_adjust = "holm") |> 
  display()
```




## Bayesian models

```{r}
#| label: rtb

# Bayesian model with default priors and bayes factors
puptreat_rtb <- stan_glm(post_happy ~ pre_happy + dose, data = puptreat_tib)

model_parameters(puptreat_rtb, test = c("pd", "bf"), null = 0) |> 
  display()
```


