Functional programming in R with Purrr

When you first started in R you likely were writing simple code to generate one outcome.

print("Hello world!")
## [1] "Hello world!"
5 * 6
## [1] 30
x <- c(1, 2, 3, 4, 5)
## [1] 1 2 3 4 5

This is great, you are learning about strings, math, and vectors in R!

Then you get started with some basic analyses. You want to see if you can find the mean of some numbers.

employee <- c('John Doe','Peter Gynn','Jolie Hope')
salary <- c(21000, 23400, 26800)
startdate <- as.Date(c('2010-11-1','2008-3-25','2007-3-14'))

# form dataframe and take mean of salary column
employ_data <- data.frame(employee, salary, startdate)
mean(employ_data$salary)
## [1] 23733.33

Eventually you hopefully get exposed to the tidyverse, and you find how this “opinionated collection of R packages designed for data science” makes data analysis in R easier and more readable!

mtcars %>% 
  group_by(cyl) %>% 
  summarize(mean(mpg))
## # A tibble: 3 x 2
##     cyl `mean(mpg)`
##   <dbl>       <dbl>
## 1  4.00        26.7
## 2  6.00        19.7
## 3  8.00        15.1

Everything is going great! You’ve likely replaced Excel at this point, and potentially SPSS or some other statistical software suite! But then you run into a problem where you need to use a function repeatedly.

You could use something like the following code to calculate one-way ANOVAs for some dependent variables and a set independent variable:

aov_mpg <- aov(mpg ~ factor(cyl), data = mtcars)
summary(aov_mpg)

aov_disp <- aov(disp ~ factor(cyll), data = mtcars)
summary(aov_disp)

aov_hp <- aov(hp ~ factor(cyl), data = mrcars)
summry(aov_hpp)

aov_wt <- aov(wt ~ factor(cyl), datas = mtcars)
summary(aov_wt)

But you copy-pasted code 3x, and oops you made some minor misspelling mistakes which throws an error! (The above code leads to errors!)

Also, what if you realized that you wanted to actually run these ANOVAs for number of gears instead of number of cylinders? You would have to go back and change the factor(cyl) call to factor(gear) 4x! This is not very efficient, and you’re more likely to end up with mistakes as you have to type everything multiple times!


How about another example.

Let’s calculate the R-squared values for the linear relationship between Weight and Miles per Gallon, according to the number of Cylinders.

I have written code below that does this for 4 cylinder cars from the mtcars dataset. This is a worst case scenario, you know some dplyr code (dplyr::filter), but are not comfortable with the pipe. That’s fine, you accomplish your goal but a lot of coding! You would have to duplicate this code for 6 cylinder and 8 cylinder cars, for even more code…

library(tidyverse)
# create df for 4 cylinder cars
cyl_4 <- filter(mtcars, cyl == 4)
## Warning: package 'bindrcpp' was built under R version 3.4.1
# create a linear model on 4 cyl cars
lm_4 <- lm(mpg ~ wt, data = cyl_4)

# get the summ
lm_4_summary <- summary(lm_4)

# get the r.squared value
lm_4cyl_r_squared <- lm_4_summary["r.squared"]

# check the value
lm_4cyl_r_squared
## $r.squared
## [1] 0.5086326

Alternatively, you could do the same thing with the pipe. A lot less typing, but to do this for all 3 subsets means we have to copy paste multiple times, so if you end up wanting to do this as a linear model of mpg ~ disp in addition to mpg ~ wt, you would have to duplicate the code 3 more times and change it 3 more times. This may not seem like a big deal, but eventually is a huge deal once you start to scale up the code (say 10+ times or 100+ times, etc).

# piped analysis
lm_4cyl_rsquared <- mtcars %>% 
  filter(cyl == 4) %>%
  lm(mpg ~ wt, data = .) %>% 
  summary() %>% 
  .$"r.squared"

#check output
lm_4cyl_r_squared
## $r.squared
## [1] 0.5086326

To solve this issue of minimizing repetition with further replication, we can dive straight into purrr! To read more about purrr Hadley Wickham recommends the iteration chapter from “R for Data Science” or alternatively you can look at the purrr documentation. Lastly, Jenny Bryan has a great purrr tutorial here. You can load purrr by itself, but it is also loaded as part of the tidyverse library.

I used to be all meep—meep—PANIC about purrr!!

I used to be all meep—meep—PANIC about purrr!!

now I’m all like map %>% map %>% PARTY!

now I’m all like map %>% map %>% PARTY!

purrr allows you to map functions to data. Appropriately the basic function in purrr is called map()! The map functions transform their input by applying a function to each element and returning a vector the same length as the input.

The base arguments for map() are: .x - list or atomic vector (logical, integer, double/numeric, and character) .f - function, formula, or atomic vector

Basically map() takes a function (.f) and applies it to data (.x).

Going back to our example of grabbing the R-squared from a linear model, we use the following code with purrr.

mtcars %>%
  split(.$cyl) %>%
  map(~ lm(mpg ~ wt, data = .)) %>%
  map(summary) %>%
  map_dbl("r.squared")
##         4         6         8 
## 0.5086326 0.4645102 0.4229655

