Split-Apply-Combine

Outline

Today, we’re going to be covering two really important concepts

1. Piping
2. Split-Apply-Combine

Piping was introduced to R with the magrittr package (get it? ), and over the past 4 years has dramatically altered how programming in R works. Split-Apply-Combine is a data analysis technique that you may have already used, but is streamlined in the dplyr and purrr packages.

Setup

~2 minute setup

Create and save your R notebook for today. Clear the workspace in case there’s anything left over from last time. Load the important packages for today’s work:

install.packages("broom")

library(lsa2017)
library(tidyverse)
library(magrittr)
library(broom)

Also, download and unzip the file called joe_vowel_files.zip from Canvas.

---

Piping

This is the new functor we’ll be using today:

%>%

We’ll pronounce that “pipe.”

So far, every time we’ve used a function, like nrow(), we’ve passed it its arguments by putting them inside the parentheses:

nrow(iy_ah)

[1] 44818

What %>% does is take everything to its left, and injects it into the first argument of the function to its right.

iy_ah %>% nrow()

[1] 44818

We’ll pronounce this “piping iy_ah into nrow().”

This might not seem so earth shattering right away, but it really shines when you want to apply a function to the output of the previous function. Let’s say you wanted to check how many rows of “iy” data there is. Using normal functional application, you’d have to do this:

nrow(filter(iy_ah, plt_vclass == "iy"))

[1] 23699

Here, the first function you apply is the most deeply embedded, and the last function you apply is the outermost. Some syntacticians feel very comfortable with bracket notation like this, but it rapidly gets unreadable, and mis-matching closing parentheses is a very common kind of error. Here’s how these same exact operations look like with the pipes.

iy_ah %>%
  filter(plt_vclass == "iy") %>%
  nrow()

[1] 23699

This is much more readable, and the functions are visually ordered in the same order they’re applied.

The benefit of piping becomes most apparent if we revisit our data cleaning pipeline from our last meeting. Here’s iy_ah wide again.

iy_ah_wide

Last time, we made this data wide to long by creating intermediate data frames. However, we could do it all in a single shot by nesting functions like so:

spread(
  separate(
    gather(iy_ah_wide, 
           key = "vowel_formant", 
           value = "hz", 
           ah_F1:iy_F2), 
    "vowel_formant", 
    into = c("vowel", "formant"), 
    sep = "_"), 
  formant, 
  hz)

But good luck writing that twice, or sharing your code with anyone. Notice how the additional arguments for spread() are maximally far away from the actual call to the function at the top?

Piping makes it much better.

iy_ah_wide %>%
  gather(key = "vowel_formant", value = "hz", ah_F1:iy_F2) %>%
  separate("vowel_formant", into  = c("vowel", "formant"), sep = "_") %>%
  spread(formant, hz)

Another nice thing about these chains is that they can be built up incrementally. That is, just the first two lines will work:

iy_ah_wide %>%
  gather(key = "vowel_formant", value = "hz", ah_F1:iy_F2)

And so will just the fist three:

iy_ah_wide %>%
  gather(key = "vowel_formant", value = "hz", ah_F1:iy_F2) %>%
  separate("vowel_formant", into  = c("vowel", "formant"), sep = "_")

The same is not true of the embedded functional application, because it depends on matching parentheses between the beginning and end.

Idiom

I am going to recommend that you put a new line after every pipe, %>%, and indent the following line with two spaces when chaining together functions.

~5 minute activity

Re-create the fruit data frame:

fruit <- data.frame(person = c("Oakley", "Charlie"),
                 apples_bought = c(5, 3),
                 apples_ate = c(1, 2),
                 oranges_bought = c(5, 4),
                 oranges_ate = c(3, 3))

Using pipes, recreate the gather, separate and spread workflow from the last meeting to produce this table (you can double check the code from last time).

person	fruit	ate	bought
Charlie	apples	2	3
Charlie	oranges	3	4
Oakley	apples	1	5
Oakley	oranges	3	5

Split-Apply-Comine

When doing data analysis, you’re going to find yourself doing these following steps a lot:

