lec04-dplyr.Rmd

---
title: Data wrangling and visualization in the tidyverse
author: Joel Östblom
---

## Lesson preamble

> ### Learning Objectives
>
> - Understand the split-apply-combine concept for data analysis.
> - Use `summarize`, `group_by`, and `tally` to split a data frame into groups
> of observations, apply a summary statistics for each group, and then combine
> the results.
> - Produce scatter plots, line plots, and histograms using ggplot.
> - Set universal plot settings.
> 
> ### Lesson outline
>
> - Split-apply-combine techniques in **`dplyr`** (25 min)
> - Using `tally` to summarize categorical data (15 min)
> - Plotting with **`ggplot2`** (20 min)
> - Building plots iteratively (25 min)

-----

## Setting up

Start by loading the required packages. Both **`ggplot2`** and **`dplyr`** are
included in the **`tidyverse`** package collection.

```{r}
# Install if needed
# install.packages('tidyverse')
library(tidyverse)
```

Load the data we saved in the previous lesson.

```{r, eval=FALSE}
# Download if needed
# download.file("https://ndownloader.figshare.com/files/2292169", "data/portal_data.csv")
surveys <- read_csv('portal_data.csv')
```

```{r, echo=FALSE}
surveys <- read_csv('data/portal_data.csv')
```

```{r}
surveys
```


## Split-apply-combine techniques in dplyr

Many data analysis tasks can be approached using the *split-apply-combine*
paradigm: split the data into groups, apply some analysis to each group, and
then combine the results.

**`dplyr`** facilitates this workflow through the use of `group_by()`
to split data and `summarize()`, which collapses each group into a single-row
summary of that group. The arguments to `group_by()` are the column names that
contain the **categorical** variables for which you want to calculate the
summary statistics. Let's view the mean `weight` by sex.

```{r}
surveys %>%
    group_by(sex) %>%
    summarize(mean_weight = mean(weight))
```

The mean weights become `NA` since there are individual observations that are
`NA`. Let's remove those observations.

```{r}
surveys %>%
    filter(!is.na(weight)) %>%
    group_by(sex) %>%
    summarize(mean_weight = mean(weight))
```

There is one row here that is neither male nor female, these are observations
where the animal escaped before the sex could not be determined. Let's remove
those as well.

```{r}
surveys %>%
    filter(!is.na(weight) & !is.na(sex)) %>%
    group_by(sex) %>%
    summarize(mean_weight = mean(weight))
```

You can also group by multiple columns:

```{r}
surveys %>%
    filter(!is.na(weight) & !is.na(sex)) %>%
    group_by(genus, sex) %>%
    summarize(mean_weight = mean(weight))
```

Since we will use the same filtered and grouped data frame in multiple code
chunks below, we could assign this subset of the data to a new variable and use
this variable in the subsequent code chunks instead of typing out the functions
each time.

```{r}
filtered_surveys <- surveys %>%
    filter(!is.na(weight) & !is.na(sex)) %>%
    group_by(genus, sex)
```

If you want to display more data, you can use the `print()` function at the end
of your chain with the argument `n` specifying the number of rows to display.

```{r}
filtered_surveys %>%
    summarize(mean_weight = mean(weight)) %>%
    print(n = 15) # Will change the knitted output, not the notebook
```

Once the data are grouped, you can also summarize multiple variables at the same
time. For instance, we could add a column indicating the minimum weight for each
species for each sex:

```{r}
filtered_surveys %>%
    summarize(mean_weight = mean(weight),
              min_weight = min(weight))
```

#### Challenge

1. Use `group_by()` and `summarize()` to find the mean, min, and max hindfoot
length for each species.

2. What was the heaviest animal measured in each year? Return the columns `year`,
`genus`, `species`, and `weight`.

```{r, include=FALSE}
## Answer 1
surveys %>%
    filter(!is.na(hindfoot_length)) %>%
    group_by(species) %>%
    summarize(
        mean_hindfoot_length = mean(hindfoot_length),
        min_hindfoot_length = min(hindfoot_length),
        max_hindfoot_length = max(hindfoot_length)
    )
## Answer 2
surveys %>%
    filter(!is.na(weight)) %>%
    group_by(year) %>%
    filter(weight == max(weight)) %>% # This is going to compare to the max weight within each group
    select(year, genus, species, weight) %>%
    arrange(year)
```


### Using tally to summarize categorical data

When working with data, it is also common to want to know the number of
observations found for each factor or combination of factors. For this, **`dplyr`**
provides `tally()`. For example, if we want to group by taxa and find the
number of observations for each taxa, we would do:

```{r}
surveys %>%
    group_by(taxa) %>%
    tally()
```

