---
title: "Chapter 9: Categorical predictors"
author: "Andy Field"
date: "2025-04-14"
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)

# create data
cloak_tib <- tibble::tibble(
    id = c("Kiara", "Anupama", "Collin", "Ian", "Vanessa", "Darrell", "Tyler", "Steven", "Katheryn", "Kirkpatrick", "Melissa", "Kinaana", "Shajee'a", "Sage", "Jorge", "Conan","Tamara","Man", "Joseph", "Israa", "Kathryn", "Luis", "Devante", "Jerry"),
    cloak = gl(2, 12, labels = c("No cloak", "Cloak")),
    mischief = c(3, 1, 5, 4, 6, 4, 6, 2, 0, 5, 4, 5, 4, 3, 6, 6, 8, 5, 5, 4, 2, 5, 7, 5)
  )

# check levels
levels(cloak_tib$cloak)
```


## Descriptive statistics

```{r}
#| label: datawizardry

cloak_tib |> 
  group_by(cloak) |> 
  describe_distribution(ci = 0.95, iterations = 1000) |> 
  display()
```


## Fittting a linear model

```{r}
#| label: cloak_lm

# fit the model
cloak_lm <- lm(mischief ~ cloak, data = cloak_tib)
# extract the parameter estimates to render nicely
parameters(cloak_lm) |> 
  display()
```



## Comparing means without the raw scores

```{r}
#| label: compare_means_no_data

# a function to compare means

t_from_means <- function(m1, m2, sd1, sd2, n1, n2, p_digits = 3, digits = 2){
  df <- n1 + n2 - 2
	poolvar <- (((n1-1)*sd1^2)+((n2-1)*sd2^2))/df
	t <- (m1-m2)/sqrt(poolvar*((1/n1)+(1/n2)))
	sig <- 2*(1-(pt(abs(t),df)))
	paste0("t(", df, ") = ", round(t, digits), ", p = ", round(sig, p_digits))
}

# Try it out

t_from_means(m1 = 5, m2 = 3.75, sd1 = 1.651, sd2 = 1.913, n1 = 12, n2 = 12)
```

## Exploring data and testing assumptions

```{r}
#| label: inv_plot

ggplot(cloak_tib, aes(cloak, mischief)) +
  geom_violin() +
  stat_summary(fun.data = "mean_cl_normal") +
  labs(x = "Cloak group", y = "Acts of mischief") +
  theme_minimal()
```

## Independent *t*-test

```{r}
#| label: indie_t

# regular t-test
cloak_mod <- t.test(mischief ~ cloak, data = cloak_tib)
model_parameters(cloak_mod) |> 
  display(column_names = c("Outcome", "Predictor", "No cloak", "Cloak", "Mean difference", "95% CI", "*t*", "*p*"))

# robust t-test based on trimmed means
cloak_yn  <- WRS2::yuen(mischief ~ cloak, data = cloak_tib)
model_parameters(cloak_yn) |> 
  display()

# robust t-test based on trimmed means and bootstrap
cloak_ybt <- WRS2::yuenbt(mischief ~ cloak, data = cloak_tib) 

# robust t-test based on M-estimator
cloak_pb <- WRS2::pb2gen(mischief ~ cloak, data = cloak_tib) 

# Effect size using the control group SD
glass_delta(mischief ~ cloak, data = cloak_tib, reference = "No cloak", adjust = FALSE) |>
  display()

# re-do but adjusting for small sample bias
glass_delta(mischief ~ cloak, data = cloak_tib, reference = "No cloak", adjust = TRUE) |>
  display()

# Effect size using the pooled SD

cohens_d(mischief ~ cloak, data = cloak_tib, adjust = TRUE) |> 
  display()

# Bayes factor

cloak_bf <- BayesFactor::ttestBF(formula = mischief ~ cloak, data = cloak_tib, rscale = "medium")
model_parameters(cloak_bf) |> 
  display()
```

## Repeated measures designs

```{r}
#| label: depen_t

# enter the data

cloak_rm_tib <- tibble::tibble(
    id = rep(c("Alia", " Elisia", "Dominika", "Noah", "Haashid", "Zainab", "Reyna", "Victoria", "Najaat", "Jordyn", "Zaki", "Rajohn"), 2),
    cloak = gl(2, 12, labels = c("No cloak", "Cloak")),
    mischief = c(3, 1, 5, 4, 6, 4, 6, 2, 0, 5, 4, 5, 4, 3, 6, 6, 8, 5, 5, 4, 2, 5, 7, 5)
  )

# make it wide and compute the difference between scores

cloak_messy_tib <- cloak_rm_tib |>
  tidyr::pivot_wider(
    id_cols = id,
    values_from = mischief,
    names_from = cloak
    ) |>
  rename_with(.data = _, ~ tolower(str_replace(.x, pattern = " ", replacement = "_"))) |> 
  mutate(
    diff = cloak - no_cloak
  )

# qqplot

ggplot(cloak_messy_tib, aes(sample = diff)) +
    qqplotr::stat_qq_band(fill = "#5c97bf", alpha = 0.3) +
    qqplotr::stat_qq_line(colour = "#5c97bf") +
    qqplotr::stat_qq_point() +
    labs(x = "Theoretical quantiles", y = "Sample quantiles") +
    theme_minimal()


# the t-test
cloak_rm_mod  <- t.test(Pair(cloak, no_cloak) ~ 1, data = cloak_messy_tib)
model_parameters(cloak_rm_mod) |> 
  display()


# robust t-test
cloak_rm_rob  <- WRS2::yuend(x = cloak_messy_tib$cloak, y = cloak_messy_tib$no_cloak)
model_parameters(cloak_rm_rob) |> 
  display()

# effect size that factors in dependency
repeated_measures_d(Pair(cloak, no_cloak) ~ 1, data = cloak_messy_tib, method = "rm") |> 
  display()
```

## Bayesian models


```{r}
#| label: bayes


# prior for intercept 
ggplot(tibble(x = c(-1, 5)), aes(x=x)) +
  geom_function(fun = dnorm, args = list(mean = 2, sd = 0.7), colour = "#88CCEE", linewidth = 1) +
  scale_x_continuous(breaks = seq(-1, 5, 1)) +
  geom_vline(xintercept = 0, colour = "#999933", linetype = 2, linewidth = 1) +
  geom_vline(xintercept = 2, colour = "#CC6677", linewidth = 1) +
  ggtitle('Prior for intercept') +
  theme_minimal()

# prior for cloak 
ggplot(tibble(x = c(-3, 7)), aes(x=x)) +
  geom_function(fun = dnorm, args = list(mean = 2, sd = 1.2), colour = "#88CCEE", linewidth = 1) +
  scale_x_continuous(breaks = seq(-3, 7, 1)) +
  geom_vline(xintercept = 0, colour = "#999933", linetype = 2, linewidth = 1) +
  geom_vline(xintercept = 2, colour = "#CC6677", linewidth = 1) +
  ggtitle('Prior for cloak') +
  theme_minimal()


# Bayesian model with informative priors and bayes factors

cloak_rstn <- stan_glm(mischief ~ cloak,
                      data = cloak_tib,                      
                      prior = normal(location = 2, scale = 1.2),
                      prior_intercept = normal(location = 2, scale = 0.7))

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


