R-data-manipulation.Rmd

---
title: "Data manipulation and analysis using tidyverse R"
author: "Ben Bond-Lamberty"
date: "`r Sys.Date()`"
output: slidy_presentation
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

library(dplyr)

library(ggplot2)
theme_set(theme_bw())
```

# Introduction {.bigger}

## The plan

A short workshop covering reproducibility and code/data management; data filtering, reshaping, and joining; and summarizing.

* Reproducible research and data management (~15 minutes)
* Filtering and cleaning data (~30 minutes; the `gapminder` dataset)
* Reshaping and joining data (~30 minutes)
* Summarizing and manipulating data (~60 minutes; the `babynames` dataset)

Slides: https://rpubs.com/bpbond/874167

## Requirements

**This workshop assumes a basic to intermediate knowledge of R.**

If you want to do the hands-on exercises (encouraged!), make sure up-to-date
versions of the following packages are installed:

- `dplyr` - fast, flexible tool for working with data frames
- `tidyr` - reshaping and cleaning data
- `ggplot2` - popular package for visualizing data
- `gapminder` - life expectancy, GDP per capita, and population for 142 countries
- `babynames` - names provided to the SSA 1880-2017

## A note about 'tidyverse'

What we're working with today--the "tidyverse"--constitute a [particular and popular dialect of R](http://tidyverse.org),
but the principles we'll go over are broadly applicable. 

Like anything else, `dplyr` has advantages and disadvantages. In particular, **I do not recommend using it if you don't already have experience with base R**.

I will point out base R equivalents as we go.

More information: https://github.com/matloff/TidyverseSkeptic

**The `data.table` package** is also worth checking out for its speed,
stability, and zero dependencies.


# Reproducibility and data management {.bigger}

## Reproducible research...

We are in the era of collaborative 'big data', but even if you work by yourself
with 'little data' you have to have some skills to deal with those data.

<div class='left' style='float:left;width:48%'>
**Most fundamentally, your results have to be reproducible.**

>Your most important collaborator is your future self. It’s important to make a workflow that you can use time and time again, and even pass on to others in such a way that you don’t have to be there to walk them through it. [Source](http://berkeleysciencereview.com/reproducible-collaborative-data-science/)
</div>
<div class='right' style='float:right;width:48%'>
<img src="data-manipulation-images/future_self.png" width="250" />

 <font size="2"> https://xkcd.com/1421/</font> 

</div>

## ...is the future

<div class='left' style='float:left;width:48%'>
**Prepare yourself for the future**. Funders, 
journals, governments, colleagues are all pushing for more reproducibility and 
openness: open-access journals, open code,
and [data deposition](https://www.agu.org/-/media/Files/Publications/JGR-Biogeo-Data-Sharing-Guidance.pdf) 
at the peer review stage and beyond.

It's a slow but steady ratchet.

Reproducibility generally means *scripts* tied to *open source software* with 
effective *data management* and *archiving*.
</div>
<div class='right' style='float:right;width:48%'>
<img src="data-manipulation-images/Ratchet_example.gif" width="400" />

 <font size="2"> By ZabMilenko at English Wikipedia, CC BY 3.0, https://commons.wikimedia.org/w/index.php?curid=49717967</font> 

</div>

## You can't reproduce

...what you've lost. What if you need access to a file as it existed 1, 10, or 100, or 1000 days ago?

<div class='left' style='float:left;width:48%'>
- Incremental backups (minimum)
- Don't depend on yourself! Has to be _automatic_.
- Version control. A *repository* holds files and tracks changes: what, by whom, why
</div>
<div class='right' style='float:right;width:48%'>
<img src="data-manipulation-images/tardis.jpg" width="400" />
</div>

## Version control

**Git** (and website [GitHub](https://github.com)) are the most popular version
control tools for use with R, and many other languages. They offer:

<div class='left' style='float:left;width:48%'>
- Version control
- Sharing work with collaborators in a *repository*
- Issue tracking
- Public or private code
</div>
<div class='right' style='float:right;width:48%'>
<img src="data-manipulation-images/git_2x.png" width="300" />

 <font size="2"> https://xkcd.com/1597/ </font> 

</div>

## Scripts for data analysis

Version control and scripts address two of the biggest problems with 
managing code and data:

* Tracking *changes over time*
* Understanding and reproducing the *analytical steps*.

**Ideally, _every_ step in your analysis is programmatic**. This means it is performed by a script that can be read and understood in the future.

Manual steps are hard to reproduce, laborious, and error-prone.

## Reproducibility is a process

*Don't let the perfect be the enemy of the good.* 

Upgrade and improve your workflow and skills over time:

>Organizing analyses so that they are reproducible is not easy. It requires diligence and a considerable investment of time: to learn new computational tools, and to organize and document analyses as you go.

>But partially reproducible is better than not at all reproducible. Just try to make your next paper or project better organized than the last.

A great and practical guide: http://kbroman.org/steps2rr/

## Reproducible research example

A typical project/paper directory for me, slightly idealized:
```
1-download.R
2-prep_data.R
3-analyze_data.R
4-manuscript_report.Rmd
logs/
output/
rawdata/
```

This directory contains *scripts* that are backed up both *locally* and *remotely*.
It is under *version control*, so it's easy to track changes over time.

There's also [targets](https://docs.ropensci.org/targets/) and [make](https://en.wikipedia.org/wiki/Make_(software)), i.e. fully 
automated analysis workflows...but that's a topic for another day.


# Filtering and cleaning data {.bigger}

## gapminder

In honor of the late [Hans Rosling](https://en.wikipedia.org/wiki/Hans_Rosling), we'll use the `gapminder` dataset.

```{r, message=FALSE}
library(gapminder)
gapminder
```

Read the help page about `gapminder`.

## Pipes and pipelines

The `magrittr` package (used by both `dplyr` and `tidyr`) provides the `%>%` operator, which allows us to _pipe_ an object forward into a function or call expression.

Note that `x %>% f` is _usually_ equivalent to `f(x)`.

```{r, eval=FALSE}
print(gapminder)