1. Splitting the data up into subsets.
2. Applying some kind of function to those subsets.
3. Combining the results back together

Here is some fake data from an experiment of speakers’ fundamental frequency. Go ahead and run this code:

  pitch <- data_frame(speaker = rep(c("Oakley", "Charlie", "Azaria"), 
                                    times =  c(3,3,2)),
                       F0 = c(124, 109, 125, 103, 121, 123, 181, 206))

speaker	F0
Oakley	124
Oakley	109
Oakley	125
Charlie	103
Charlie	121
Charlie	123
Azaria	181
Azaria	206

One thing we might want to know is what the average F0 is for each speaker. Here’s how we do that, conceptually, with the Split-Apply-Combine workflow.

Split the data up

First, split the data up into subsets based on the “speaker” column:

speaker	F0
Azaria	181
Azaria	206

speaker	F0
Charlie	103
Charlie	121
Charlie	123

speaker	F0
Oakley	124
Oakley	109
Oakley	125

Apply some function to the data

In each subset, calculate the average F.

speaker	mean_F0
Azaria	193.5

speaker	mean_F0
Charlie	115.6667

speaker	mean_F0
Oakley	119.3333

Combine the result

Combine these results into a new table.

speaker	mean_F0
Azaria	193.5000
Charlie	115.6667
Oakley	119.3333

Split-Apply-Combine in dplyr

The dplyr package was written specifically to streamline the Split-Apply-Combine workflow. All dplyr “verbs” take a data frame as input, does something to them, and returns another data frame. Let’s start with the summarise() function.

summarise(data,
new_col1 = summary_function1,
new_col2 = summary_function2)

• data
  • The data frame you want to summarise
• summary functions
  • This isn’t a named argument. Here, you give summarise() a function that can take many values as input, and returns a single value as output.

For example, let’s calculate the mean F0 across all of the data in the pitch data frame.

pitch %>%
  summarise(mean_F0 = mean(F0))

This doesn’t look like much, because it’s just a one row data frame with one column called mean_F0 with one value, which is the mean of all of the F0 measurements from each speaker.

It’s also possible to do multiple summaries at once by just adding more summary functions to the summarise() verb.

pitch %>%
  summarise(mean_F0 = mean(F0),
            median_F0 = median(F0),
            sd_F0 = sd(F0),
            mad_F0 = mad(F0),
            range_F0 = max(F0)-min(F0),
            ndata = length(speaker))

The summarise() function is what we want to apply to our data, but for the Split-Apply-Combine workflow, we want to split the data according to the speaker, apply the summarise() function to each subset, then combine the results back together.

To split up the data, we first pass the pitch data frame to group_by().

group_by(data, groupA, groupB)

• data
  • The data frame you want to apply Split-Apply-Combine to.
• groupings
  • These are the names of columns in the data frame that you want to split the data up by. group_by() will create subsets of every unique combination of values in these columns.

When you apply group_by() to a data frame by itself, nothing really special happens:

pitch %>% 
  group_by(speaker)

But when you apply summarise() after group_by(), you get the mean F0 for each speaker:

pitch %>% 
  group_by(speaker) %>%
  summarise(mean_F0 = mean(F0))

Here’s what just happened.

1. group_by() to split the data frame up according to the unique values of speaker.
2. summarise() calculated the mean F0 within each subset.
3. They were then automatically combined back together.

~5 minute activity

Just to refresh your memories, the iy_ah data frame contains all data from the vowels /i:/ and /ɑ/ from the PNC. The column idstring is a unique ID for each speaker, and plt_vclass codes the vowel class.

Calculate the average F1 and F2 for /i:/ and /ɑ/ for each speaker.

Additional dplyr verbs

Here are a few more dplyr verbs that you can string together.

verb	description
filter()	We’ve already used this to subset dataframes.
summarise()	This takes a data frame, and outputs a new data frame based on the summary you asked for
mutate()	This takes a data frame, and adds additional columns based on the formula you give it
select()	This takes a data frame, and returns only the columns you ask for
arrange()	Reorders the rows of the data frame

The one function here that requires some additional instruction is mutate().

