tidyverse-tips

Tidy notes I took from Dave Robinson's screencasts #TidyTuesday

View on GitHub

Motivation

Below are some notes I have taken on David Robinson’s screencasts, with tips and tricks I use for my own R peregrinations in the Tidyverse framework. Hopefully, these notes will be useful to others. I assume you know the basics of the Tidyverse. If not, check out the quick and dirty introduction to the Tidyverse I wrote some time ago.

It took me some time to fully transition from base R to the Tidyverse. I really clicked when I started watching screencasts by David Robinson on his Youtube channel.

Each week, he goes live for an hour or so, and does exploratory analyses in R with data he has never seen before (!). The data come from the #TidyTuesday project brought to us by the R for Data Science Online Learning Community. This is an awesome initiative, check out the beautiful #tidytuesday visualisations on Twitter.

David Robinson shares the #RStats code he writes during his screencasts, and provides annotations for us to grasp at once what the screencast is about, and what the main steps of the analyses are (including the R functions he uses).

Below are my personal notes, with an obvious bias towards the tricks I find useful for my own work. Check out David Robinson’s talk Ten Tremendous Tricks for Tidyverse. You will also be interested in Emily Robinson’s talk The Lesser Known Stars of the Tidyverse.

Let’s get to it. First data wrangling, second data visualisation. I will use some of the awesome illustrations by Allison Horst throughout, check out her work here.

Setting the scene

Load the tidyverse suite of packages:

library(tidyverse)
## ── Attaching packages ──────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.2     ✓ purrr   0.3.4
## ✓ tibble  3.0.3     ✓ dplyr   1.0.1
## ✓ tidyr   1.1.1     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.5.0
## ── Conflicts ─────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

By the way, if you need to explore the functions of a package, a trick is to use the autocompletion in RStudio. Just type in ?package_name:: and you should see the functions of the package_name package in a pull-down menu.

Now set the theme you like for data visualisation with ggplot2:

theme_set(theme_light())

For illustration, I will use the starwars dataset that comes with the tidyverse collection of packages:

data("starwars")
starwars_raw <- starwars

Data wrangling

Inspect the data

Use the viewer View() to inspect your data:

starwars_raw %>%
  View()

The tibble format is great to double check the format of your columns:

starwars_raw
## # A tibble: 87 x 14
##    name  height  mass hair_color skin_color eye_color birth_year sex   gender
##    <chr>  <int> <dbl> <chr>      <chr>      <chr>          <dbl> <chr> <chr> 
##  1 Luke…    172    77 blond      fair       blue            19   male  mascu…
##  2 C-3PO    167    75 <NA>       gold       yellow         112   none  mascu…
##  3 R2-D2     96    32 <NA>       white, bl… red             33   none  mascu…
##  4 Dart…    202   136 none       white      yellow          41.9 male  mascu…
##  5 Leia…    150    49 brown      light      brown           19   fema… femin…
##  6 Owen…    178   120 brown, gr… light      blue            52   male  mascu…
##  7 Beru…    165    75 brown      light      blue            47   fema… femin…
##  8 R5-D4     97    32 <NA>       white, red red             NA   none  mascu…
##  9 Bigg…    183    84 black      light      brown           24   male  mascu…
## 10 Obi-…    182    77 auburn, w… fair       blue-gray       57   male  mascu…
## # … with 77 more rows, and 5 more variables: homeworld <chr>, species <chr>,
## #   films <list>, vehicles <list>, starships <list>

chr, int, dbl, fct and list are for characters, integers, doubles, factors and lists respectively.

We are gonna use only a few columns of the original dataset:

starwars <- starwars_raw %>%
  select(name, gender, species, mass, height)

Count, count, count

Crude counting

To get a sense of the data, count stuff and use sort = TRUE to arrange things by decreasing counts:

starwars %>% 
  count(name, sort = TRUE) # does a character appear only once
