---
title: "Split-Apply-Combine and Plyr"
author: "Statistical Computing, 36-350"
date: "Tuesday April 3, 2018"
---

Last week: Fitting models to data
===

- Fitting models is critical to both statistical inference and prediction
- Exploratory data analysis is a very good first step and gives you a sense of what you're dealing with before you start modeling
- Linear regression is the most basic modeling tool of all, and one of the most ubiquitous
- `lm()` allows you to fit a linear model by specifying a formula, in terms of column names of a given data frame
- Utility functions `coef()`, `fitted()`, `residuals()`, `summary()`, `plot()`, `predict()` are very handy and should be used over manual access tricks
- Logistic regression is the natural extension of linear regression to binary data; use `glm()` with `family="binomial"` and all the same utility functions
- Generalized additive models add a level of flexibility in that they allow the predictors to have nonlinear effects; use `gam()` and utility functions 

<span style="color: blue">Part I</span>
===

<span style="color: blue">*Split-apply-combine*</span>

Reminder: iterating in R without `for()`
===

We've learned some tools in R for iteration without explicit `for()` loops:

- Indexing with conditionals + vectorization 
- `apply()`: apply a function to rows or columns of a matrix or data frame
- `lapply()`: apply a function to elements of a list or vector
- `sapply()`: same as the above, but simplify the output (if possible)
- `tapply()`: apply a function to levels of a factor vector

Split-apply-combine
===

Today we will learn a general strategy for data analysis/management that can be summmarized in three conceptual steps:

- **Split** whatever data object we have into meaningful chunks
- **Apply** the function of interest to each element in this division
- **Combine** the results into a new object of the desired structure

These are conceptual steps; often the apply and combine steps can be performed for us by a single call to the appropriate function from the `apply()` family

Simple but powerful
===

Does split-apply-combine sound simple? It is, but it's very powerful when combined with the right data structures

- As usual, compared to explicit `for()` loops, often requires far less code
- Fits nicely with our previous idea of top-down function design, simply at a broader scope
- Sets you in the right direction towards learning how to use MapReduce/Hadoop for really, really big data sets

Strikes data set 
===

Data set on 18 countries over 35 years (compiled by Bruce Western, in the Sociology Department at Harvard University). The measured variables:

- `country`, `year`: country and year of data collection
- `strike.volume`: days on strike per 1000 workers
- `unemployment`: unemployment rate
- `inflation`: inflation rate
- `left.parliament`: leftwing share of the goverment
- `centralization`: centralization of unions
- `density`: density of unions

---

```{r}
strikes.df = 
  read.csv("http://www.stat.cmu.edu/~ryantibs/statcomp-S18/data/strikes.csv")
dim(strikes.df) # Since 18 × 35 = 630, some years missing from some countries
head(strikes.df)
```

An interesting question
===

Is there a relationship between a country's ruling party alignment (left versus right) and the volume of strikes?

