---
title: "Chapter 6: The Beast of bias"
author: "Andy Field"
date: "2025-03-27"
date-modified: today
format: 
  html:
    theme: journal
    embed-resources: true
    toc: true
execute: 
    message: false
    warning: false
---

```{r}
#| label: setup

library(easystats)
library(tidyverse)

download_tib <- here::here("data/download_festival.csv") |> 
  read_csv() |>  
  mutate(
    ticket_no = as.character(ticket_no),
    gender = as_factor(gender) |> fct_relevel("Male", "Female", "Non-binary")
  )

# download_tib <- discovr::download

```

## Restructuring data

### Tidying messy data

```{r}
#| label: tdy_data

# Turn the messy data to tidy
download_tidy_tib <- download_tib |> 
  pivot_longer(
  cols = day_1:day_3,
  names_to = "day",
  values_to = "hygiene"
)

# print the data
download_tidy_tib

# tidy up the labels  for day
download_tidy_tib <- download_tidy_tib |> 
   mutate(
    day = str_to_sentence(day) |> str_replace("_", " ")
  )

# print the data
download_tidy_tib


# Make the data messy again

# Turn the messy data to tidy
download_tib <- download_tidy_tib |> 
  mutate(
    day = str_to_lower(day) |> str_replace(" ", "_")
  ) |> 
  pivot_wider(
    id_cols = c(ticket_no, gender),
    names_from = "day",
    values_from = "hygiene"
    )


```

## Spotting outliers

```{r}
#| label: plot_outliers

#Scatterplot

ggplot(download_tib, aes(x = day_1, y = day_2)) +
  geom_point(fill = "#882255", colour = "#CC6677", alpha = 0.7) +
  scale_y_continuous(breaks = seq(0, 5, 1)) +
  scale_x_continuous(breaks = seq(0, 20, 2)) +
  labs(x = "Hygiene scores, 0-5 (day 1)", y = "Hygiene scores, 0-5 (day 2)") +
  theme_minimal()

# Histogram
ggplot(download_tib, aes(day_1)) +
  geom_histogram(binwidth = 0.2, fill = "#882255", colour = "#CC6677", alpha = 0.2) +
  labs(y = "Frequency", x = "Hygiene scores (0-5)") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

# Boxplot
ggplot(download_tib, aes(x = "Day 1", y = day_1)) +
  geom_boxplot(fill = "#CC6677", alpha = 0.7) +
  scale_y_continuous(breaks = seq(0, 20, 2)) +
  labs(x = "Day of festival", y = "Hygiene scores (0-5)") +
  theme_minimal()

```

```{r}
#| label: catch_the_pigeon

# filter the data to find the unusual score
download_tib |> 
  filter(day_1 > 4)

# Correct the mis-typed score

download_tib <- download_tib |>
  mutate(
    day_1 = recode_values(day_1, list(`2.02` = 20.02))
  )

# recreate the tidy version (otherwise it will also have the error)

download_tidy_tib <- download_tib |> 
  pivot_longer(
    cols = day_1:day_3,
    names_to = "day",
    values_to = "hygiene") |> 
   mutate(
    day = str_to_sentence(day) |> str_replace("_", " ")
    )

```

### Using z-scores

```{r}
#| label: z_scores

 
zdownload_tib <- standardise(x = download_tib,
                          select = day_1:day_3)

# list z-scores greater or equal to 1.96 on day 1
zdownload_tib |> 
  filter(abs(day_1) >= 1.96) |> 
  arrange(day_1)

# list z-scores greater or equal to 2.58 on any day
zdownload_tib |> 
  filter(
    if_any(day_1:day_3, ~ . >= 2.58)
    )
```

```{r}
#| label: z_prop

# write our own function to count the proportion of z-scores
 get_cum_percent <- function(var,  cut_off = 1.96){
  ecdf_var <- abs(var) |> ecdf()
  100*(1 - ecdf_var(cut_off))
 }
 
# use the function
 get_cum_percent(zdownload_tib$day_1)
 
# use other thresholds

zdownload_tib |> 
  summarize(
    z_1.96 = get_cum_percent(day_1),
    z_2.58 = get_cum_percent(day_1,cut_off = 2.58),
    z_3.29 = get_cum_percent(day_1,cut_off = 3.29)
  )

 
# apply to multiple variables and use other thresholds

zdownload_tidy_tib <- zdownload_tib |> 
  pivot_longer(
    cols = day_1:day_3,
    names_to = "day",
    values_to = 'zhygiene'
  ) 

zdownload_tidy_tib |> 
  group_by(day) |> 
  summarize(
    z_1.96 = get_cum_percent(zhygiene),
    z_2.58 = get_cum_percent(zhygiene,cut_off = 2.58),
    z_3.29 = get_cum_percent(zhygiene,cut_off = 3.29)
  )



```

## qqplots

```{r}
#| label: qqplots

# create a qqplot for each day
ggplot(download_tidy_tib, aes(sample = hygiene)) +
  qqplotr::stat_qq_band(fill = "#DDCC77", alpha = 0.3) +
  qqplotr::stat_qq_line(colour = "#DDCC77") +
  qqplotr::stat_qq_point(colour = "#CC6677", alpha = 0.2, size = 1) +
  labs(x = "Theoretical quantiles", y = "Sample quantiles") +
  facet_wrap(~day, ncol = 1) +
  theme_minimal()

#create a detrended plot of day 2

ggplot(download_tib, aes(sample = day_2)) +
  qqplotr::stat_qq_band(fill = "#DDCC77", alpha = 0.3, detrend = TRUE) +
  qqplotr::stat_qq_line(colour = "#DDCC77", detrend = TRUE) +
  qqplotr::stat_qq_point(colour = "#CC6677", alpha = 0.2, size = 1, detrend = TRUE) +
  labs(x = "Theoretical quantiles", y = "Sample quantiles") +
  theme_minimal()


```

## Descriptive statistics

```{r}
#| label: datawizardry


download_tidy_tib |> 
  group_by(day) |> 
  describe_distribution(hygiene) |>
  display(digits = 2)
```

## Spotting hetroscedasticity

```{r}
#| label: boxplots

download_tidy_tib |> 
  filter(day == "Day 1") |> 
ggplot(aes(x = gender, y = hygiene, colour = gender, fill = gender)) +
  geom_boxplot(alpha = 0.5) +
  scale_y_continuous(breaks = seq(0, 4, 1)) +
  scale_colour_viridis_d() +
  scale_fill_viridis_d() +
  labs(x = "Gender of concert goer", y = "Hygiene scores (0-5)") +
  theme_minimal() +
  theme(legend.position = "none")


download_tidy_tib |> 
  filter(day == "Day 1") |> 
  group_by(gender) |> 
  summarize(
    variance = var(hygiene)
  )

```
