---
title: "R code in Chapter 6"
author: "Andy P. Field"
date: 'Last updated `r format(Sys.Date(), "%d %B %Y")`'
output: html_document
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  warning = FALSE,
  message = FALSE
)

library(magrittr)
library(dplyr)
library(ggplot2)
```

***
## Usage

This file contains code relevant to chapter 1 of [Discovering Statistics Using R and RStudio](https://www.discovr.rocks) by [Andy Field](https://www.discoveringstatistics.com/about). These files contain abridged sections from the book so there are some copyright considerations but I offer them under a [Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License](http://creativecommons.org/licenses/by-nc-nd/4.0/).^[Basically you can use this material for teaching and non-profit activities but do not meddle with it or claim it as your own work.]

***

## Load general packages

Remember to load the tidyverse:

```{r, eval = F}
library(tidyverse)
```

## Load the data

Remember to install the package for the book with `library(discovr)`. After which you can load data into a tibble by executing:

`name_of_tib <- discovr::name_of_data`

For example, execute these lines to create the tibbles referred to in the chapter:

```{r}
download_tib <- discovr::download
```

If you want to read the file from the CSV and you have set up your project folder as I suggest in Chapter 1, then the general format of the code you would use is:

```{r, eval = FALSE}
tibble_name <- here::here("../data/name_of_file.csv") %>%
  readr::read_csv() %>%
  dplyr::mutate(
    ...
    code to convert character variables to factors
    ...
  )
```

In which you'd replace `tibble_name` with the name you want to assign to the tibble and change `name_of_file.csv` to the name of the file that you're trying to load. You can use mutate to convert categorical variables to factors. For example, for the `download` data you would load it by executing:

```{r, eval = FALSE}
library(here)

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

## Self-test

The original data: 1, 3, 4, 3, 2. The mean and sum of squares are:

```{r}
friends <- c(1, 3, 4, 3, 2)
mean_old <- mean(friends)
ss_old <- (friends-mean_old)^2 %>% sum(.)

print(paste0("Original mean: ", mean_old))
print(paste0("Original SS: ", ss_old))
```

The data with the 4 replaced with a 10: 1, 3, 10, 3, and 2. The mean and sum of squares are:

```{r}
friends_out <- c(1, 3, 10, 3, 2)
mean_out <- mean(friends_out)
ss_out <- (friends_out-mean_out)^2 %>% sum(.)

print(paste0("Mean with outlier: ", mean_out))
print(paste0("SS with outlier: ", ss_out))
```

## plot histogram

```{r}
ggplot2::ggplot(download_tib, aes(day_1)) +
  geom_histogram(binwidth = 1, fill = "#56B4E9", colour = "#336c8b", alpha = 0.2) +
  labs(y = "Frequency", x = "Hygiene scores (0-5)", title = "Hygiene scores on day 1") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))
```

## Boxplot day 1

```{r}
ggplot2::ggplot(download_tib, aes(x = "Day 1", y = day_1)) +
  geom_boxplot(fill = "#5c97bf", alpha = 0.7) +
  scale_y_continuous(breaks = seq(0, 20, 2)) +
  labs(x = "Day of festival", y = "Hygiene scores (0-5)") +
  theme_minimal()
```

## Filter the data

```{r}
download_tib %>% 
  dplyr::filter(day_1 > 4)
```

## Replace the incorrect value

```{r}
download_tib <- download_tib %>%
  dplyr::mutate(
    day_1 = dplyr::recode(day_1, `20.02` = 2.02)
  )
```

## Replot the histogram and boxplot

```{r}
ggplot2::ggplot(download_tib, aes(day_1)) +
  geom_histogram(binwidth = 1, fill = "#56B4E9", colour = "#336c8b", alpha = 0.2) +
  labs(y = "Frequency", x = "Hygiene scores (0-5)", title = "Hygiene scores on day 1") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

ggplot2::ggplot(download_tib, aes(x = "Day 1", y = day_1)) +
  geom_boxplot(fill = "#5c97bf", alpha = 0.7) +
  scale_y_continuous(breaks = seq(0, 4, 1)) +
  labs(x = "Day of festival", y = "Hygiene scores (0-5)") +
  theme_minimal()
```

## Tranforming between messy and tidy data formats

```{r}
download_tidy_tib <- download_tib %>% 
  tidyr::pivot_longer(
  cols = day_1:day_3,
  names_to = "day",
  values_to = "hygiene",
) 
```

```{r, eval = FALSE}
# If you want to tidy up the labels of day

download_tidy_tib <- download_tidy_tib %>% 
  dplyr::mutate(
    day = stringr::str_to_sentence(day) %>%
      stringr::str_replace(., "_", " ")
  )
```



Converting back to messy format:

```{r}
download_tib <- download_tidy_tib %>% 
  tidyr::pivot_wider(
  id_cols = c(ticket_no, gender),
  names_from = "day",
  values_from = "hygiene",
) 
```

splitting multiple variables

```{r}

# Generate some random 'messy' format data

set.seed(666)
wetwipe_tib <- tibble::tibble(
  id = 1:10,
  wetwipe_day1 = round(rnorm(10, 1.7, 0.6), 2),
  wetwipe_day2 = round(rnorm(10, 1.5, 0.6), 2),
  wetwipe_day3 = round(rnorm(10, 1.5, 0.6), 2),
  none_day1 = round(rnorm(10, 1.7, 0.6), 2),
  none_day2 = round(rnorm(10, 0.8, 0.6), 2),
  none_day3 = round(rnorm(10, 0.8, 0.6), 2)
)

# Tidy that shit up ..

wetwipe_tidy_tib <- wetwipe_tib %>% 
  tidyr::pivot_longer(
    cols = wetwipe_day1:none_day3,
    names_to = c("wetwipe", "day"),
    names_sep = "_",
    values_to = "hygiene"
  )

# Reverse that shit ....

wetwipe_tib <- wetwipe_tidy_tib %>% 
  tidyr::pivot_wider(
    id_cols = id,
    names_from = c("wetwipe", "day"),
    names_sep = "_",
    values_from = "hygiene"
  )
```

