3.1-Combining_data.qmd

---
title: "Constructing Complex Plots"
subtitle: "Combining Data Sources"
author:
  - Elizabeth King
  - Kevin Middleton
format:
  revealjs:
    theme: [default, custom.scss]
    standalone: true
    self-contained: true
    logo: QMLS_Logo.png
    slide-number: true
    show-slide-number: all
    link-external-newwindow: true
---

## This week

```{r}
#| label: setup
#| echo: false
#| warning: false
#| message: false

library(tidyverse)
library(cowplot)
ggplot2::theme_set(theme_cowplot())
library(paletteer)

library(Data4Ecologists)
library(abdData)
```

1. Combining data sources
2. Adding elements
3. Highlighting elements
4. Reusable themes


## "Complex plots"

- Combining `geom`s
- Data from different sources (`data.frame`s, regressions, etc.)
- Non-standard elements or axis labels (text with math symbols)
- Custom colors or gradients
- Highlighting different elements with color, labels or arrows


## Example

![](Images/CGCS.png)


## Example

![](Images/RNAi_DEseq.jpg)


## Example

![](Images/Climate.png)


## Zooplankton diversity {.scrollable}

Diversity of zooplankton (`Diversity`) prey in each of 5 replicate `Block`s of three `Treatment` levels.^[Svanbäck, R and DI Bolnick. 2007. Intraspecific competition drives increased resource use diversity within a natural population. *Proc R Soc London Ser B, Biol Sci* 274: 839-844.]

```{r}
ZP <- Zooplankton |> 
  rename(Treatment = treatment,
         Diversity = zooplankton,
         Block = block) |> 
  mutate(Treatment = case_when(
    Treatment == "control" ~ "Control",
    Treatment == "low" ~ "Low",
    Treatment == "high" ~ "High"
  ),
  Treatment = factor(Treatment, levels = c("Low", "Control", "High")),
  Block = factor(Block)) |> 
  as_tibble()

ZP
```


## Points by groups

```{r}
#| echo: true
#| output-location: slide
ggplot(ZP, aes(x = Treatment, y = Diversity)) +
  geom_point(size = 2,
             position = position_jitter(width = 0.05,
                                        seed = 6734747))
```


## Boxplot

- Median
- 1st and 3rd quartiles (25th and 75th percentiles)
- Whiskers calculated variously (here ±1.5 * IQR)
- "Outliers" outside the whiskers

```{r}
#| echo: true
#| output-location: slide
ggplot(ZP, aes(x = Treatment, y = Diversity)) +
  geom_boxplot()
```


## Boxplot with points

- Add points on top of the boxplot (`ggplot` adds `geom`s sequentially)
- Turn off outlier highlighting (`outlier.shape = NA`)

```{r}
#| echo: true
#| output-location: slide
ggplot(ZP, aes(x = Treatment, y = Diversity)) +
  geom_boxplot(outlier.shape = NA) +
  geom_point(size = 2,
             position = position_jitter(width = 0.05,
                                        seed = 6734747))
```


## Precomputing quantities

```{r}
#| echo: true
SEM <- function(x) {
  return(sd(x) / sqrt(length(x)))
}

ZP_means <- ZP |> 
  group_by(Treatment) |> 
  summarise(mean_Diversity = mean(Diversity),
            upper_bound = mean_Diversity + 2 * SEM(Diversity),
            lower_bound = mean_Diversity - 2 * SEM(Diversity),
            .groups = "drop")

ZP_means
```


## Plots from multiple data sources

- Start with an empty `ggplot()`
- Use `data` argument to `geom_`...

```{r}
#| echo: true
#| output-location: slide
ggplot() +
  geom_point(data = ZP, aes(x = Treatment, y = Diversity),
             size = 2,
             position = position_jitter(width = 0.05, seed = 6734747)) +
  geom_point(data = ZP_means,
             aes(x = Treatment, y = mean_Diversity),
             size = 4,
             color = "gray50") +
  geom_errorbar(data = ZP_means,
                aes(x = Treatment, ymin = lower_bound, ymax = upper_bound),
                width = 0.1, size = 1,
                color = "gray50")
```


## Reaction norms for paired data