library(dplyr)
gapminder %>% print()
gagminder %>% head()
gapminder %>% head(n=20)
gapminder %>%
  print() %>%
  summary()    # what is non-piped equivalent?
summary(print(gapminder))
```

RStudio has a [keyboard shortcut](https://support.rstudio.com/hc/en-us/articles/200711853-Keyboard-Shortcuts-in-the-RStudio-IDE) for this.

R now (as of v4.1) has its own built-in pipe operator: `|>`.

## dplyr

The `dplyr` package uses _verbs_ (functions) to operate on _tibbles_ (data frames).
When using pipes, this looks like:

```{r, eval=FALSE}
some_data_frame %>%

  do_something() %>%

  do_something_else() %>%

  getting_dizzy_from_so_much_doing()
```

As noted above, **pipelines have tradeoffs** like any other tool.
But they do provide easy-to-read algorithmic steps in many cases.

Let's go over some of those possible `do_something` steps.

## Filter

Very commonly used, as it lets you _filter_ a dataset by one or more conditions.

```{r, eval=FALSE}
gapminder %>% filter(country == "Egypt")
gapminder %>% filter(country == "Egypt", year > 2000) # AND
gapminder %>% filter(country == "Egypt" | year > 2000) # OR
gapminder %>% 
    filter(country %in% c("Egypt", "New Zealand", "Chad"))
```

This is equivalent to base R's `subset` function.

```{r, eval=FALSE}
subset(gapminder, country == "Egypt")
```

## Filter

```{r}
library(ggplot2)
gapminder %>%
  filter(year == 1997) %>%
  ggplot(aes(gdpPercap, lifeExp, size = pop, color = continent)) +
  geom_point() +
  scale_x_log10()
