---
title: "Chapter 11: Comparing several independent means"
author: "Andy Field"
date: "2025-04-22"
date-modified: today
format: 
  html:
    theme: journal
    embed-resources: true
    toc: true
    code-overflow: wrap
execute: 
    message: false
    warning: false
---

```{r}
#| label: setup

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

# load data
puppy_tib <- here::here("data/puppies.csv") |>
  read_csv() |>
  mutate(
    dose = as_factor(dose)
  )

# or load from the discovr package
# puppy_tib <- discovr::puppies

```


A good idea to check that the levels of **dose** are in the order No puppies, 15 minutes, 30 minutes by executing:

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

# check the order of factor levels
levels(puppy_tib$dose)

# If they're not in the correct order, re-order them:
puppy_tib <- puppy_tib |>
  mutate(
    dose = fct_relevel(dose, "No puppies", "15 mins", "30 mins")
  )
```


## Self-test

```{r}
#| label: dummy_mod

puppy_tib <- puppy_tib |>
  mutate(
    long = ifelse(dose == "30 mins", 1, 0),
    short = ifelse(dose == "15 mins", 1, 0)
  )

puppy_lm <- lm(happiness ~ long + short, data = puppy_tib)
test_wald(puppy_lm) |> display()
model_parameters(puppy_lm) |> display()
```

## Self-test

Create the dummy variables:

```{r}
#| label: create_contrasts

puppy_tib <- puppy_tib |>
  mutate(
    contrast1 = ifelse(dose == "No puppies", -2/3, 1/3),
    contrast2 = case_when(
      dose == "No puppies" ~ 0, 
      dose == "15 mins" ~ -0.5,
      dose == "30 mins" ~ 0.5)
  )
```

Fit and inspect the model:

```{r}
#| label: fit_contrast_mod

puppy_con_lm <- lm(happiness ~ contrast1 + contrast2, data = puppy_tib)
test_wald(puppy_con_lm) |> display()
model_parameters(puppy_con_lm) |> display()
```

## Summary statistics and plot



```{r}
#| label: descriptives

# a violin plot
ggplot(puppy_tib, aes(dose, happiness)) + 
  geom_violin(fill = "#CC6677", colour = "#CC6677", alpha = 0.2) +
  stat_summary(fun.data = "mean_cl_boot", colour = "#882255") +
  labs(x = "Dose of puppies", y = "Happiness (0-10)") +
  scale_y_continuous(breaks = 1:7) +
  theme_minimal()

puppy_tib |> 
  group_by(dose) |> 
  describe_distribution(select = "happiness", ci = 0.95, iterations = 1000) |> 
  display()
```



## Set contrasts

This code sets a treatment contrast to the variable **dose** that compares each category to the last.

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

contrasts(puppy_tib$dose) <- contr.treatment(3, base = 3)

# We can look at the contrast by executing:
contrasts(puppy_tib$dose)
```


Let's set the actual planned contrasts from the chapter:

```{r}
#| label: set_contrasts

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

contrasts(puppy_tib$dose) <- cbind(puppy_vs_none, long_vs_short)
contrasts(puppy_tib$dose) # check the contrasts are set correctly
```


## Fit the model
### Assumptions met

```{r}
#| label: puppy_mod

puppy_lm <- lm(happiness ~ dose, data = puppy_tib, na.action = na.exclude)

# overall test of fit and effect sizes
anova(puppy_lm) |> 
  model_parameters(es_type = "omega", ci = 0.95) |> 
  display()

# view model parameters
model_parameters(puppy_lm) |> 
  display()

# plot residuals
check_model(puppy_lm)

# check Cooks

check_outliers(puppy_lm,
               method = "cook",
               threshold = list('cook' = 0.5))

# Bonferroni adjusted post hoc tests
estimate_contrasts(puppy_lm, p_adjust = "bonferroni", ci = 0.95) |>  
  display()
```

### Trend analysis

```{r}
#| label: trend

# set the contrast
contrasts(puppy_tib$dose) <- contr.poly(3)
# fit the model
puppy_trend <- lm(happiness ~ dose, data = puppy_tib, na.action = na.exclude)
# inspect parameters
model_parameters(puppy_trend) |> 
  display()
```

### Heteroscedasticity

Because of doing the trend analysis above, where we reset the contrast for dose, let's set it back to how it was for the main example.

```{r}
#| label: welch

# reset the contrasts
contrasts(puppy_tib$dose) <- cbind(puppy_vs_none, long_vs_short)

# Welch's F
oneway.test(happiness ~ dose, data = puppy_tib) |>
  model_parameters() |> 
  display()
# parameter estimates CI and tests using heteroscedasticity consistent SEs
model_parameters(puppy_lm, vcov = "HC4") |> 
  display()

```

### Robust parameter estimates:

```{r}
#| label: rob_pe

# estimation based on M-estimators
puppy_rob <- robust::lmRob(happiness ~ dose, data = puppy_tib)
summary(puppy_rob)

# ----WRS2 package----
# trimmed means and post hoc tests
WRS2::t1way(happiness ~ dose, data = puppy_tib, nboot = 1000) |> 
  model_parameters() |> 
  display()

WRS2::lincon(happiness ~ dose, data = puppy_tib) |>  
  model_parameters() |>
  format_table() |> 
  display()
  

# bootstrap method
WRS2::t1waybt(happiness ~ dose, data = puppy_tib, nboot = 1000)
WRS2::mcppb20(happiness ~ dose, data = puppy_tib, nboot = 1000)
```


### Self-test

```{r}
#| label: wrs_st
#|

# trimmed means and post hoc tests
WRS2::t1way(happiness ~ dose, data = puppy_tib, tr = 0.1, nboot = 1000) |> 
model_parameters() |> 
  display()

WRS2::lincon(happiness ~ dose, data = puppy_tib, tr = 0.1) |> 
  model_parameters() |> 
  display()
```

## Bayesian model

```{r}
#| label: bfs


# Bayesian model with default priors
puppy_dp <- stan_glm(happiness ~ dose, data = puppy_tib)
# get summary and bayes factors
model_parameters(puppy_dp, test = c("pd", "bf"), null = 0) |> 
  display()



x_tib <- tibble(x = c(-5, 10))
# prior for therapy vs no puppies
ggplot(x_tib, aes(x=x)) +
  geom_function(fun = dnorm, args = list(mean = 2, sd = 2)) +
  scale_x_continuous(breaks = seq(-5, 10, 1)) +
  geom_vline(xintercept = 0, colour = "red") +
  theme_minimal()

# prior for therapy vs no puppies
ggplot(x_tib, aes(x=x)) +
  geom_function(fun = dnorm, args = list(mean = 1, sd = 0.5)) +
  scale_x_continuous(breaks = seq(-5, 10, 1)) +
  geom_vline(xintercept = 0, colour = "red") +
  theme_minimal()

# prior for intercept
ggplot(x_tib, aes(x=x)) +
  geom_function(fun = dnorm, args = list(mean = 4, sd = 1)) +
  scale_x_continuous(breaks = seq(-2, 10, 1)) +
  geom_vline(xintercept = 0, colour = "red") +
  theme_minimal()


puppy_rtb <- stan_glm(happiness ~ dose,
                      data = puppy_tib,
                      prior = normal(location = c(2, 1), scale = c(2, 0.5)),
                      prior_intercept = normal(location = 4, scale = 1))

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


