---
title: "Chapter 7: Correlation"
author: "Andy Field"
date: "2025-04-02"
date-modified: today
format: 
  html:
    theme: journal
    embed-resources: true
    toc: true
execute: 
    message: false
    warning: false
---

```{r}
#| label: setup

library(easystats)
library(tidyverse)

# load exam data
exam_tib <- here::here("data/exam_anxiety.csv") |> 
  read_csv() |> 
  mutate(
    sex = as_factor(sex)
  )

# load liar data

liar_tib <- here::here("data/biggest_liar.csv") |> 
  read_csv() |> 
  mutate(
    novice = as_factor(novice)
  )

# load cat data

cats_tib <- here::here("data/roaming_cats.csv") |> 
  read_csv() |> 
  mutate(
    sex = as_factor(sex)
  )

# or load from the dsicovr package
# exam_tib <- discovr::exam_anxiety
# liar_tib <- discovr::biggest_liar
# cats_tib <- discovr::roaming_cats

```

## Pearson's r

```{r}
#| label: exam_anxiety

#plot the data

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

# Pearson correlation for two variables

correlation(data = exam_tib,
            select = c("exam_grade", "revise"),
            rename = c("Exam grade (%)", "Time revising")) |> 
  display()

# Pearson correlation for multiple variables
correlation(data = exam_tib, select = c("exam_grade", "revise", "anxiety")) |> 
   display()


# output in grid format
exam_r <- correlation(data = exam_tib, select = c("exam_grade", "revise", "anxiety")) |> 
  summary()
display(exam_r)

# use across to calculate coefficients of determination
exam_r |> 
  mutate(
    across(where(is.numeric), \(x) x^2)
    ) |> 
  display(digits = 3)
```

## Robust correlations

```{r}
#| label: robr

# Winsorized correlations

correlation(data = exam_tib,
            select = c("exam_grade", "revise", "anxiety"),
            winsorize = 0.2) |> 
  display()

# Percentage bend correlation

correlation(data = exam_tib,
            select = c("exam_grade", "revise", "anxiety"),
            method = "percentage",
            beta = 0.2) |> 
  display()

```

## Spearman's rho

```{r}
#| label: spearman

correlation(data = liar_tib, select = c("position", "creativity"), method = "spearman") |> 
  display()


```

## Kendall's tau

```{r}
#| label: kendall

correlation(data = liar_tib, select = c("position", "creativity"), method = "kendall") |> 
  display()


```

## Biserial and point-biserial correlations

```{r}
#| label: biserial

# create a binary variable from sex (and also a reverse coded version for pedagogic reasons)
cats_tib <- cats_tib |> 
  mutate(
    sex_bin = ifelse(sex == "Male", 0, 1),
    sex_bin_recode = ifelse(sex == "Male", 1, 0)
  )

#point-biserial correlation
correlation(data = cats_tib, select = c("time", "sex_bin")) |> 
  display()

# biserial correlation
correlation(data = cats_tib,
            select = c("time", "sex_bin"),
            method = "biserial") |> 
  display()
```

## Partial correlation

```{r}
#| label: partialr

# partial correlation
correlation(data = exam_tib,
            select = c("exam_grade", "revise", "anxiety"),
            partial = TRUE) |> 
  display()
```

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

# semi-partial correlation (at the time of writing this doesn't work)
correlation(data = exam_tib,
            select = c("exam_grade", "revise", "anxiety"),
            partial = "semi") |> 
  display()

```

## Comparing correlations

```{r}
#| label: compare_r

# get the Pearson's r for males and females
exam_tib |> 
  group_by(sex) |> 
  correlation(select = c("exam_grade", "anxiety")) |> 
  display(digits = 3)


#split into male and female tibbles

exam_male <- exam_tib |>  
  filter(sex == "Male")

exam_female <- exam_tib |> 
  filter(sex == "Female")

# compare correlations

WRS2::twopcor(x1 = exam_male$exam_grade,
              y1 = exam_male$anxiety,
              x2 = exam_female$exam_grade,
              y2 = exam_female$anxiety)

```