Antibody levels in Red-winged blackbirds before and after experimental manipulation of testosterone^[Hasselquist D, JA Marsh, PW Sherman, JC Wingfield. 1999. Is avian immunocompetence suppressed by testosterone? *Behav Ecol Sociobiol* 45: 167-175.]

```{r}
BB <- Blackbirds |> 
  mutate(ID = factor(1:n())) |> 
  rename(Before = before,
         After = after) |> 
  select(ID, Before, After)

BB
```


## Pivoting

```{r}
#| echo: true
BB_long <- BB |> 
  pivot_longer(cols = -ID,
               names_to = "Timepoint",
               values_to = "Antibody_level")

BB_long
```


## Plotting reaction norms

- `group = ID` tells `geom_line()` to associate lines with IDs.

```{r}
#| echo: true
#| output-location: slide
ggplot(BB_long, aes(x = Timepoint, y = Antibody_level, group = ID)) +
  geom_point() +
  geom_line()
```


## Ordering factors

Also many options in the `forcats` package (`fct_` functions) when you already have a factor.

```{r}
#| echo: true
#| output-location: slide
BB_long <- BB_long |> 
  mutate(Timepoint = factor(Timepoint, levels = c("Before", "After")))

ggplot(BB_long, aes(x = Timepoint, y = Antibody_level, group = ID)) +
  geom_point() +
  geom_line()
```


## Adding a mean and error bars

- `stat_summary()` can add the output of summary functions.
- `mean_cl_boot` give the mean and bootstrapped 95% CI. (requires `Hmisc`)
- `aes(group = -1)` ungroups the data so that the statistics will calculate for each `Timepoint`

```{r}
#| echo: true
#| output-location: slide
ggplot(BB_long, aes(x = Timepoint, y = Antibody_level, group = ID)) +
  geom_point() +
  geom_line() +
  stat_summary(aes(group = -1),
               fun.data = "mean_cl_boot",
               color = "magenta",
               size = 1)
```


## Horizontal, Vertical, and `ab` lines

```{r}
#| echo: true
#| output-location: slide
ggplot() +
  geom_hline(yintercept = -5:5, color = "orangered") +
  geom_vline(xintercept = -5:5, color = "orangered") +
  geom_abline(slope = seq(1, 10, by = 0.5), intercept = 0)
```


## Horizontal, Vertical, and `ab` lines: refined

```{r}
#| echo: true
#| output-location: slide
ggplot() +
  geom_hline(yintercept = -5:5, color = "orangered") +
  geom_vline(xintercept = -5:5, color = "orangered") +
  geom_abline(slope = seq(1, 10, by = 0.5), intercept = 0) +
  coord_equal()
```


## Adding regression lines

1. `geom_smooth()`
    - Can handle a range of models from `lm()` to GLM and GAM
    - Splits by aesthetics automatically
2. Compute regression and use `predict()` or another helper function
    - You handle all the prediction
    - Necessary for more complex models (e.g., mixed/multilevel)


## Bird species richness {.scrollable}

Bird species richness in different habitat patches sampled in Jamaica.^[Kennedy CM, Marra PP, Fagan WF, Neel MC. (2010). Landscape matrix and species traits mediate responses of Neotropical resident birds to forest fragmentation in Jamaica. *Ecol Monogr* 80(4): 651-669.]

```{r}
BSR <- birds |> 
  rename(Patch_ID = patch,
         n_Species = S,
         Landscape_type = landscape,
         Area = area,
         log_Area = log.area.,
         Year = year) |> 
  filter(Year == 2005) |> 
  select(-Year) |> 
  mutate(Landscape_type = factor(Landscape_type)) |> 
  drop_na()
BSR
```

## `geom_smooth()`

Defaults to a loess smoother with confidence interval

```{r}
#| echo: true
#| output-location: slide
ggplot(BSR, aes(x = log_Area, y = n_Species, color = Landscape_type)) +
  geom_point() +
  geom_smooth()
```


## `geom_smooth(formula = y ~ x, method = "lm")`

"Standard" linear regression

```{r}
#| echo: true
#| output-location: slide
ggplot(BSR, aes(x = log_Area, y = n_Species, color = Landscape_type)) +
  geom_point() +
  geom_smooth(formula = y ~ x, method = "lm")
```


## Removing the confidence interval

`se = FALSE`

