---
title: "Chapter 5: Describing and visualising data"
author: "Andy Field"
date: "2025-03-20"
date-modified: today
format: 
  html:
    theme: journal
    embed-resources: true
    toc: true
    code-overflow: wrap
execute: 
    message: false
    warning: false
---

```{r}
#| label: setup

library(easystats)
library(tidyverse)

eddiefy <- discovr::eddiefy

```

## Describing data

### Central tendency

```{r}
#| label: median

# Create data
insta_tib <- tibble(
    followers = c(57, 40, 103, 234, 93, 53, 116, 98, 108, 121, 22)
   ) 

# Print the median as a single value
median(insta_tib$followers)
```

```{r}
#| label: mean

# Print the mean with no trimming
mean(insta_tib$followers)

# Print the 10% trimmed mean 
mean(insta_tib$followers, trim = 0.1)
```

```{r}
#| label: ct_table

insta_tib |>
  summarize(
    median =  median(followers),
    mean =  mean(followers),
    trimmed_mean =  mean(followers, trim = 0.1)
  )
```

### Dispersion

```{r}
#| label: dispersion

# Compute range
max(insta_tib$followers) - min(insta_tib$followers)

# Compute quartiles
quantile(insta_tib$followers, probs = c(0.25, 0.5, 0.75))

# Compute lowqer quartile
quantile(insta_tib$followers, probs = 0.25)

# Compute lowqer quartile
quantile(insta_tib$followers, probs = 0.75)

# Compute the interquartile range:
IQR(insta_tib$followers)

# variance and standard deviation by executing:
var(insta_tib$followers)
sd(insta_tib$followers)
```

Put everything into a table

```{r}
#| label: summary_table

insta_tib |>
  summarize(
    median =  median(followers),
    mean =  mean(followers),
    trimmed_mean =  mean(followers, trim = 0.1),
    range = max(followers) - min(followers),
    lower_quartile = quantile(followers, probs = 0.25),
    upper_quartile = quantile(followers, probs = 0.75),
    IQR = IQR(followers),
    var = var(followers),
    sd = sd(followers)
    )

```

Create a nicer table for rendering using `display() as shown in`@tbl-summary.

```{r}
#| label: tbl-summary
#| tbl-cap: Summary statistics for instagram data

insta_tbl <- insta_tib |>
  summarize(
    median =  median(followers),
    mean =  mean(followers),
    trimmed_mean =  mean(followers, trim = 0.1),
    range = max(followers) - min(followers),
    lower_quartile = quantile(followers, probs = 0.25),
    upper_quartile = quantile(followers, probs = 0.75),
    IQR = IQR(followers),
    var = var(followers),
    sd = sd(followers)
    )

display(insta_tbl,
        digits = 2,
        column_names = c("Median", "Mean", "5% trimmed Mean", "Range", "*Q~L~*", "*Q~U~*", "IQR", "*s*^2^", "*s*"))
```

### Shape

We can use the `datawizard` package to compute skew and kurtosis. If we have loaded `easystats` this package will be available to us and we can directly use the functions `skewness()` and `kurtosis()`.

```{r}
#| label: shape

# Skewness
skewness(insta_tib$followers)

# Kurtosis
kurtosis(insta_tib$followers)
```

A detour into classes:

```{r}
#| label: demo_class

#create a skewness object called 'insta_skew'
insta_skew <- skewness(insta_tib$followers)
insta_skew$Skewness
insta_skew$SE

```

Adding skew and kurtosis to a table:

```{r}
#| label: skew_tbl


insta_tib |>
  summarize( 
    skew = skewness(followers)$Skewness,
    kurtosis = kurtosis(followers)$Kurtosis
    )
```

### Confidence intervals

```{r}
#| label: ci

# get a normal 95% CI for a mean
mean_cl_normal(insta_tib$followers)

# get a normal 90% CI for a mean
mean_cl_normal(insta_tib$followers, conf.int = 0.90)

# extract the lower bound of a 95% CI
mean_cl_normal(insta_tib$followers)$ymin

```

Tabulating the confidence interval

```{r}
#| label: ci_tbl

insta_tib |>
  summarize(
    mean =  mean_cl_normal(followers)$y,
    ci_lower = mean_cl_normal(followers)$ymin,
    ci_upper = mean_cl_normal(followers)$ymax
    )

```

### Descriptives with datawizard

```{r}
#| label: datawizardry

describe_distribution(x = insta_tib, ci = 0.95, iterations = 500)

```

### Nicer tables with `display()`

To get nicer rendered tables use `display()`

```{r}
#| label: no_display

describe_distribution(x = insta_tib)
```