## # A tibble: 87 x 2
##    name                    n
##    <chr>               <int>
##  1 Ackbar                  1
##  2 Adi Gallia              1
##  3 Anakin Skywalker        1
##  4 Arvel Crynyd            1
##  5 Ayla Secura             1
##  6 Bail Prestor Organa     1
##  7 Barriss Offee           1
##  8 BB8                     1
##  9 Ben Quadinaros          1
## 10 Beru Whitesun lars      1
## # … with 77 more rows
starwars %>% 
  count(gender, sort = TRUE) # 
## # A tibble: 3 x 2
##   gender        n
##   <chr>     <int>
## 1 masculine    66
## 2 feminine     17
## 3 <NA>          4

Filter out missing values before counting:

starwars %>% 
  filter(!is.na(gender)) %>% # filter out missing values
  count(gender, sort = TRUE)
## # A tibble: 2 x 2
##   gender        n
##   <chr>     <int>
## 1 masculine    66
## 2 feminine     17

Count with respect to more than one variable:

starwars %>% 
  count(species, gender, sort = TRUE) # 
## # A tibble: 42 x 3
##    species  gender        n
##    <chr>    <chr>     <int>
##  1 Human    masculine    26
##  2 Human    feminine      9
##  3 Droid    masculine     5
##  4 <NA>     <NA>          4
##  5 Gungan   masculine     3
##  6 Mirialan feminine      2
##  7 Wookiee  masculine     2
##  8 Zabrak   masculine     2
##  9 Aleena   masculine     1
## 10 Besalisk masculine     1
## # … with 32 more rows

Weighted counting

You may want to use the wt argument to get a count weighted by another variable. Compare the call to count() with and without the wt argument in the example below. When wt = mass is used, we compute sum(mass) for each species, otherwise we compute the number of rows in each species:

starwars %>%
  count(species, wt = mass)
## # A tibble: 38 x 2
##    species       n
##    <chr>     <dbl>
##  1 Aleena       15
##  2 Besalisk    102
##  3 Cerean       82
##  4 Chagrian      0
##  5 Clawdite     55
##  6 Droid       279
##  7 Dug          40
##  8 Ewok         20
##  9 Geonosian    80
## 10 Gungan      148
## # … with 28 more rows
starwars %>%
  count(species)
## # A tibble: 38 x 2
##    species       n
##    <chr>     <int>
##  1 Aleena        1
##  2 Besalisk      1
##  3 Cerean        1
##  4 Chagrian      1
##  5 Clawdite      1
##  6 Droid         6
##  7 Dug           1
##  8 Ewok          1
##  9 Geonosian     1
## 10 Gungan        3
## # … with 28 more rows

You can add the counts to your tibble by using add_count(). Check out the n column:

starwars %>%
  add_count(species, wt = mass)
## # A tibble: 87 x 6
##    name               gender    species  mass height     n
##    <chr>              <chr>     <chr>   <dbl>  <int> <dbl>
##  1 Luke Skywalker     masculine Human      77    172 1821.
##  2 C-3PO              masculine Droid      75    167  279 
##  3 R2-D2              masculine Droid      32     96  279 
##  4 Darth Vader        masculine Human     136    202 1821.
##  5 Leia Organa        feminine  Human      49    150 1821.
##  6 Owen Lars          masculine Human     120    178 1821.
##  7 Beru Whitesun lars feminine  Human      75    165 1821.
##  8 R5-D4              masculine Droid      32     97  279 
##  9 Biggs Darklighter  masculine Human      84    183 1821.
## 10 Obi-Wan Kenobi     masculine Human      77    182 1821.
## # … with 77 more rows

Complete counting

When you count with respect to two or more variables, it may happen that you don’t have all combinations. No worries, you’re covered with complete(). Compare the example below without and with the call to complete(). The Aleena species has no feminine representative, but this info is implicit. When complete() is used, a row is added with a NA for feminine gender:

starwars %>% 
  filter(species %in% c('Aleena','Droid')) %>%
  count(species, gender)
## # A tibble: 3 x 3
##   species gender        n
##   <chr>   <chr>     <int>
## 1 Aleena  masculine     1
## 2 Droid   feminine      1
## 3 Droid   masculine     5
starwars %>% 
  filter(species %in% c('Aleena','Droid')) %>%
  count(species, gender) %>%
  complete(species, gender)