This generates an output from all 3 of our linear models according to number of cylinders in 5 lines of code! This is the beauty of purrr, efficient scaling of functions!

Let’s break down our linear model R-squared code.

We take the mtcars dataset, split it into data subsets according to the number of cylinders, apply a linear model of mpg by wt to each subset of data, apply a summary function and then pull out the r.squared value. However, while purrr is readable, we need to cover a few quirks of using it.

mtcars %>%
  split(.$cyl) %>%
  map(~ lm(mpg ~ wt, data = .)) %>%
  map(summary) %>%
  map_dbl("r.squared")
##         4         6         8 
## 0.5086326 0.4645102 0.4229655

For our code here you may have noticed we have a “.” placed twice within the code. This is a placeholder for the data, we can see this below. The “.” indicate the left-hand side data, or in this case mtcars. Our split call splits the mtcars dataframe into 3 dataframes, each stored within a list. This may seem odd, but it allows map to cycle through our 3 dataframes and replicate the lm() function on each of them individually.

# piped version
mtcars %>% 
  split(.$cyl)
# base R version
split(mtcars, mtcars$cyl)
## $`4`
##                 mpg cyl  disp  hp drat    wt  qsec vs am gear carb
## Datsun 710     22.8   4 108.0  93 3.85 2.320 18.61  1  1    4    1
## Merc 240D      24.4   4 146.7  62 3.69 3.190 20.00  1  0    4    2
## Merc 230       22.8   4 140.8  95 3.92 3.150 22.90  1  0    4    2
## Fiat 128       32.4   4  78.7  66 4.08 2.200 19.47  1  1    4    1
## Honda Civic    30.4   4  75.7  52 4.93 1.615 18.52  1  1    4    2
## Toyota Corolla 33.9   4  71.1  65 4.22 1.835 19.90  1  1    4    1
## Toyota Corona  21.5   4 120.1  97 3.70 2.465 20.01  1  0    3    1
## Fiat X1-9      27.3   4  79.0  66 4.08 1.935 18.90  1  1    4    1
## Porsche 914-2  26.0   4 120.3  91 4.43 2.140 16.70  0  1    5    2
## Lotus Europa   30.4   4  95.1 113 3.77 1.513 16.90  1  1    5    2
## Volvo 142E     21.4   4 121.0 109 4.11 2.780 18.60  1  1    4    2
## 
## $`6`
##                 mpg cyl  disp  hp drat    wt  qsec vs am gear carb
## Mazda RX4      21.0   6 160.0 110 3.90 2.620 16.46  0  1    4    4
## Mazda RX4 Wag  21.0   6 160.0 110 3.90 2.875 17.02  0  1    4    4
## Hornet 4 Drive 21.4   6 258.0 110 3.08 3.215 19.44  1  0    3    1
## Valiant        18.1   6 225.0 105 2.76 3.460 20.22  1  0    3    1
## Merc 280       19.2   6 167.6 123 3.92 3.440 18.30  1  0    4    4
## Merc 280C      17.8   6 167.6 123 3.92 3.440 18.90  1  0    4    4
## Ferrari Dino   19.7   6 145.0 175 3.62 2.770 15.50  0  1    5    6
## 
## $`8`
##                      mpg cyl  disp  hp drat    wt  qsec vs am gear carb
## Hornet Sportabout   18.7   8 360.0 175 3.15 3.440 17.02  0  0    3    2
## Duster 360          14.3   8 360.0 245 3.21 3.570 15.84  0  0    3    4
## Merc 450SE          16.4   8 275.8 180 3.07 4.070 17.40  0  0    3    3
## Merc 450SL          17.3   8 275.8 180 3.07 3.730 17.60  0  0    3    3
## Merc 450SLC         15.2   8 275.8 180 3.07 3.780 18.00  0  0    3    3
## Cadillac Fleetwood  10.4   8 472.0 205 2.93 5.250 17.98  0  0    3    4
## Lincoln Continental 10.4   8 460.0 215 3.00 5.424 17.82  0  0    3    4
## Chrysler Imperial   14.7   8 440.0 230 3.23 5.345 17.42  0  0    3    4
## Dodge Challenger    15.5   8 318.0 150 2.76 3.520 16.87  0  0    3    2
## AMC Javelin         15.2   8 304.0 150 3.15 3.435 17.30  0  0    3    2
## Camaro Z28          13.3   8 350.0 245 3.73 3.840 15.41  0  0    3    4
## Pontiac Firebird    19.2   8 400.0 175 3.08 3.845 17.05  0  0    3    2
## Ford Pantera L      15.8   8 351.0 264 4.22 3.170 14.50  0  1    5    4
## Maserati Bora       15.0   8 301.0 335 3.54 3.570 14.60  0  1    5    8

Similarily, the “.” in or first map call is a placeholder for data, but in this case it will cycle through our list of 3 dataframes generated by the previous pipe. You can see that we get a list of 3 lm() outputs, we need to map a summary call to each of these to get access to R-squared.

mtcars %>%
  split(.$cyl) %>%
  map(~ lm(mpg ~ wt, data = .))