```{r}
#| label: display

describe_distribution(x = insta_tib) |> 
  display(column_names = c("Variable", "*M*", "*s*", "IQR", "Range", "Skewness", "Kurtosis", "*n*", "Missing"))

```

## ggplot2

### An example to illustrate key ideas

```{r}
#| label: eddiedat

eddiefy_mini <- eddiefy |> 
  filter(album_name %in% c("Killers", "Piece of Mind", "Powerslave", "The Final Frontier"))

```

@fig-eddiegg shows the distribution and mean Spotify valence scores for four Iron Maiden albums.

```{r}
#| label: fig-eddiegg
#| fig-cap: Valence scores for four Iron Maiden albums


ggplot(eddiefy_mini, aes(x = album_name, y = valence, fill = album_name)) +
  theme_minimal() +
  geom_violin() +
  stat_summary(fun.data = "mean_cl_normal") +
  labs(x = "Album name", y = "Valence (0-1)") +
  scale_y_continuous(breaks = seq(0, 0.8, 0.1))

```

### geoms and stats

```{r}
#| label: fig-eddiehist
#| fig-cap: Histogram of Iron Maiden valence scores

ggplot(eddiefy, aes(x = valence)) +
  geom_histogram()
```

### Adding a theme

```{r}
#| label: fig-eddiedarkhist
#| fig-cap: Histogram of Iron Maiden valence scores

ggplot(eddiefy, aes(x = valence)) +
  geom_histogram() +
  theme_dark()

```

### 

## Histograms

```{r}
#| label: setup_for_histogram

wish_tib <- here::here("data/jiminy_cricket.csv") |> 
  read_csv() |> 
  mutate(
    strategy = as_factor(strategy),
    time = as_factor(time) |> fct_relevel("Baseline")
  )


# wish_tib <- discovr::jiminy_cricket
```

```{r}
#| label: exploring_histograms

# A default histogram
ggplot(wish_tib, aes(x = success)) +
  geom_histogram()

# change the binwidth
ggplot(wish_tib, aes(x = success)) +
  geom_histogram(binwidth = 5)

# change the bars to blue defined by hex code "#56B4E9":
ggplot(wish_tib, aes(x = success)) +
  geom_histogram(binwidth = 5, fill = "#56B4E9")

# change the colour of the lines around the bars using hex code #336c8b:
ggplot(wish_tib, aes(x = success)) +
  geom_histogram(binwidth = 5, fill = "#56B4E9", colour = "#336c8b")

# change opacity: histogram bars with 20% opacity:
ggplot(wish_tib, aes(x = success)) +
  geom_histogram(binwidth = 5, fill = "#56B4E9", colour = "#336c8b", alpha = 0.2)

# replace axis labels

ggplot(wish_tib, aes(x = success)) +
  geom_histogram(binwidth = 5, fill = "#56B4E9", colour = "#336c8b", alpha = 0.2) +
  labs(y = "Frequency", x = "Success scores (0-100)")

# apply a different theme
ggplot(wish_tib, aes(x = success)) +
 geom_histogram(binwidth = 5, fill = "#56B4E9", colour = "#336c8b", alpha = 0.2) +
 labs(y = "Frequency", x = "Success scores (0-100)") +
 theme_minimal()

# apply facets
ggplot(wish_tib, aes(x = success)) +
 geom_histogram(binwidth = 5, fill = "#56B4E9", colour = "#336c8b", alpha = 0.2) +
 labs(y = "Frequency", x = "Success scores (0-100)") +
 facet_grid(strategy~time) +
 theme_minimal()



```

## Frequency polygons

```{r}
#| label: freq_poly

ggplot(wish_tib, aes(x = success)) +
 geom_freqpoly(binwidth = 5, colour = "#336c8b") +
 labs(y = "Frequency", x = "Success scores (0-100)") +
 facet_grid(strategy~time) +
 theme_minimal()

# adjust the line style
ggplot(wish_tib, aes(x = success)) +
 geom_freqpoly(binwidth = 5, colour = "#336c8b", linetype = 2) +
 labs(y = "Frequency", x = "Success scores (0-100)") +
 facet_grid(strategy~time) +
 theme_minimal()

# adjust the line thickness
ggplot(wish_tib, aes(x = success)) +
 geom_freqpoly(binwidth = 5, colour = "#336c8b", size = 2) +
 labs(y = "Frequency", x = "Success scores (0-100)") +
 facet_grid(strategy~time) +
 theme_minimal()

```

## Boxplots