```

## Select

Also extremely useful. Note the different notations for selecting columns:

```{r, eval=FALSE}
select(gapminder, pop, year)
gapminder %>% select(pop, year)
# Base R equivalent: gapminder[c("pop", "year")]
gapminder %>% select(-lifeExp, -gdpPercap)
gapminder %>% select(-1)
```

There are lots of other cool ways to select columns--see `?select`.

```{r, eval=FALSE}
gapminder %>% select(starts_with("c"))
gapminder %>% select(where(is.numeric))
```

In base R:

```{r, eval=FALSE}
gapminder[c("pop", "year"),]
```

## <span style="color: red;">Exercise: filtering and selecting</span>

Let's focus on a single country's data for a bit. Write a pipeline that picks 
out Egypt data only, removes the continent and country columns, and assigns 
the result to a variable `Egypt`. How many rows does the resulting dataset have?

## <span style="color: red;">Exercise: filtering and selecting</span>

Let's focus on a single country's data for a bit. Write a pipeline that picks 
out Egypt data only, removes the continent and country columns, and assigns 
the result to a variable `Egypt`. How many rows does the resulting dataset have?

```{r}
gapminder %>%
  filter(country == "Egypt") %>%
  select(-continent, -country) ->
  Egypt
```

```{r, echo=FALSE}
Egypt
```

## Uniting and separating

These functions can be very useful.

```{r}
library(tidyr)
gapminder %>% unite(coco, country, continent)
gapminder %>%
  unite(coco, country, continent) %>%
  separate(coco,
           into = c("country", "continent"),
           sep = "_")
```

## Renaming and mutating

The `mutate` function in particular is used a _lot_ in `dplyr` pipelines.

```{r}
Egypt %>%
    rename(population = pop) %>% 
    mutate(logpop = log(population)) 

# base R:
# gapminder$logpop <- log(pop)
# colnames(gapminder)[colnames(gapminder) == "pop"] <- "population"
```

## Mutating

Several important notes:

1. You can have multiple assignments within a single `mutate`; this is more efficient than separate `mutate` calls
2. New variables overwrite existing variables of the same name

```{r, eval=FALSE}
# This will not work correctly
gapminder %>% 
    mutate(pop = mean(pop),
           pop_sd = sd(pop))
```

4. Variables can be removed by setting their value to NULL
5. There's also `transmute`, which adds new variables and drops existing ones

## Mutating

Several important notes:

6. _Row-wise operations_ require special handling:

```{r}
df <- tibble(x = 1:2, y = 3:4, z = 5:6)
df %>% mutate(m = mean(c(x, y, z)))
df %>% rowwise() %>% mutate(m = mean(c(x, y, z)))
```

**Note** that for a large dataset, base R's `rowMeans` function will be _much_ 
faster than this, however. Use the best tool for the job!

## Data cleaning

Two handy `tidyr` functions allow you to fill data based on adjacent values, 
and fill in `NA` values.

```{r}
people <- tibble(who = c("Alice", "Bob", "Carol"), age = c(25, NA, 45))
people %>% fill(age, .direction = "down")
people %>% replace_na(list(age = 99))
```


# Reshaping and joining data {.bigger}

## Reshaping ("pivoting") data

Let's start with the simple `people` data frame, but adding another column:

```{r}
people$height <- c(160, 170, 180)
people
```

Note that there is one column of _metadata_ (name) and two of _data_ 
(age and height).

## Make data longer

Let's say we want to put `people` into long (or _tidy_) format--where every row
is a different observation. In other words, we need to collapse multiple
columns of data into a single one.

For this we use `tidyr::pivot_longer`, which at a minimum needs to know: 
what's the data source, and what column(s) do we want to pivot?

```{r}
people %>% pivot_longer(c("age", "height"))
```

(Note base R has a `reshape` function, but I find it difficult to use.)

## Make data longer

We could also say "pivot everything except for the 'who' column", or provide
custom column names in the resulting data frame:

```{r}
people %>% pivot_longer(-who)
people %>% pivot_longer(-who, names_to = "variable", values_to = "datum")
```

## Make data wider

To transform data into a wider form, `pivot_wider` needs to know at least
two things: which column holds the new column names, and which the data.

```{r}
people %>% 
    pivot_longer(-who) %>%
    pivot_wider(names_from = "name")