## # A tibble: 4 x 3
##   species gender        n
##   <chr>   <chr>     <int>
## 1 Aleena  feminine     NA
## 2 Aleena  masculine     1
## 3 Droid   feminine      1
## 4 Droid   masculine     5

If you don’t like NA, you may use a 0 instead with the fill = list() argument:

starwars %>% 
  filter(species %in% c('Aleena','Droid')) %>%
  count(species, gender) %>%
  complete(species, gender, fill = list(n = 0))
## # A tibble: 4 x 3
##   species gender        n
##   <chr>   <chr>     <dbl>
## 1 Aleena  feminine      0
## 2 Aleena  masculine     1
## 3 Droid   feminine      1
## 4 Droid   masculine     5

Computations across columns

There is a new function in dplyr which allows applying functions to columns of a dataset. In the example below, I compute the mean and standard deviation (storing them in a list) across all columns which have a numeric format (where(is.numeric)), while taking care of the missing values (na.rm = TRUE):

starwars %>%
  summarize(across(where(is.numeric), 
                   list(mean = ~mean(.x, na.rm = TRUE), 
                        sd = ~sd(.x, na.rm = TRUE))))
## # A tibble: 1 x 4
##   mass_mean mass_sd height_mean height_sd
##       <dbl>   <dbl>       <dbl>     <dbl>
## 1      97.3    169.        174.      34.8

Just use list(mean = mean, sd = sd) if you have no missing values. Also check out across(starts_with()) and across(everywhere()) for even more flexibility.

The decade trick

Probably the best of David Robinson’s tricks, he uses the integer division %/% to count things at coarse levels. For example, below is how you would count the number of characters in classes of height of width 10 centimeters.

starwars %>%
  count(height_classes = 10 * (height %/% 10))
## # A tibble: 18 x 2
##    height_classes     n
##             <dbl> <int>
##  1             60     1
##  2             70     1
##  3             80     1
##  4             90     4
##  5            110     1
##  6            120     1
##  7            130     1
##  8            150     3
##  9            160    10
## 10            170    15
## 11            180    20
## 12            190    12
## 13            200     4
## 14            210     2
## 15            220     3
## 16            230     1
## 17            260     1
## 18             NA     6

This trick also works for years count(year = 10 * (year %/% 10)), think of decades for example, hence the ‘decade trick’.

You may tune a bit the output by changing the name of the n column:

starwars %>%
  count(height_classes = 10 * (height %/% 10), 
        name = "class_size")
## # A tibble: 18 x 2
##    height_classes class_size
##             <dbl>      <int>
##  1             60          1
##  2             70          1
##  3             80          1
##  4             90          4
##  5            110          1
##  6            120          1
##  7            130          1
##  8            150          3
##  9            160         10
## 10            170         15
## 11            180         20
## 12            190         12
## 13            200          4
## 14            210          2
## 15            220          3
## 16            230          1
## 17            260          1
## 18             NA          6

Lumping levels together

It has always been a nightmare to deal with factors in R. The package forcats is a game changer. For example, the function fct_lump() is useful to manipulate factors with many levels. You just lump levels together by preserving the n most common values, the other levels are lumped in a new Other level:

starwars %>%
  filter(!is.na(species)) %>%
  count(species = fct_lump(f = species, n = 3))
## # A tibble: 4 x 2
##   species     n
##   <fct>   <int>
## 1 Droid       6
## 2 Gungan      3
## 3 Human      35
## 4 Other      39

If you’d like to get rid of Other, just pipe a filter(species != 'Other').

Pipe a (G)LM

Statistical analyses can easily be piped in your tidyverse workflow. You will find two simple examples below. Check out the broom package which cleans up the messy output of built-in R functions like lm() or glm().

Also, recommended is tidymodels a collection of packages for doing statistical analyses in the machine learning spirit. If you want to learn tidymodels by examples, I strongly recommend Julia Silge’s screencasts and her supervised machine learning course.

Simple linear regression

starwars %>%
  lm(mass ~ height, data = .) %>%
  anova()
## Analysis of Variance Table
## 
## Response: mass
##           Df  Sum Sq Mean Sq F value Pr(>F)
## height     1   29854   29854  1.0404  0.312
## Residuals 57 1635658   28696