```{r}
#| label: boxplots

# a basic boxplot
ggplot(wish_tib, aes(x = time, y = success)) +
  geom_boxplot()

# let's funk it up a bit
ggplot(wish_tib, aes(x = time, y = success)) +
  geom_boxplot(fill = "#5c97bf", alpha = 0.7) +
  scale_y_continuous(breaks = seq(0, 100, 10)) +
  labs(x = "Time", y = "Success (out of 100)") +
  theme_minimal()

# using facets
ggplot(wish_tib, aes(x = time, y = success)) +
  geom_boxplot(fill = "#5c97bf", alpha = 0.7) +
  scale_y_continuous(breaks = seq(0, 100, 10)) +
  labs(x = "Time", y = "Success (out of 100)") +
  facet_wrap(~ strategy) +
  theme_minimal()

# let's add colour

ggplot(wish_tib, aes(x = time, y = success, fill = strategy)) +
  geom_boxplot(alpha = 0.7) +
  scale_y_continuous(breaks = seq(0, 100, 10)) +
  labs(x = "Time", y = "Success (out of 100)", fill = "Success strategy") +
  theme_minimal()


```

## Plotting means

```{r}
#| label: setup_for_means

notebook_tib <- here::here("data/notebook.csv") |> 
  read_csv() |> 
  mutate(
    gender_identity = as_factor(gender_identity),
    film = as_factor(film)
    )


# notebook_tib <- discovr::notebook
```

```{r}
#| label: mean_plots

# a basic plo of means
ggplot(notebook_tib, aes(x = film, y = arousal)) +
  stat_summary(fun = "mean", geom = "point", size = 3)

# add some labels and a theme
ggplot(notebook_tib, aes(x = film, y = arousal)) +
  stat_summary(fun = "mean", geom = "point", size = 3) +
  labs(x = "Type of film", y = "Arousal (out of 100)") +
  theme_minimal()

# rescale the y-axis and add senseible breaks for the y-axis ticks
ggplot(notebook_tib, aes(x = film, y = arousal)) +
  stat_summary(fun = "mean", geom = "point", size = 3) +
  coord_cartesian(ylim = c(0, 50)) +
  scale_y_continuous(breaks = seq(0, 50, 5)) +
  labs(x = "Type of film", y = "Arousal (out of 100)") +
  theme_minimal()

# include the 95% CI for each mean

ggplot(notebook_tib, aes(x = film, y = arousal)) +
  stat_summary(fun.data = "mean_cl_normal", geom = "pointrange") +
  coord_cartesian(ylim = c(0, 50)) +
  scale_y_continuous(breaks = seq(0, 50, 5)) +
  labs(x = "Type of film", y = "Arousal (out of 100)") +
  theme_minimal()


# include the robust 95% CI for each mean

ggplot(notebook_tib, aes(x = film, y = arousal)) +
  stat_summary(fun.data = "mean_cl_boot", geom = "pointrange") +
  coord_cartesian(ylim = c(0, 50)) +
  scale_y_continuous(breaks = seq(0, 50, 5)) +
  labs(x = "Type of film", y = "Arousal (out of 100)") +
  theme_minimal()

# split the plot by gender identity using facets
ggplot(notebook_tib, aes(x = film, y = arousal)) +
  stat_summary(fun.data = "mean_cl_normal", geom = "pointrange") +
  coord_cartesian(ylim = c(0, 50)) +
  scale_y_continuous(breaks = seq(0, 50, 5)) +
  labs(x = "Type of film", y = "Arousal (out of 100)") +
  facet_wrap(~gender_identity) +
  theme_minimal()

# split the plot by gender identity using colour
ggplot(notebook_tib, aes(x = film, y = arousal, colour = gender_identity)) +
  stat_summary(fun.data = "mean_cl_normal", geom = "pointrange") +
  coord_cartesian(ylim = c(0, 50)) +
  scale_y_continuous(breaks = seq(0, 50, 5)) +
  labs(x = "Type of film", y = "Arousal (out of 100)", colour = "Gender identity") +
  theme_minimal()

# use position dodge

ggplot(notebook_tib, aes(x = film, y = arousal, colour = gender_identity)) +
  stat_summary(fun.data = "mean_cl_normal", geom = "pointrange", position = position_dodge(width = 0.5)) +
  coord_cartesian(ylim = c(0, 50)) +
  scale_y_continuous(breaks = seq(0, 50, 5)) +
  labs(x = "Type of film", y = "Arousal (out of 100)", colour = "Gender identity") +
  theme_minimal()


```

## Violin plots