mutate(data, new_col = value)

• data
  • the data that you want to mutate
• new columns
  • Here, you name new columns, and provide some kind of value, usually using some other function.

The simplest and least useful use of mutate() adds a new column with a single value:

pitch %>%
  mutate(jawn = "foo")

However, a much more useful use of mutate() creates transformed columns. For example, sometimes when working with pitch, you might estimate pitch differences in terms of semitones. Every doubling of pitch is 12 semitones. It’s the case that if you calculate the difference between the log2 transform of any two numbers that are doubles of eachother, you get 1.

# 1 octave differences
log2(6) - log2(3)

[1] 1

log2(18) - log2(9)

[1] 1

# 2 octave differences
log2(12) - log2(3)

[1] 2

If you wanted to estimate, for each speaker, their pitch range in terms of semitones, step 1 would be to convert their F0 in Hz into log2(Hz).

pitch %>%
  mutate(log2_F0 = log2(F0))

Then, for every speaker, subtract their max log2_F0 from their min log2_F0.

pitch %>%
  mutate(log2_F0 = log2(F0)) %>%
  group_by(speaker)%>%
  summarise(octave_range = max(log2_F0) - min(log2_F0),
            semitone_range = octave_range * 12)

~5 minute activity

Try the following:

pitch %>%
  mutate(number = n())

What happened? What do you think n() does?

Now try this:

pitch %>%
  group_by(speaker)%>%
  mutate(number = n())

What’s different this time?

Split-Apply-Combine for other data types

Data frames are not the only kind of data that you might want to use the Split-Apply-Combine workflow on. For example, the joe_vowel_files directory (that you unzipped earlier) contains one .csv for each vowel class from when I was an interviewer once. This is an unconventional data organization scheme, but not unheard of. For example, you might have a separate file of vowel data for each speaker.

We can use a variant of the Split-Apply-Combine workflow to load and combine all of these files at once.

~2 minute activity

This part may require some assistance. The Sys.glob() function searchers for all files that match the pattern you specify. The command below says “List all files that end in csv that are in the ../data/joe_vowel_files/ directory”. The preceeding .. in the path means that it should start the search from one directory above where the R notebook is saved. If you have set up your course directory as recommended, this should work. If you don’t get output that looks like this, I’ll come around and try to figure out why.

Sys.glob("../data/joe_vowel_files/*csv")

 [1] "../data/joe_vowel_files/*hr.csv"
 [2] "../data/joe_vowel_files/Tuw.csv"
 [3] "../data/joe_vowel_files/ae.csv" 
 [4] "../data/joe_vowel_files/ahr.csv"
 [5] "../data/joe_vowel_files/aw.csv" 
 [6] "../data/joe_vowel_files/ay.csv" 
 [7] "../data/joe_vowel_files/ay0.csv"
 [8] "../data/joe_vowel_files/e.csv"  
 [9] "../data/joe_vowel_files/ey.csv" 
[10] "../data/joe_vowel_files/eyF.csv"
[11] "../data/joe_vowel_files/i.csv"  
[12] "../data/joe_vowel_files/iy.csv" 
[13] "../data/joe_vowel_files/iyF.csv"
[14] "../data/joe_vowel_files/iyr.csv"
[15] "../data/joe_vowel_files/o.csv"  
[16] "../data/joe_vowel_files/oh.csv" 
[17] "../data/joe_vowel_files/ow.csv" 
[18] "../data/joe_vowel_files/owF.csv"
[19] "../data/joe_vowel_files/owr.csv"
[20] "../data/joe_vowel_files/uh.csv" 
[21] "../data/joe_vowel_files/uw.csv" 

Assign this vector of file names to the variable vowel_files.

vowel_files <- Sys.glob("../data/joe_vowel_files/*csv")

If we apply read.csv() to the first vowel file name, it will read it in as a data frame.

read.csv(vowel_files[1])

But there are 21 vowel files, and what we don’t want to do is write out something like vowel_2 <- read_csv(vowel_files[2]) 21 times! Even copy-pasting will involve errors.