Here is the tidy version of it:

starwars %>%
  lm(mass ~ height, data = .) %>%
  broom::tidy()
## # A tibble: 2 x 5
##   term        estimate std.error statistic p.value
##   <chr>          <dbl>     <dbl>     <dbl>   <dbl>
## 1 (Intercept)  -13.8     111.       -0.124   0.902
## 2 height         0.639     0.626     1.02    0.312

Get several summary statistics with function glance(), including the $R^2$ and the AIC:

starwars %>%
  lm(mass ~ height, data = .) %>%
  broom::glance()
## # A tibble: 1 x 12
##   r.squared adj.r.squared sigma statistic p.value    df logLik   AIC   BIC
##       <dbl>         <dbl> <dbl>     <dbl>   <dbl> <dbl>  <dbl> <dbl> <dbl>
## 1    0.0179      0.000696  169.      1.04   0.312     1  -386.  777.  783.
## # … with 3 more variables: deviance <dbl>, df.residual <int>, nobs <int>

Logistic regression

starwars %>%
  mutate(human = if_else(species == 'Human', 1, 0)) %>%
  glm(human ~ height, data = ., family = "binomial") %>%
  summary()
## 
## Call:
## glm(formula = human ~ height, family = "binomial", data = .)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -1.134  -1.025  -0.929   1.345   1.397  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.025166   1.198297  -0.856    0.392
## height       0.003489   0.006715   0.520    0.603
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 104.83  on 77  degrees of freedom
## Residual deviance: 104.55  on 76  degrees of freedom
##   (9 observations deleted due to missingness)
## AIC: 108.55
## 
## Number of Fisher Scoring iterations: 4

Here is the tidy version of it:

starwars %>%
  mutate(human = if_else(species == 'Human', 1, 0)) %>%
  glm(human ~ height, data = ., family = "binomial") %>%
  broom::tidy()
## # A tibble: 2 x 5
##   term        estimate std.error statistic p.value
##   <chr>          <dbl>     <dbl>     <dbl>   <dbl>
## 1 (Intercept) -1.03      1.20       -0.856   0.392
## 2 height       0.00349   0.00672     0.520   0.603

Let’s have a look to the fitted values and residuals with the function augment():

starwars %>%
  mutate(human = if_else(species == 'Human', 1, 0)) %>%
  glm(human ~ height, data = ., family = "binomial") %>%
  broom::augment()
## # A tibble: 78 x 9
##    .rownames human height .fitted .resid .std.resid   .hat .sigma .cooksd
##    <chr>     <dbl>  <int>   <dbl>  <dbl>      <dbl>  <dbl>  <dbl>   <dbl>
##  1 1             1    172  -0.425  1.36       1.37  0.0129   1.17 0.0102 
##  2 2             0    167  -0.443 -0.996     -1.00  0.0135   1.17 0.00445
##  3 3             0     96  -0.690 -0.902     -0.937 0.0747   1.18 0.0219 
##  4 4             1    202  -0.320  1.32       1.33  0.0211   1.17 0.0151 
##  5 5             1    150  -0.502  1.40       1.41  0.0193   1.17 0.0166 
##  6 6             1    178  -0.404  1.35       1.36  0.0130   1.17 0.00999
##  7 7             1    165  -0.450  1.37       1.38  0.0139   1.17 0.0112 
##  8 8             0     97  -0.687 -0.903     -0.938 0.0732   1.18 0.0214 
##  9 9             1    183  -0.387  1.35       1.35  0.0136   1.17 0.0103 
## 10 10            1    182  -0.390  1.35       1.36  0.0135   1.17 0.0102 
## # … with 68 more rows

Miscealleanous

?parse_number
parse_number("$1000")
## [1] 1000
parse_number("1,234,567.78")
## [1] 1234568

Data visualisation

Ordering bars in a bar plot

Compare:

starwars %>%
  filter(!is.na(species)) %>%
  count(species = fct_lump(species, 3)) %>%
  ggplot(aes(x = n, y = species)) + 
  geom_col()