```

Here `tidyr` has guessed that if the names of the new columns are coming from
"name", then the values are probably coming from the numeric column "value".

And so we get our original `people` data frame back.

What happens if we take column names from the "who" column?

## Make data wider

Pivoting data wider may result in `NA` values if not every
row-column combination is present in the dataset:

```{r}
people %>% 
    pivot_wider(names_from = "who", values_from = "age")
people %>% 
    pivot_wider(names_from = "who", values_from = "age", values_fill = -1)
```

## Reshaping gapminder data

Time for something a bit more complex. Back to our gapminder-derived `Egypt`
dataset!

```{r}
Egypt %>% 
    pivot_longer(c(lifeExp, pop, gdpPercap))
```

## Reshaping gapminder data

```{r}
Egypt %>%
  pivot_longer(-year) %>%
  ggplot(aes(year, value)) + geom_line() +
   facet_wrap(~name, scales = "free")
```

Experiment. Why do these do what they do?

```{r, eval=FALSE}
Egypt %>% pivot_longer(lifeExp)
Egypt %>% pivot_longer(-lifeExp)
```

## <span style="color: red;">Exercise: reshaping</span>

Our `Egypt` dataset has one row per year, and three data columns. Write code
to transform it to have one row per datum, and each year in a
different column. Do this for post-2000 data only.

Hint: start by filtering and then pivoting it into _long_ form.

## <span style="color: red;">Exercise: reshaping</span>

Our `Egypt` dataset has one row per year, and three data columns. Write code
to transform it to have one row per datum, and each year in a
different column. Do this for post-2000 data only.

Hint: start by filtering and then pivoting it into _long_ form.

```{r}
Egypt %>% 
    filter(year > 2000) %>% 
    pivot_longer(-year) %>% 
    pivot_wider(names_from = "year")
```

## Joins

Frequently we have more than one data source:

```{r}
people
hometowns <- tibble(who = c("Alice", "Bob", "Dave"),
                    where = c("Abilene", "Billsville", "Darien"))
hometowns
```

Here, we'd like to `merge` (the base R function) or _join_ these together.

## Left joins

Probably the most commonly used join is `left_join`.
This works by including all rows in the _first_ data frame, and any
matches it finds from the _second_:

```{r}
people %>% left_join(hometowns)
```

In other words, it's a _lookup_.

Note the message. **It's good practice to specify the common columns** 
that joins should use:

```{r, eval=FALSE}
people %>% left_join(hometowns, by = "who")
```

## Left joins

The order matters. Why do these produce different results?

```{r}
people %>% left_join(hometowns, by = "who")
hometowns %>% left_join(people, by = "who")
```

## Other joins

There are many different join operations available.

```{r}
people %>% right_join(hometowns, by = "who")
people %>% inner_join(hometowns, by = "who")
people %>% full_join(hometowns, by = "who")
# anti_join returns all rows from x without a match in y;
# semi_join returns all rows from x with a match
people %>% anti_join(hometowns, by = "who")
```


# Summarizing data {.bigger}

## Summarizing and manipulating data

Thinking back to the typical data pipeline, we often want to summarize data by groups as an intermediate or final step. For example, for each subgroup we might want to:

* Compute mean, max, min, etc. (`n` → 1)
* Compute rolling mean and other window functions (`n` → `n`)
* Fit models and extract their parameters, goodness of fit, etc.

Specific examples:

* `gapminder`: what's the year of maximum GDP for each country?
* `babynames`: what's the most common name over time?

## Split-apply-combine

These are generally known as *split-apply-combine* problems.

<img src="data-manipulation-images/split_apply_combine.png" width="600" />

 <font size="2"> https://github.com/ramnathv/rblocks/issues/8</font> 

## Why we're focusing on dplyr

<div class='left' style='float:left;width:48%'>
**There are many ways to tackle split-apply-combine in R.**

The `dplyr` package specializes in data frames, but also allows you 
to work with remote, out-of-memory databases, using exactly the same tools, 
because it abstracts away *how* your data is stored.

`dplyr` is very fast for most, though not all, operations on _data frames_ (tabular data). But again, like any tool, it has some disadvantages too.
</div>
<div class='right' style='float:right;width:48%'>
<img src="data-manipulation-images/dplyr_logo.png" width="300" />

 <font size="2"> https://github.com/tidyverse/dplyr </font> 

</div>

## Grouping

`dplyr` verbs (functions) become particularly powerful when used in conjunction 
with *groups* we define in the dataset. Grouping doesn't change the data, 
but instead _groups_ it in preparation for the next operation we perform.

```{r}
gapminder %>%
  group_by(country)