To Split-Apply-Combine read_csv() to the vector of file names vowel_files, we’ll use map().

map(list, function, …)

• list
  • the list of things that you want to split up
• function
  • the function that you want to apply to every item in the list
• ...
  • any extra arguments you need to pass to the function.

It’s easiest to see this at work:

all_vowels <- map(vowel_files, read.csv)

all_vowels[1]

[[1]]

NA

all_vowels[21]

[[1]]

NA

Here’s what has happened:

• map() applied read.csv() to each file name.
• read.csv() read each file in as a dataframe.
• map() returned a list of dataframes that read.csv() produced.

There’s one last step we need combine all of these dataframes together. bind_rows() will combine all of these dataframes into one big dataframe.

vowels_df <- bind_rows(all_vowels)

head(vowels_df)

Of course, this can be chained together with pipes:

vowel_files %>%
  map(read.csv)%>%
  bind_rows()

Unequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorUnequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vector

Recursive Data Frames

There is one last use of map() that is really useful, but might be a bit mind-bending for right now, but will hopefully start to make more sense as we work through the course, or when you re-read these notes. I wouldn’t necessarilly file this under “required core R competency”, but it is really very useful.

Let’s say we were interested in doing a correlation test between F1 and F2 for all of my vowels (this isn’t an especially interesting thing to do, it’s just for illustration). We could do that like so:

cor.test(vowels_df$F1, vowels_df$F2)

    Pearson's product-moment correlation

data:  vowels_df$F1 and vowels_df$F2
t = -1.9254, df = 327, p-value = 0.05505
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.211579986  0.002279927
sample estimates:
       cor 
-0.1058742 

What if we wanted to do a correlation test, not on all of my vowels together, but within each vowel? Step 1 would be to generalize the line that ran cor.test() above into its own function.

#' @name vowel_cor_test
#' @description conducts a correlation test on the data frame (df) using the columns F1 and F2
vowel_cor_test <- function(df){
  cor.test(df$F1, df$F2)
}

vowel_cor_test(vowels_df)

    Pearson's product-moment correlation

data:  df$F1 and df$F2
t = -1.9254, df = 327, p-value = 0.05505
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.211579986  0.002279927
sample estimates:
       cor 
-0.1058742 

The vowel_cor_test() function we just created is what we want to apply. Next we need to figure out how to split up the data, which we’rem going to do by creating a recursive dataframe:

vowels_nested <- vowels_df %>%
                    group_by(plt_vclass)%>%
                    filter(n() > 5)%>%
                    nest()
vowels_nested

This is the mind bending part. In vowels_nested, the first column is the vowel class, for which there is one row each. In the new data column, the value in each cell is actually a whole data frame.

vowels_nested$data[1]

[[1]]

NA

This means we can now use map() to apply the vowel_cor_test() function to each of these data frames.

vowels_test <- vowels_nested %>%
                  mutate(cor = map(data, vowel_cor_test))
vowels_test

Each item in the cor column is now a correlation test result.

vowels_test$cor[1]

[[1]]

    Pearson's product-moment correlation

data:  df$F1 and df$F2
t = -0.39131, df = 23, p-value = 0.6992
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.4616203  0.3242266
sample estimates:
        cor 
-0.08132266 

We’ve successfully fit a correlation test for every vowel, but the results are kind of trapped away in these strange print outs. We can liberate them with the function tidy() from the broom package. tidy() will convert many statistical tests like this into a dataframe representation.

tidy(vowels_test$cor[[1]])

We actually need to map the tidy() function to the cor column for this to work:

vowels_results <- vowels_test %>%
                    mutate(cor_df = map(cor, tidy))%>%
                    unnest(cor_df)
vowels_results %>%
  arrange(-abs(estimate))

And, as usual, this whole procedure could be chained together into one big analysis.

vowels_df %>%
  group_by(plt_vclass)%>%
  filter(n() > 5)%>%
  nest() %>%
  mutate(cor = map(data, vowel_cor_test),
         cor_df = map(cor, tidy))%>%
  unnest(cor_df)%>%
  arrange(-abs(estimate))