with:

starwars %>%
  filter(!is.na(species)) %>%
  count(species = fct_lump(species, 3)) %>%
  mutate(species = fct_reorder(species, n)) %>%
  ggplot(aes(x = n, y = species)) + 
  geom_col()

We have reordered the species factor by n which is created by the call to count(), and the Droid and Gungan species get ordered adequately. Neat isn’t it?! The function fct_reorder() is another nice feature of the forcats package.

Ordering bars within bar plots

When we use facet_wrap() to get a barplot for each level of a factor, the reordering should apply within the levels of this factor. The fct_reorder() call above won’t work:

starwars %>%
  filter(!is.na(species)) %>%
  count(species = fct_lump(species, 5),
        gender) %>%
  mutate(species = fct_reorder(species, n)) %>%
  ggplot(aes(x = n, y = species)) + 
  geom_col() +
  facet_wrap(vars(gender))

There is exactly what we need in the tidtext package. This is the function reorder_within() which works with scale_y_reordered:

library(tidytext)
starwars %>%
  filter(!is.na(species)) %>%
  count(species = fct_lump(species, 5),
        gender) %>%
  mutate(species = reorder_within(species, n, gender)) %>%
  ggplot(aes(x = n, y = species)) + 
  geom_col() +
  scale_y_reordered() + 
  facet_wrap(vars(gender))

Little tricks I always forget about

Free the scales

The argument scales = "free" is useful when using facet_wrap(). It allows the X and Y axes to have their own scale in each panel. You can choose to have a free scale on the X axis only with scales = "free_x", same thing for the Y axis with scales = "free_y".

Flip coordinates

We used to add a coord_flip() following geom_col() to improve the reading of a bar plot by having the categories on the Y axis. This extra line of code is no longer needed as we can simply permute the variables in the aes().

Titles too long

Also, in a facet_wrap(), the title of each panel might be too long so that it doesn’t read properly. There are two ways to fix that. Either you decrease the font size with a theme(strip.text = element_text(size = 6)) or you truncate the title with a mutate(tr_title_text = str_trunc(title_text, 25).

Log scale

It often makes sense to plot your data using log scales. It is very easy to do in ggplot2 by piping a scale_x_log10() or a scale_y_log10().

Axes format

To improve the reading of your figure, it might be useful to represent the unit of an axis in percentage or display numbers with commas. The scales package is what you need. For example, pipe a scale_y_continuous(labels = scales::percent) to have your Y axis in percentages, or scale_x_continuous(labels = scales::comma) to add commas to the numbers of your X axis.

Legend of a map

To draw maps, Dave Robinson uses the sf package. Check out the official site. In case you’re interested in, I wrote an introduction to GIS and mapping in R using the sf package.

When I use maps, I always forget how to tune the legend. The scale_fill_gradient2() function is great in that it allows defining colors for low and high ends of the gradient (e.g. low = "brown", high = "darkgreen"), and the midpoint of the diverging scale (e.g. midpoint = .1).

Stuff I need to learn

David Robinson uses a bunch of tools that he has perfectly internalised so that he doesn’t have to go back and forth between RStudio and the internet to find how to wrangle this or plot that. I am faaaaaaaaaaar from having his level of expertise. I have identified a few things I need to learn to improve myself.

Regular expressions

I find them boring but regular expressions for describing patterns in strings are very useful when you have to filter rows based on some patterns (str_detect()), remove characters (str_remove()) or separate rows (separate_rows()). Good resources are this book chapter dedicated to strings and the vignette of the stringr package.

Web scraping

When you analyse data, sometimes you need some extra information that can be extracted from the internet. The rvest package helps you scrape information from web pages. The official site is a good starting point.

Animation

It is relatively easy to animate a plot (whatever you have made: boxplots, scatter plots, maps, etc) using the gganimate package. Check out the official site.

Networks

Another great package by Thomas Lin Pedersen is ggraph to visualise data that are structured in relations like networks and trees. Check out the official site.

Dimensional modeling

The cross by functions look terrific to aggregate data by calendar dates. I need to understand how to work with them. Check out the official site and this talk by Kaelen Medeiros.