## $`4`
## 
## Call:
## lm(formula = mpg ~ wt, data = .)
## 
## Coefficients:
## (Intercept)           wt  
##      39.571       -5.647  
## 
## 
## $`6`
## 
## Call:
## lm(formula = mpg ~ wt, data = .)
## 
## Coefficients:
## (Intercept)           wt  
##       28.41        -2.78  
## 
## 
## $`8`
## 
## Call:
## lm(formula = mpg ~ wt, data = .)
## 
## Coefficients:
## (Intercept)           wt  
##      23.868       -2.192

We next map our summary function to each of the list items to get cleaner outputs with R-squared values. We now have the rest of our statistical output, including p values and R-squared.

mtcars %>%
  split(.$cyl) %>%
  map(~ lm(mpg ~ wt, data = .)) %>%
  map(summary)
## $`4`
## 
## Call:
## lm(formula = mpg ~ wt, data = .)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.1513 -1.9795 -0.6272  1.9299  5.2523 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   39.571      4.347   9.104 7.77e-06 ***
## wt            -5.647      1.850  -3.052   0.0137 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.332 on 9 degrees of freedom
## Multiple R-squared:  0.5086, Adjusted R-squared:  0.454 
## F-statistic: 9.316 on 1 and 9 DF,  p-value: 0.01374
## 
## 
## $`6`
## 
## Call:
## lm(formula = mpg ~ wt, data = .)
## 
## Residuals:
##      Mazda RX4  Mazda RX4 Wag Hornet 4 Drive        Valiant       Merc 280 
##        -0.1250         0.5840         1.9292        -0.6897         0.3547 
##      Merc 280C   Ferrari Dino 
##        -1.0453        -1.0080 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   28.409      4.184   6.789  0.00105 **
## wt            -2.780      1.335  -2.083  0.09176 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.165 on 5 degrees of freedom
## Multiple R-squared:  0.4645, Adjusted R-squared:  0.3574 
## F-statistic: 4.337 on 1 and 5 DF,  p-value: 0.09176
## 
## 
## $`8`
## 
## Call:
## lm(formula = mpg ~ wt, data = .)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.1491 -1.4664 -0.8458  1.5711  3.7619 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  23.8680     3.0055   7.942 4.05e-06 ***
## wt           -2.1924     0.7392  -2.966   0.0118 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.024 on 12 degrees of freedom
## Multiple R-squared:  0.423,  Adjusted R-squared:  0.3749 
## F-statistic: 8.796 on 1 and 12 DF,  p-value: 0.01179

Our last map is a bit different. You can see we use map_dbl this time. This indicates we want our output to be a dbl or numeric outcome. We get nice named numbers!

mtcars %>%
  split(.$cyl) %>%
  map(~ lm(mpg ~ wt, data = .)) %>%
  map(summary) %>%
  map_dbl("r.squared")
##         4         6         8 
## 0.5086326 0.4645102 0.4229655

If we had not indicated map_dbl, but instead used map we would get a list of the same outcome.

mtcars %>%
  split(.$cyl) %>% 
  map(~ lm(mpg ~ wt, data = .)) %>%
  map(summary) %>%
  map("r.squared")
## $`4`
## [1] 0.5086326
## 
## $`6`
## [1] 0.4645102
## 
## $`8`
## [1] 0.4229655

You could also use map_dfr which binds the outputs into rows of a dataframe.

mtcars %>%
  split(.$cyl) %>% 
  map(~ lm(mpg ~ wt, data = .)) %>%
  map(summary) %>%
  map_dfr("r.squared")
## # A tibble: 1 x 3
##     `4`   `6`   `8`
##   <dbl> <dbl> <dbl>
## 1 0.509 0.465 0.423

There are limitless applications of purrr and other functions within purrr that greatly empower your functional programming in R. I hope that this guide motivates you to add purrr to your toolbox and explore this useful tidyverse package!

As a brief teaser to some more applications of purrr, I’ll leave you with this example. I mentioned calculating ANOVAs across multiple variables at the beginning. Break down this example on your own and see what you think! (You can copy paste this code into R, but need to load the tidyverse and broom packages first).

mtcars %>%
  mutate(cyl = factor(cyl)) %>%
  select(mpg, disp, hp) %>%
  map(~ aov(.x ~ cyl, data = mtcars)) %>%
  map_dfr(~ broom::tidy(.), .id = 'source') %>%
  mutate(p.value = round(p.value, 5))
##   source      term df       sumsq       meansq statistic p.value
## 1    mpg       cyl  1    817.7130    817.71295  79.56103       0
## 2    mpg Residuals 30    308.3342     10.27781        NA      NA
## 3   disp       cyl  1 387454.0926 387454.09261 130.99888       0
## 4   disp Residuals 30  88730.7021   2957.69007        NA      NA
## 5     hp       cyl  1 100984.1721 100984.17209  67.70993       0
## 6     hp Residuals 30  44742.7029   1491.42343        NA      NA

In closing, I’d like to thank several #r4ds Online Learning Community members for their help in my personal understanding of purrr: Frank Farach, Michael Kuehn, and Kent Johnson.

If you are interested in joining this community led by Jesse Maegan check out her post here!

Related