![washington](http://www.stat.cmu.edu/~ryantibs/statcomp-S18/lectures/washington63.jpg) 
![madison](http://www.stat.cmu.edu/~ryantibs/statcomp-S18/lectures/madison11.jpg)

---

How could we approach this?

- Worst way: by hand, write 18 separate code blocks
- Bad way: explicit `for()` loop, where we loop over countries 
- Best way: split appropriately, then use `sapply()`

Work with just one chunk of data
===

**Step 0**: design a function that will work on one chunk of data

So let's write code to do regression on the data from (say) just Italy

```{r}
strikes.df.italy = strikes.df[strikes.df$country=="Italy", ] # Data for Italy
dim(strikes.df.italy)
head(strikes.df.italy)
```

---

```{r}
italy.lm = lm(strike.volume ~ left.parliament, data=strikes.df.italy)
summary(italy.lm)
```

---

```{r}
plot(strikes.df.italy$left.parliament, strikes.df.italy$strike.volume,
     main="Italy strike volume versus leftwing alignment",
     ylab="Strike volume", xlab="Leftwing alignment")
abline(coef(italy.lm), col=2)
```

Functionalization
===

Now let's turn this into a function

```{r}
my.strike.lm = function(country.df) {
  coef(lm(strike.volume ~ left.parliament, data=country.df))
}
my.strike.lm(strikes.df.italy) # New way
coef(italy.lm) # Old way, same result
```

Split our data into appropriate chunks
===

**Step 1**: split our data into appropriate chunks, each of which can be handled by our function. Here, the function `split()` is often helpful: `split(df, f=my.factor)` splits a data frame `df` into several data frames, defined by constant levels of the factor `my.factor`

So we want to split `strikes.df` into 18 smaller data frames, each of which has the data for just one country

```{r}
strikes.by.country = split(strikes.df, f=strikes.df$country)
class(strikes.by.country) # It's a list
names(strikes.by.country) # It has one element for each country
head(strikes.by.country$Italy) # Same as what we saw before
```

Apply our function and combine the results
===

**Steps 2 and 3**: apply our function to each chunk of data, and combine the results. Here, the functions `lapply()` or `sapply()` are often helpful 

So we want to apply `my.strikes.lm()` to each data frame in `strikes.by.country`. Think about what the output will be from each function call: vector of length 2 (intercept and slope), so we can use `sapply()`

```{r}
strikes.coefs = sapply(strikes.by.country, FUN=my.strike.lm)
strikes.coefs 
```

---

```{r}
# We don't care about the intercepts, only the slopes (2nd row).
# Some are positive, some are negative! Let's plot them:
plot(1:ncol(strikes.coefs), strikes.coefs[2,], xaxt="n",
     xlab="", ylab="Regression coefficient", 
     main="Countrywise labor activity by leftwing score")
axis(side=1, at=1:ncol(strikes.coefs),
     labels=colnames(strikes.coefs), las=2, cex.axis=0.5)
abline(h=0, col="grey")
```

<span style="color: blue">Part II</span>
===

<span style="color: blue">*Plyr: `a*ply()` and `l*aply()`*</span>

The `plyr` package
===

Most popular R package of all time (most downloads): `plyr`. This is true for a good reason! Provides extremely useful family of apply-like functions

- Advantage over the built-in `apply()` family is its consistency
- All `plyr` functions are of the form `**ply()`
- Replace `**` with characters denoting types:
    * First character: input type, one of `a`, `d`, `l` 
    * Second character: output type, one of `a`, `d`, `l`, or `_` (drop)

`a*ply()`: input is an array
===

The signature for all `a*ply()` functions is:

```
a*ply(.data, .margins, .fun, ...)
```

- `.data` : an array
- `.margins` : index (or indices) to split the array by
- `.fun` : the function to be applied to each piece
- `...` : additional arguments to be passed to the function

Note that this resembles:

```
apply(X, MARGIN, FUN, ...)
```

Examples of `a*ply()`
===

```{r}
my.array = array(1:27, c(3,3,3))
rownames(my.array) = c("Curly", "Larry", "Moe")
colnames(my.array) = c("Groucho", "Harpo", "Zeppo")
dimnames(my.array)[[3]] = c("Bart", "Lisa", "Maggie")
my.array
```

---

```{r}
library(plyr)
aaply(my.array, 1, sum) # Get back an array
adply(my.array, 1, sum) # Get back a data frame
alply(my.array, 1, sum) # Get back a list
```

---

```{r}
aaply(my.array, 2:3, sum) # Get back a 3 x 3 array
adply(my.array, 2:3, sum) # Get back a data frame
alply(my.array, 2:3, sum) # Get back a list
```

`l*ply()`: input is a list
===

The signature for all `l*ply()` functions is:

```
l*ply(.data, .fun, ...)
```

- `.data` : a list
- `.fun` : the function to be applied to each element
- `...` : additional arguments to be passed to the function

Note that this resembles:

```
lapply(X, FUN, ...)
```

Examples of `l*ply()`
===

```{r}
my.list = list(nums=rnorm(1000), lets=letters, pops=state.x77[,"Population"])
laply(my.list, range) # Get back an array
ldply(my.list, range) # Get back a data frame
llply(my.list, range) # Get back a list
```

---

```{r, error=TRUE}
laply(my.list, summary) # Doesn't work! Outputs have different types/lengths
ldply(my.list, summary) # Doesn't work! Outputs have different types/lengths
llply(my.list, summary) # Works just fine
```

The fourth option for `*`
===

The fourth option for `*` is `_`: the function `a_ply()` (or `l*ply()`) has no explicit return object, but still runs the given function over the given array (or list), possibly producing side effects

```{r}
par(mfrow=c(3,3), mar=c(4,4,1,1))
a_ply(my.array, 2:3, plot, ylim=range(my.array), pch=19, col=6)
```

<span style="color: blue">Part III</span>
===

<span style="color: blue">*Plyr: `d*ply()`*</span>


`d*ply()` : the input is a data frame
===

The signature for all `d*ply()` functions is:

```
d*ply(.data, .variables, .fun, ...)
```

- `.data` : a data frame
- `.variables` : variable (or variables) to split the data frame by
- `.fun` : the function to be applied to each piece
- `...` : additional arguments to be passed to the function

Note that this resembles:

```
tapply(X, INDEX, FUN, ...)
```

Strikes data set, revisited
===

```{r}
# Regression coefficients separately for each country, old way:
strikes.list = split(strikes.df, f=strikes.df$country)
strikes.coefs = sapply(strikes.list, my.strike.lm)
head(strikes.coefs)
```

---

```{r}
# Getting regression coefficient separately for each country, new way:
strikes.coefs.a = daply(strikes.df, .(country), my.strike.lm)
head(strikes.coefs.a) # Get back an array, note the difference to sapply()
```

---

```{r}
strikes.coefs.d = ddply(strikes.df, .(country), my.strike.lm)
head(strikes.coefs.d) # Get back a data frame
```
  
---

```{r}
strikes.coefs.l = dlply(strikes.df, .(country), my.strike.lm)
head(strikes.coefs.l) # Get back a list
```

Splitting on two (or more) variables
===

The function `d*ply()` makes it very easy to split on two (or more) variables: we just specify them, separated by a "," in the `.variables` argument

```{r}
# First create a variable that indicates whether the year is pre 1975, and add
# it to the data frame
strikes.df$yearPre1975 = strikes.df$year <= 1975
# Then use (say) ddply() to compute regression coefficients for each country 
# pre and post 1975
strikes.coefs.1975 = ddply(strikes.df, .(country, yearPre1975), my.strike.lm)
dim(strikes.coefs.1975) # Note that there are 18 x 2 = 36 rows
head(strikes.coefs.1975)
```

---

```{r}
# We can also create factor variables on-the-fly with I(), as we've seen before
strikes.coefs.1975 = ddply(strikes.df, .(country, I(year<=1975)), my.strike.lm)
dim(strikes.coefs.1975) # Again, there are 18 x 2 = 36 rows
head(strikes.coefs.1975)
```

Parallelization
===

What happens if we have a really large data set and we want to use split-apply-combine?

If the individual tasks are unrelated, then we should be speed up the computation by performing them **in parallel**

The `plyr` functions make this quite easy: let's take a look at the full signature for `daply()`:

```
daply(.data, .variables, .fun = NULL, ..., .progress = "none",
  .inform = FALSE, .drop_i = TRUE, .drop_o = TRUE, .parallel = FALSE,
  .paropts = NULL)
```

The second to last argument `.parallel` (default FALSE) is for parallelization. If set to TRUE, then it performs the individual tasks in parallel, using the `foreach` package

The last argument `.paropts` is for more advanced parallelization, these are additional arguments to be passed to `foreach`

---

For more, read the `foreach` package first. May take some time to set up the parallel backend (this is often system specific)

But once set up, parallelization is simple and beautiful with `**ply()`! The difference is just, e.g., 

```
daply(strikes.df, .(country), my.strike.lm)
```

versus 

```
daply(strikes.df, .(country), my.strike.lm, .parallel=TRUE)
```

Summary
===

- The split-apply-combine framework is very helpful for breaking up large data analysis tasks, on ragged data, into manageable chunks 
    * Step 0: design a function to work on one chunk
    * Step 1: split your data into appropriate chunks
    * Steps 2 and 3: apply your function on each chunk and combine the results (these steps can often be done together)
- Apply functions are very important in the context of the split-apply-combine paradigm
- The `plyr` package makes life easier by making the input and outputs types (array, data frame, or list?) as explicit as possible
    * `a*ply()`: input is an array, `*` specifies output type
    * `l*ply()`: input is a list, `*` specifies output type
    * `d*ply()`: input is a data frame, `*` specifies output type