We can also use `tally()` when grouping on multiple variables:

```{r}
surveys %>%
    group_by(taxa, sex) %>%
    tally()
```

Here, `tally()` is the action applied to the groups created by `group_by()` and
counts the total number of records for each category.

If there are many groups, `tally()` is not that useful on its own. For example,
when we want to view the five most abundant species among the observations:

```{r}
surveys %>%
    group_by(species) %>%
    tally()
```

Since there are 40 rows in this output, we would like to order the table to
display the most abundant species first. In `dplyr`, we say that we want to
`arrange()` the data.

```{r}
surveys %>%
    group_by(species) %>%
    tally() %>%
    arrange(n)
```

Still not that useful. Since we are interested in the most abundant species, we
want to display those with the highest count first, in other words, we want to
arrange the column `n` in descending order:

```{r}
surveys %>%
    group_by(species) %>%
    tally() %>%
    arrange(desc(n)) %>%
    head(5)
```

If we want to include more attributes about these species, we can include these
in the call to `group_by()`:

```{r}
surveys %>%
    group_by(species, taxa, genus) %>%
    tally() %>%
    arrange(desc(n)) %>%
    head(5)
```

Be careful not to include anything that would split the group into subgroups,
such as `sex`, `year` etc.

#### Challenge

1. How many individuals were caught in each `plot_type` surveyed?

2. You saw above how to count the number of individuals of each `sex` using a
combination of `group_by()` and `tally()`. How could you get the same result
using `group_by()` and `summarize()`? Hint: see `?n`.


```{r, include=FALSE}
## Answer 1
surveys %>%
    group_by(plot_type) %>%
    tally()

## Answer 2
surveys %>%
  group_by(sex) %>%
  summarize(n = n())
```


## Plotting with ggplot2

**`ggplot2`** is a plotting package that makes it simple to create complex plots
from data frames. The name **`ggplot2`** comes from its inspiration, the book "A
grammar of graphics", and the main goal is to allow coders to express
their desired outcome on a high level instead of telling the computer every
detail about what will happen. For example, you would say "color my data by
species" instead of "go through this data frame and plot any observations of
species1 in blue, any observations of species2 in red, etc". Thanks to this
functional way of interfaces with data, only minimal changes are required if the
underlying data change or to change the type of plot. This helps in thinking
about the data and creating publication quality plots with minimal amounts of
adjustments and tweaking.

ggplot graphics are built step by step by adding new elements, or layers. Adding layers in
this fashion allows for extensive flexibility and customization of plots. To
build a ggplot, we need to:

1. Use the `ggplot()` function and bind the plot to a specific data frame using the
      `data` argument

```{r}
ggplot(data = surveys)
```

Remember, if the arguments are provided in the right order then the names of the
arguments can be omitted.

```{r}
ggplot(surveys)
```

2. Define aesthetics (`aes`), by selecting the variables to be plotted and the
   variables to define the presentation such as plotting size, shape color, etc.

```{r}
ggplot(surveys, aes(x = weight, y = hindfoot_length))
```

3. Add `geoms` -- geometrical objects as a graphical representation of the data
in the plot (points, lines, bars). **`ggplot2`** offers many different geoms; we
will use a few common ones today, including:
      * `geom_point()` for scatter plots, dot plots, etc.
      * `geom_line()` for trend lines, time-series, etc.
      * `geom_histogram()` for histograms

To add a geom to the plot use `+` operator. Because we have two continuous
variables, let's use `geom_point()` first:

```{r}
# If this takes way too long on your machine, create a subset from a random
# sample of a suitable size and continue working with this instead of `survey`.
#survey_subset <- sample_n(surveys, size = 5000)

ggplot(surveys, aes(x = weight, y = hindfoot_length)) +
  geom_point()
```

The `+` in the **`ggplot2`** package is particularly useful because it allows you
to modify existing `ggplot` objects. This means you can easily set up plot
"templates" and conveniently explore different types of plots, so the above
plot can also be generated with code like this:

```{r, first-ggplot-with-plus}
# Assign plot to a variable
surveys_plot <- ggplot(surveys, aes(x = weight, y = hindfoot_length))

# Draw the plot
surveys_plot + geom_point()
```

Notes:

- Anything you put in the `ggplot()` function can be seen by any geom layers
  that you add (i.e., these are universal plot settings). This includes the x and
  y axis you set up in `aes()`.
- You can also specify aesthetics for a given geom independently of the
  aesthetics defined globally in the `ggplot()` function.
- The `+` sign used to add layers must be placed at the end of each line containing
a layer. If, instead, the `+` sign is added in the line before the other layer,
**`ggplot2`** will not add the new layer and R will return an error message.


### Building plots iteratively