```

In base R:

```{r, eval=FALSE}
split(gapminder, ~ pop + year) # but this creates a new object
```

## Summarising

Most frequently, that "next operation" is a 

* `summarise` (`n` → 1), or 
* `mutate` (`n` → `n`):

```{r}
gapminder %>%
  group_by(country) %>%
  summarise(maxpop = max(pop))

```

In base R:

```{r, eval=FALSE}
aggregate(. ~ country, data = gapminder, max)
```

## Summarising and columns

```{r}
gapminder %>%
  group_by(continent) %>%
  summarise(meanLifeExp = mean(lifeExp),
            maxpop = max(pop))
```

Notice that the only columns left, after summarising, are:

* The grouping variable(s), and 
* The summarised ones.

## Summarising and groups

By default, every summary operation removes the _last_ grouping variable,
while leaving the others intact:

```{r}
gapminder %>%
  group_by(continent, country) %>%
  summarise(max_LifeExp = max(lifeExp))
```

This makes it easy to do certain kinds of multi-step operations; for example,
if we want to know what the _mean maximum country life expectancy_ is by continent:

```{r}
gapminder %>%
  group_by(continent, country) %>%
  summarise(max_LifeExp = max(lifeExp), .groups = "drop_last") %>% 
  summarise(mean(max_LifeExp))
```

You can also re-group as needed--See the `.groups` argument--or `ungroup`.

## Operating on multiple columns and/or functions

We can apply a function to multiple columns...

```{r}
gapminder %>%
  group_by(country) %>%
  summarise(across(lifeExp:gdpPercap, max))
```

...or multiple functions to one or more columns:

```{r}
gapminder %>%
  group_by(country) %>%
  summarise(across(pop, c(min, max)))
```

## More complex summaries

We now have the tools to build up a long pipeline to, e.g., 
compute the min and max for all numeric variables and produce
a table with continent as columns headers,
and variable (gdpPercap, lifeExp, pop) and statistic as rows.

```{r,}
# We define this function to end up with nice column names
min_max <- list(
  min = ~min(.x, na.rm = TRUE), 
  max = ~max(.x, na.rm = TRUE)
)

gapminder %>% 
    filter(continent %in% c("Africa", "Americas", "Asia")) %>% 
    group_by(continent) %>% 
    # compute min and max for the numeric columns
    summarise(across(where(is.numeric), min_max)) %>% 
    pivot_longer(-continent) %>% 
    separate(name, into = c("variable", "stat")) %>%
    pivot_wider(names_from = "continent", values_from = "value")
```

## Introducing `babynames`

Explore `babynames` a bit. How many rows, columns does it have? How many unique names?

```{r}
library(babynames)
babynames
```

## Summarizing babynames

What does this calculate?

```{r, message=FALSE}
babynames %>%
    group_by(year, sex) %>%
    slice_max(prop)