```{r}
#| label: violins

ggplot(notebook_tib, aes(x = film, y = arousal, fill = gender_identity)) +
  geom_violin() +
  stat_summary(fun.data = "mean_cl_normal", geom = "pointrange", position = position_dodge(width = 0.5)) +
  coord_cartesian(ylim = c(0, 50)) +
  scale_y_continuous(breaks = seq(0, 50, 5)) +
  labs(x = "Type of film", y = "Arousal (out of 100)", fill = "Gender identity") +
  theme_minimal()


# adjust the dodge

ggplot(notebook_tib, aes(x = film, y = arousal, fill = gender_identity)) +
  geom_violin(alpha = 0.7) +
  stat_summary(fun.data = "mean_cl_normal", geom = "pointrange", position = position_dodge(width = 0.9)) +
  coord_cartesian(ylim = c(0, 50)) +
  scale_y_continuous(breaks = seq(0, 50, 5)) +
  labs(x = "Type of film", y = "Arousal (out of 100)", fill = "Gender identity") +
  theme_minimal()

# add boxplots rather than pointranges

ggplot(notebook_tib, aes(x = film, y = arousal, fill = gender_identity)) +
  geom_violin() +
  geom_boxplot() +
  coord_cartesian(ylim = c(0, 50)) +
  scale_y_continuous(breaks = seq(0, 50, 5)) +
  labs(x = "Type of film", y = "Arousal (out of 100)", fill = "Gender identity") +
  theme_minimal()

#adjust the boxplots

ggplot(notebook_tib, aes(x = film, y = arousal, fill = gender_identity)) +
  geom_violin() +
  geom_boxplot(alpha = 0, width = 0.2, position = position_dodge(width = 0.9)) +
  coord_cartesian(ylim = c(0, 50)) +
  scale_y_continuous(breaks = seq(0, 50, 5)) +
  labs(x = "Type of film", y = "Arousal (out of 100)", fill = "Gender identity") +
  theme_minimal()

```

\

```{r}
#| label: order_matters

#Placing geom_violin after stat_summary()
ggplot(notebook_tib, aes(film, arousal, fill = gender_identity)) +
  stat_summary(fun.data = "mean_cl_normal", geom = "pointrange", position = position_dodge(width = 0.9)) +
  geom_violin() +
  labs(x = "Type of film", y = "Arousal (out of 100)", fill = "Gender identity") +
  theme_minimal()

# now place it before

ggplot(notebook_tib, aes(film, arousal, fill = gender_identity)) +
  geom_violin() +
  stat_summary(fun.data = "mean_cl_normal", geom = "pointrange", position = position_dodge(width = 0.9)) +
  labs(x = "Type of film", y = "Arousal (out of 100)", fill = "Gender identity") +
  theme_minimal()
```

## Repeated measures

```{r}
#| label: rm_designs

hiccups_tib <- here::here("data/hiccups.csv") |> 
  read_csv() |> 
  mutate(
    intervention = as_factor(intervention)
    )

# hiccups_tib <- discovr::hiccups

# view the data ordered by id


hiccups_tib |> 
  arrange(id)
```

```{r}
#| label: rm_plot

ggplot(hiccups_tib, aes(x = intervention, y = hiccups)) +
  stat_summary(fun.data = "mean_cl_normal", geom = "pointrange") +
  labs(x = "Intervention", y = "Number of hiccups") +
  coord_cartesian(ylim = c(0, 20)) +
  scale_y_continuous(breaks = seq(0, 20, 5)) +
  theme_minimal()
```

## Mixed designs

```{r}
#| label: grammar_data


grammar_tib <- here::here("data/social_media.csv") |> 
  read_csv() |> 
  mutate(
    time = as_factor(time),
    media_use = as_factor(media_use)
    )

# grammar_tib <- discovr::social_media

# View the ordered tibble

grammar_tib |> 
  arrange(id)
```

```{r}
#| label: grammar_plot

ggplot(grammar_tib, aes(x = time, y = grammar, colour = media_use)) +
  stat_summary(fun.data = "mean_cl_normal", geom = "pointrange", position = position_dodge(width = 0.05)) +
  labs(x = "Time", y = "Grammar score (%)", colour = "Media use") +
  coord_cartesian(ylim = c(0, 75)) +
  scale_y_continuous(breaks = seq(0, 75, 5)) +
  theme_minimal()

# now create the plot but with a line connecting the dots

ggplot(grammar_tib, aes(x = time, y = grammar, colour = media_use)) +
  stat_summary(fun = "mean", geom = "line", position = position_dodge(width = 0.05), aes(group = media_use)) +
  stat_summary(fun.data = "mean_cl_normal", geom = "pointrange", position = position_dodge(width = 0.05)) +
  labs(x = "Time", y = "Grammar score (%)", colour = "Media use") +
  coord_cartesian(ylim = c(0, 75)) +
  scale_y_continuous(breaks = seq(0, 75, 5)) +
  theme_minimal()
```