Building plots with ggplot is typically an iterative process. We start by
defining the dataset we'll use, lay the axes, and choose a geom:

```{r}
ggplot(surveys, aes(x = weight, y = hindfoot_length)) +
    geom_point()
```

Then, we start modifying this plot to extract more information from it. For
instance, we can add transparency (`alpha`) to reduce overplotting:


```{r}
ggplot(data = surveys, aes(x = weight, y = hindfoot_length)) +
    geom_point(alpha = 0.2)
```

Based on the hindfoot length and the weights, there appears to be 4-5 clusters
in this data. Potentially, one of the categorical variables we have in the data
could explain this pattern. Coloring the data points according to a
categorical variable is an easy way to find out if there seems to be
correlation. Let's try this with `plot_type`.

```{r}
ggplot(surveys, aes(x = weight, y = hindfoot_length, color = plot_type)) +
    geom_point(alpha = 0.2)
```

It seems like the type of plot the animal was captured on correlates well with
some of these clusters, but there are still many that are quite mixed. Let's try
to do better! This time, the information about the data can provide some clues
to which variable to look at. The plot above suggests that there might be 4-5
clusters, so a variable with 4-5 values is a good guess for what could explain
the observed pattern in the scatter plot.

```{r}
surveys %>%
    summarize_all(n_distinct)
```

Remember that there are still `NA` values here, that's why there appears to be
three sexes although there is only male and female. There are four taxa so that
could be a good candidate, let's see which those are.

```{r}
surveys %>%
    distinct(taxa)
```

It seems reasonable that these taxa contain animals different enough to have
diverse weights and length of their feet. Lets use this categorical variable to
color the scatter plot.

```{r}
ggplot(surveys, aes(x = weight, y = hindfoot_length, color = taxa)) +
    geom_point(alpha = 0.2)
```

Only rodents? That was unexpected... Let's check what's going on.

```{r}
surveys %>%
    group_by(taxa) %>%
    tally()
```

There is definitely mostly rodents in our data set...

```{r}
surveys %>%
    filter(!is.na(hindfoot_length)) %>% # control by removing `!`
    group_by(taxa) %>%
    tally()
```

...and it turns out that only rodents, have had their hindfeet measured!

Let's remove all animals that did not have their hindfeet measured, including
those rodents that did not. Animals without their weight measured will also be
removed.

```{r}
surveys_hf_wt <- surveys %>%
    filter(!is.na(hindfoot_length) & !is.na(weight))

surveys_hf_wt %>%
    summarize_all(n_distinct)
```

Maybe the genus can explain what we are seeing.

```{r}
ggplot(surveys_hf_wt, aes(x = weight, y = hindfoot_length, color = genus)) +
    geom_point(alpha = 0.2)
```

Now this looks good! There is a clear separation between different genus, but
also significant spread within genus, for example in the weight of the green
Neotoma observations. There are also two clearly separate clusters that are both
colored in olive green (Dipodomys). Maybe separating the observations into
different species would be better?

```{r}
ggplot(surveys_hf_wt, aes(x = weight, y = hindfoot_length, color = species)) +
    geom_point(alpha = 0.2)
```

Great! Together with the genus plot, this definitely seem to explain most of the
variance we see in the hindfoot length and weight measurements. It is still a
bit messy as it appears like we have around 5 clusters, but there are 21 species
in the legend.

```{r}
surveys %>%
    filter(!is.na(hindfoot_length) & !is.na(weight)) %>%
    group_by(species) %>%
    tally() %>%
    arrange(desc(n))
```

There is a big drop from 838 to 159, let's include only those with more than 800
observations.

```{r}
surveys_abun_species <- surveys %>%
    filter(!is.na(hindfoot_length) & !is.na(weight)) %>%
    group_by(species) %>%
    mutate(n = n()) %>% # add count value to each row
    filter(n > 800) %>%
    select(-n)

surveys_abun_species
```

Still has almost 25k observations, so only 10k was removed.

```{r}
ggplot(surveys_abun_species, aes(x = weight, y = hindfoot_length, color = species)) +
    geom_point(alpha = 0.2)
```


#### Challenge

Create a scatter plot of `hindfoot_length` over `species` with the `weight` showing in different colors.
Is there any problem with this plot? *Hint: think about how many observations there are*

```{r, include=FALSE}
ggplot(surveys_abun_species, aes(x = weight, y = species, color = hindfoot_length)) +
    geom_point(size = 0.1, position = 'jitter')
```


*Parts of this lesson material were taken and modified from [Data
Carpentry](https://datacarpentry.org) under their CC-BY copyright license. See
their [lesson page](https://datacarpentry.org/R-ecology-lesson/03-dplyr.html)
for the original source.*