```

The various `slice_` functions are quite handy.

## Summarizing babynames

```{r, echo=FALSE}
library(ggrepel)
babynames %>%
    group_by(year, sex) %>%
    slice_max(prop) %>% 
    ungroup() %>% 
    # the one tricky bit follows: for each sex x year group, we want the 
    # first name to be plotted, but not the following ones
    arrange(sex, year) %>% 
    mutate(change = c(TRUE, name[-1] != head(name, -1)),
           name2 = if_else(change, name, "")) %>% 
    # and now plot
    ggplot(aes(year, prop, color = name, label = name2)) + 
    geom_point() + 
    geom_text_repel(nudge_y = 0.01) + 
    facet_grid(sex~.) + 
    guides(colour = "none")
```

https://en.wikipedia.org/wiki/Linda_(1946_song)

## <span style="color: red;">Exercise: the `babynames` dataset</span>

Load the dataset using `library(babynames)`.

Read its help page. Look at its structure (rows, columns, summary).

Use `dplyr` to calculate the total number of names in the SSA database for each 
year and sex, with each sex's counts in a separate column. Hint: `n()`.

Make a graph or table showing how popular YOUR name has been over time (either 
its proportion, or rank).

## <span style="color: red;">Exercise: the `babynames` dataset</span>

Load the dataset using `library(babynames)`.

Read its help page. Look at its structure (rows, columns, summary).

Use `dplyr` to calculate the total number of names in the SSA database for each 
year and sex, with each sex's counts in a separate column. Hint: `n()`.

Make a graph or table showing how popular YOUR name has been over time (either 
its proportion, or rank).

```{r, eval=FALSE}
babynames %>% 
    group_by(year, sex) %>% 
    summarise(n = n()) %>% 
    pivot_wider(names_from = "sex", values_from = "n")

babynames %>%
    filter(name == "Benjamin") %>%
    ggplot(aes(year, n, color = sex)) +
    geom_point() +
    ggtitle("Benjamin")
```

## Another example: ranking

```{r}
babynames %>%
    group_by(year, sex) %>%
    mutate(rank = row_number(desc(n))) %>%
    filter(name %in% c("Benjamin", "Rachel")) %>%
    ggplot(aes(year, rank, color = name, shape = sex)) +
    geom_point() + scale_y_log10() 
    ggtitle("Benjamin and Rachel")
```

## Model-fitting

The `diamonds` dataset is included with `ggplot2`.

```{r, message=FALSE}
ggplot(diamonds, aes(carat, price, color = cut)) + 
    geom_point() + 
    geom_smooth(method = "lm") +
    scale_x_log10() + scale_y_log10()
```

## Model-fitting

This old code now generates a warning. Leaving it here, in the slides,
to make sure to mention `dplyr`'s penchant for changing (and breaking) behavior.

```{r}
library(broom) # to sweep up -- tidy -- model outputs
diamonds %>% 
    group_by(cut) %>% 
    do(mod = lm(log(price) ~ log(carat), data = .)) %>% 
    # need to change this to `reframe` (no other change needed, thankfully)
    summarise(cut = cut,
              tidy(mod))
```

## Important things we didn't talk about

- Getting your data _into_ R
- Working with non-text/tabular data
- Lists
- Handling dates and timestamps
- Plotting (though we had many examples)
- Efficient to high performance computing
- Reproducible documents
- Building your own package

Check out slide decks at https://rpubs.com/bpbond

But the tools we've covered here can do a lot!

# The end {.bigger}

## Thank you!

<div class='left' style='float:left;width:48%'>

**Feedback welcome.**

* <a href="mailto:bondlamberty@pnnl">bondlamberty@pnnl.gov</a>
* https://github.com/bpbond/R-workshops/issues

The _slides_ for this presentation are available here: https://rpubs.com/bpbond/874167

The _repository_ with the R code that generated the slides is here: https://github.com/bpbond/R-workshops/
</div>
<div class='right' style='float:right;width:48%'>
<img src="data-manipulation-images/talk-like-a-pirate.gif" width="450" />

 <font size="2"> Source?</font> 

</div>