```{r}
#| echo: true
#| output-location: slide
p1 <- ggplot(BSR, aes(x = log_Area, y = n_Species, color = Landscape_type)) +
  geom_point() +
  geom_smooth(formula = y ~ x, method = "lm", se = FALSE)
p1
```


## Removing the grouping

One line through all the points: `aes(group = -1)`

```{r}
#| echo: true
#| output-location: slide
ggplot(BSR, aes(x = log_Area, y = n_Species, color = Landscape_type)) +
  geom_point() +
  geom_smooth(aes(group = -1), formula = y ~ x, method = "lm", se = FALSE)
```


## "Manual" prediction

```{r}
#| echo: true
#| output-location: slide
fm <- lm(n_Species ~ log_Area + Landscape_type, data = BSR)
summary(fm)
```

- Slope for `log_Area` common to all levels of `Landscape_type`
- Separate intercepts for each `Landscape_type`

## `predict()`

`predict()` has methods for many regression-like functions:

- `predict.lm()`
- `predict.glm()`
- etc.
- Type `predict.` and then Tab to see options

By default returns the predicted values for observed data.

```{r}
#| echo: true
#| output-location: slide
predict(fm)
```


## Grid construction with `expand_grid()`

`predict()` option `newdata` allows passing a set of values for which to predict.

- Generate predictors over a range of `log_Area` for each level of `Landscape_type`
- Variable names must match the right-hand side of the formula
- `tidyr` function `expand_grid()` (or `crossing()`)

```{r}
#| echo: true
#| output-location: slide
Preds <- expand_grid(
  log_Area = seq(-0.5, 2, length.out = 200),
  Landscape_type = levels(BSR$Landscape_type) |> factor()
)
Preds
```


## Generate predictions

- Pass `newdata` to `predict()`

```{r}
#| echo: true
Preds <- Preds |> 
  mutate(Predicted = predict(fm, newdata = Preds))
Preds
```


## Plot new predictions

```{r}
#| echo: true
#| output-location: slide
p2 <- ggplot() +
  geom_point(data = BSR,
             aes(x = log_Area, y = n_Species, color = Landscape_type)) +
  geom_line(data = Preds,
            aes(x = log_Area, y = Predicted, color = Landscape_type),
            size = 1)
p2
```


## Compare

```{r}
plot_grid(p1 +
            labs(title = "geom_smooth()") +
            theme(legend.position = "none"),
          p2 +
            labs(title = "predict()") +
            theme(legend.position = "none"),
          ncol = 2)

```


## Making confidence interval bands

95% prediction interval

- Interval expected to include 95% of all new observations
- Wider than the CI for the mean estimate

```{r}
#| echo: true
#| output-location: slide
Preds <- Preds |>
  bind_cols(predict(fm, newdata = Preds,
                    interval = "prediction", level = 0.95) |> 
              as_tibble())
Preds
```


## Plotting confidence interval bands

- `geom_line()` adds the mean prediction `fit`
- `geom_ribbon()` adds the shaded bands bound by `lwr` and `upr`

```{r}
#| echo: true
#| output-location: slide
ggplot() +
  geom_point(data = BSR,
             aes(x = log_Area, y = n_Species, color = Landscape_type)) +
  geom_line(data = Preds,
            aes(x = log_Area, y = fit, color = Landscape_type),
            size = 1) +
  geom_ribbon(data = Preds,
              aes(x = log_Area, ymin = lwr, ymax = upr, fill = Landscape_type))
```


## Reordering plot elements

- Hiding the legend for `geom_ribbon`
- Relabeling axes
- Changing the color palette

```{r}
#| echo: true
#| output-location: slide
ggplot() +
  geom_ribbon(data = Preds,
              aes(x = log_Area, ymin = lwr, ymax = upr, fill = Landscape_type),
              alpha = 0.1,
              show.legend = FALSE) +
  geom_line(data = Preds,
            aes(x = log_Area, y = fit, color = Landscape_type),
            size = 1) +
  geom_point(data = BSR,
             aes(x = log_Area, y = n_Species, color = Landscape_type)) +
  labs(y = "Number of Species", x = "Area (log Ha)") +
  scale_color_paletteer_d(`"rcartocolor::Bold"`, name = "Landscape Type") +
  theme(legend.position = "inside",
        legend.position.inside = c(0.8, 0.15))
```