---
title: "R code in Chapter 7"
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(ggplot2)
library(magrittr)
```

***
## 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 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}
cats_tib <- discovr::roaming_cats
exam_tib <- discovr::exam_anxiety
liar_tib <- discovr::biggest_liar
```

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 `exam_anxiety` data you would load it by executing:

```{r, eval = FALSE}
exam_tib <- here::here("data/exam_anxiety.csv") %>%
  readr::read_csv() %>%
  dplyr::mutate(
    sex = forcats::as_factor(sex)
  )
```

## Bivariate correlations

Plot the data

```{r}
GGally::ggscatmat(exam_tib, columns = c("exam_grade", "revise", "anxiety")) +
  theme_minimal()
```

## Pearson's r

A single correlation:

```{r}
exam_tib %>% 
  dplyr::select(exam_grade, revise) %>% 
  correlation::correlation()
```

Let's scale this up to include exam anxiety:

```{r}
exam_cor <- exam_tib %>% 
  dplyr::select(exam_grade, revise, anxiety) %>% 
  correlation::correlation()

exam_cor
```

### Change the number of decimal places

```{r}
exam_tib %>% 
  dplyr::select(exam_grade, revise, anxiety) %>% 
  correlation::correlation(digits = 3, ci_digits = 3)
```

### Matrix format

```{r}
exam_tib %>% 
  dplyr::select(exam_grade, revise, anxiety) %>% 
  correlation::correlation() %>% 
  summary()
```



### Coefficient of determination

```{r}
exam_cor <- exam_tib %>% 
  dplyr::select(exam_grade, revise, anxiety) %>% 
  correlation::correlation()

(exam_cor$r)^2
```

## Robust correlations

```{r}
exam_tib %>% 
  dplyr::select(exam_grade, revise, anxiety) %>% 
  WRS2::winall()
```

```{r}
exam_tib %>% 
  dplyr::select(exam_grade, revise, anxiety) %>% 
  correlation::correlation(., method = "percentage")
```


## Spearman's rho

Load the Biggest Liar data from the CSV file (assuming you've set up a project as described in chapter 1)

```{r eval = FALSE}
liar_tib <- readr::read_csv("../data/biggest_liar.csv") %>% 
  dplyr::mutate(
    novice = forcats::as_factor(novice)
  )
```

```{r}
liar_tib %>%
  dplyr::select(position, creativity) %>% 
  correlation::correlation(method = "spearman")
```

## Kendall's tau

```{r}
liar_tib %>% 
  dplyr::select(position, creativity) %>% 
  correlation::correlation(method = "kendall")
```

### Bootstrapped confidence intervals

A quick aside, let's look at the code for writing a function to print the mean:

```{r}
print_mean <- function(variable){
	mean <- sum(variable)/nrow(variable)
       cat("Mean = ", mean)
}

exam_tib %>% 
  dplyr::select(exam_grade) %>%
  print_mean()

```

Or a sillier version:

```{r}
print_mean <- function(harry_the_hungy_hippo){
	mean <- sum(harry_the_hungy_hippo)/ nrow(harry_the_hungy_hippo)
       cat("Harry the hungry hippo say that the mean = ", mean)
}

exam_tib %>% 
  dplyr::select(exam_grade) %>%
  print_mean()

```

Here's a function to bootstrap Pearson's r:

```{r}
boot_r <- function(data, i){
  cor(data[i, "exam_grade"], data[i, "revise"])
}

grade_revise_bs <- boot::boot(exam_tib, boot_r, R = 2000)

grade_revise_bs 
boot::boot.ci(grade_revise_bs)
```

### A general bootstrap function

```{r}
boot_r <- function(data, var1, var2, i){
  cor(data[i, var1], data[i, var2])
}

grade_revise_bs <- boot::boot(exam_tib, boot_r, var1 = "exam_grade", var2 = "revise", R = 2000)
grade_revise_bs
boot::boot.ci(grade_revise_bs)
```

```{r}
grade_anx_bs <- boot::boot(exam_tib, boot_r, var1 = "exam_grade", var2 = "anxiety", R = 2000)
grade_anx_bs
boot::boot.ci(grade_anx_bs)
```

```{r}
revise_anx_bs <- boot::boot(exam_tib, boot_r, var1 = "revise", var2 = "anxiety", R = 2000)
revise_anx_bs
boot::boot.ci(revise_anx_bs)

```

### Bootrapping Spearman's rho and Kendall's tau

```{r}
boot_r <- function(data, var1, var2, method = "pearson", i){
  cor(data[i, var1], data[i, var2], method = method)
}

liar_bs <- boot::boot(liar_tib, boot_r, var1 = "position", var2 = "creativity", method = "spearman", R = 2000)
boot::boot.ci(liar_bs)

```

## Point-biserial correlation

If you need to read the data in from the csv file:

```{r, eval = FALSE}
cats_tib <- readr::read_csv("../data/roaming_cats.csv") %>% 
  dplyr::mutate(
    sex = forcats::as_factor(sex)
  )
```

```{r}
cats_tib <- cats_tib %>% 
  dplyr::mutate(
    sex_bin = ifelse(sex == "Male", 0, 1),
    sex_bin_recode = ifelse(sex == "Male", 1, 0)
  )
```

```{r}
cats_tib %>% 
  dplyr::select(time, sex_bin) %>%
  correlation::correlation()
```

coefficient of determination
```{r}
r_cat <- cats_tib %>% 
  dplyr::select(time, sex_bin) %>%
  correlation::correlation()

(r_cat$r)^2
```



```{r}
cats_tib %>% 
  dplyr::select(time, sex_bin_recode) %>%
  correlation::correlation()
```

```{r}
cats_tib %>% 
  dplyr::select(time, sex_bin) %>%
  correlation::correlation(method = "biserial")
```



## Partial correlations

```{r}
exam_tib %>% 
  dplyr::select(exam_grade, revise, anxiety) %>% 
  correlation::correlation(., partial = TRUE)
```

## Comparing rs

The correlation between exam grade and anxiety for males:

```{r}
exam_tib %>% 
  dplyr::filter(sex == "Male") %>% 
  dplyr::select(exam_grade, anxiety) %>% 
  correlation::correlation()
```

The correlation between exam grade and anxiety for males:

```{r}
exam_tib %>% 
  dplyr::filter(sex == "Female") %>% 
  dplyr::select(exam_grade, anxiety) %>% 
  correlation::correlation()
```

```{r}
exam_men <- exam_tib %>% 
  dplyr::filter(sex == "Male")

exam_women <- exam_tib %>% 
  dplyr::filter(sex == "Female")

WRS2::twopcor(x1 = exam_men$exam_grade, y1 = exam_men$anxiety, x2 = exam_women$exam_grade, y2 = exam_women$anxiety)
```