## Boxplots for each day split by gender

```{r}
ggplot2::ggplot(download_tidy_tib, aes(day, hygiene, fill = gender)) +
  geom_boxplot(alpha = 0.7) +
  scale_y_continuous(breaks = seq(0, 4, 1)) +
  labs(x = "Day of festival", y = "Hygiene scores (0-5)", fill = "Gender") +
  facet_wrap(~ gender) +
  theme_minimal()
```


## Using z-scores

Manually convert to z:

```{r}

make_z <- function(x){
  (x - mean(x, na.rm = TRUE)) / sd(x, na.rm = TRUE)
}

download_tib <- download_tib %>% 
  dplyr::mutate(
    day_1_z = make_z(day_1),
    day_2_z = make_z(day_2),
    day_3_z = make_z(day_3)
  )

download_tib
```

Using `mutate_at()`:

```{r}
download_tib <- download_tib %>% 
  dplyr::mutate_at(
    vars(day_1:day_3),
    list(z = make_z)
  )
```

Filter a variable at a time

```{r}
download_tib %>% 
  dplyr::filter(abs(day_1_z) >= 1.96) %>% 
  dplyr::arrange(day_1_z)
```

Filter by all variables: 

```{r}
download_tib %>% 
  dplyr::filter_at(
    vars(day_1_z:day_3_z), any_vars(. >= 2.58)
    ) %>% 
  dplyr::select(-c(day_1:day_3))
```



```{r}
get_cum_percent <- function(var,  cut_off = 1.96){
  ecdf_var <- abs(var) %>% ecdf()
  100*(1 - ecdf_var(cut_off))
}

download_tib %>%
  dplyr::summarize(
    `z >= 1.96` = get_cum_percent(day_1_z),
    `z >= 2.58` = get_cum_percent(day_1_z,cut_off = 2.58),
    `z >= 3.29` = get_cum_percent(day_1_z,cut_off = 3.29)
  )
```


```{r}
download_tidy_tib <- download_tidy_tib %>% 
  group_by(day) %>% 
  dplyr::mutate(
    zhygiene = make_z(hygiene)
  )


download_tidy_tib %>% 
  group_by(day) %>%
  dplyr::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)
  )
```

## Q-Q plots

Standard Q-Q plot of the hygiene scores

```{r}
download_tidy_tib %>%
  ggplot2::ggplot(., aes(sample = hygiene)) +
  qqplotr::stat_qq_band(fill = "#5c97bf", alpha = 0.3) +
  qqplotr::stat_qq_line(colour = "#5c97bf") +
  qqplotr::stat_qq_point(alpha = 0.2, size = 1) +
  labs(x = "Theoretical quantiles", y = "Sample quantiles") +
  facet_wrap(~day, ncol = 1) +
  theme_minimal()
```

Detrended Q-Q plot

```{r}
download_tidy_tib %>%
  dplyr::filter(day == "Day 2") %>% 
  ggplot2::ggplot(., aes(sample = hygiene)) +
  qqplotr::stat_qq_band(fill = "#5c97bf", alpha = 0.3, detrend = TRUE) +
  qqplotr::stat_qq_line(colour = "#5c97bf", detrend = TRUE) +
  qqplotr::stat_qq_point(alpha = 0.2, size = 1, detrend = TRUE) +
  labs(x = "Theoretical quantiles", y = "Sample quantiles") +
  theme_minimal()
```


## Skew and kurtosis

```{r}
download_tidy_tib %>% 
  dplyr::group_by(day) %>% 
  dplyr::summarize(
    valid_cases = sum(!is.na(hygiene)),
    mean = ggplot2::mean_cl_normal(hygiene)$y,
    ci_lower = ggplot2::mean_cl_normal(hygiene)$ymin,
    ci_upper = ggplot2::mean_cl_normal(hygiene)$ymax,
    skew = moments::skewness(hygiene, na.rm = TRUE),
    kurtosis = moments::kurtosis(hygiene, na.rm = TRUE)
  ) 
```

## Splitting by gender

```{r}
download_tidy_tib %>% 
  dplyr::group_by(gender, day) %>% 
  dplyr::summarize(
    valid_cases = sum(!is.na(hygiene)),
    mean = ggplot2::mean_cl_normal(hygiene)$y,
    ci_lower = ggplot2::mean_cl_normal(hygiene)$ymin,
    ci_upper = ggplot2::mean_cl_normal(hygiene)$ymax,
    skew = moments::skewness(hygiene, na.rm = TRUE),
    kurtosis = moments::kurtosis(hygiene, na.rm = TRUE)
  )

```

```{r}
download_tidy_tib %>%
  ggplot2::ggplot(., aes(sample = hygiene)) +
  qqplotr::stat_qq_band(fill = "#5c97bf", alpha = 0.3) +
  qqplotr::stat_qq_line(colour = "#5c97bf") +
  qqplotr::stat_qq_point(alpha = 0.2, size = 1) +
  labs(x = "Theoretical quantiles", y = "Sample quantiles") +
  facet_grid(day~gender, scales = "free_x") +
  theme_minimal()
```