## Scatterplots

```{r}
#| label: scatter_data

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


# exam_tib  <- discovr::exam_anxiety

ong_tib <- here::here("data/ong_2011_tidy.csv") |> 
  read_csv() |> 
  mutate(
    id = as_factor(id),
    sex = as_factor(sex),
    rating_type = as_factor(rating_type)
  )

# ong_tib <- discovr::ong_tidy
```

### Bivatriate scatterplot

```{r}
#| label: biv_scat

#We'll build up the plot piece by piece for pedagogic reasons, but in realioty we'd do this in a sinfgle block of code

# A rudimentary scatterplot
ggplot(exam_tib, aes(x = anxiety, y = exam_grade)) +
  geom_point() 

# Add some more descriptive labels 
ggplot(exam_tib, aes(x = anxiety, y = exam_grade)) +
  geom_point() +
  labs(x = "Exam anxiety", y = "Exam performance (%)")

# scale the axes
ggplot(exam_tib, aes(x = anxiety, y = exam_grade)) +
  geom_point() +
  labs(x = "Exam anxiety", y = "Exam performance (%)") +
  coord_cartesian(ylim = c(0, 130)) +
  scale_y_continuous(breaks = seq(0, 130, 10)) +
  scale_x_continuous(breaks = seq(0, 100, 10))

# add a theme
ggplot(exam_tib, aes(x = anxiety, y = exam_grade)) +
  geom_point() +
  labs(x = "Exam anxiety", y = "Exam performance (%)") +
  coord_cartesian(ylim = c(0, 130)) +
  scale_y_continuous(breaks = seq(0, 130, 10)) +
  scale_x_continuous(breaks = seq(0, 100, 10)) +
  theme_minimal()

# add a summary line
ggplot(exam_tib, aes(x = anxiety, y = exam_grade)) +
  geom_smooth() +
  geom_point() +
  labs(x = "Exam anxiety", y = "Exam performance (%)") +
  coord_cartesian(ylim = c(0, 130)) +
  scale_y_continuous(breaks = seq(0, 130, 10)) +
  scale_x_continuous(breaks = seq(0, 100, 10)) +
  theme_minimal()


# make the summary line stroight and change its properties
ggplot(exam_tib, aes(x = anxiety, y = exam_grade)) +
  geom_smooth(method = "lm", colour = "#ef4836", fill = "#ef4836") +
  geom_point() +
  labs(x = "Exam anxiety", y = "Exam performance (%)") +
  coord_cartesian(ylim = c(0, 130)) +
  scale_y_continuous(breaks = seq(0, 130, 10)) +
  scale_x_continuous(breaks = seq(0, 100, 10)) +
  theme_minimal()

```

### Grouped scatterplot

```{r}
#| label: gp_scatter

# a basic plot
ggplot(exam_tib, aes(x = anxiety, y = exam_grade, colour = sex)) + 
   geom_point() +
   geom_smooth(method = "lm")

#colour the confidence intervals by sex
ggplot(exam_tib, aes(x = anxiety, y = exam_grade, colour = sex)) + 
   geom_point() +
   geom_smooth(method = "lm", aes(fill = sex), alpha = 0.1)

# add some labels
ggplot(exam_tib, aes(x = anxiety, y = exam_grade, colour = sex)) + 
   geom_point() +
   geom_smooth(method = "lm", aes(fill = sex), alpha = 0.1) +
  labs(x = "Exam anxiety", y = "Exam performance (%)", colour = "Sex assigned at birth", fill = "Sex assigned at birth")

# scale the axes

ggplot(exam_tib, aes(x = anxiety, y = exam_grade, colour = sex)) + 
   geom_point() +
   geom_smooth(method = "lm", aes(fill = sex), alpha = 0.1) +
  labs(x = "Exam anxiety", y = "Exam performance (%)", colour = "Sex assigned at birth", fill = "Sex assigned at birth") +
  coord_cartesian(ylim = c(0, 130)) +
  scale_y_continuous(breaks = seq(0, 130, 10)) +
  scale_x_continuous(breaks = seq(0, 100, 10)) +
  theme_minimal()

```

### Plotting multiple variables

```{r}
#| label: ggpairs_demo

GGally::ggpairs(ong_tib, columns = c("rating", "extraversion", "narcissism"))

# add mapping to colour by rating type
GGally::ggpairs(ong_tib,
                columns = c("rating", "extraversion", "narcissism"),
                mapping = aes(colour = rating_type, alpha = 0.1)) + 
  theme_minimal()
```
