3.4.1 Quick graphs

We can use built-in grouping features of ggplot2 to inspect the data visually without having to pre-manipulate it much.

First, let’s look at a histogram of weekly avocado sales, across all weeks, and all metropolitan areas. How many avocados are people buying in a given week? I’m also going to split this histogram to get two different filled histograms for conventional and organic avocado sales so we can see those sales separately. I’ll set the colors to… ah, avocado-themed.

avo_colors = c("sienna4", "olivedrab2")

plot_hist_avo_total_volume <- avocado %>%
  ggplot(aes(x = total_volume, fill = type)) +
  geom_histogram(position = "identity", bins = 25, alpha = 0.5) +
  scale_fill_manual(values = avo_colors) +
  theme_bw()

plot_hist_avo_total_volume

Ooh, yikes, the weekly total avocados sold looks super gamma-skewed. Most weeks it’s maybe a few tens of thousands of avocados sold, but in some weeks and some regions it gets into the millions of avocados sold (that’s a lot of avocados in a week!), so we would probably do well to log-transform the data to try to get it to take a more normal-ish distribution.

Right now, we can use one of ggplot2’s scale helper functions to set the scale to be on log-10 units.

plot_hist_avo_total_volume +
  scale_x_log10()

That looks much better! Now we can proceed.

Next, let’s break up the graph by region, so we can see the sub-histograms of total avocados sold for each metropolitan area. A super-fast way to break up your graph by each grouping unit is to facet your ggplot, or to break it up into small multiples. I like to use facet_wrap(~ grouping_unit) to break a single plot up, so that the data for each metro area is displayed on its own sub-plot, or facet. The “wrap” part of facet_wrap() tries to arrange the facets in “reading” order, from left-to-right and then top-to-bottom to fill the plot area as efficiently as possible.

plot_hist_avo_total_volume +
  # let's not forget to log-scale the x-axis!
  scale_x_log10() +
  facet_wrap(~ region)

Incidentally, graphing data at first can be a great way of identifying important variables in the data that need to be considered. For example, what if I had forgotten to break up the sales data by conventional vs. organic avocado type?

avocado %>%
  # Notice that I've removed the fill aesthetic
  # so now it's not splitting by avo type
  ggplot(aes(x = total_volume)) +
  geom_histogram(position = "identity", bins = 25, alpha = 0.5) +
  scale_x_log10() +
  theme_bw()

I might see a bimodal histogram and say “huh, that’s funny.” If I weren’t so careful, I might continue my business without accounting for this bimodality. But optimally I’d stop and think “what additional variable might be driving this bimodality?” and after some exploration, realize that I need to consider conventional vs. organic avocado type when looking at avocado sales in order to properly characterize the data. Always keep an eye out!

3.4.2 Basic summary statistics

In addition to graphs, we can generate simple one-shot summary statistics using group_by() and summarize() from the dplyr package.

As a companion to the exploratory graphs of total_volume we just generated, we can calculate numerical summaries of the weekly number of avocados sold as a quick metric of the difference in buying patterns for conventional and organic avocados. Before we calculate these summaries, though, we should create a new column for log-10-transformed total_volume, as we saw earlier in our graphs that it’s massively right-tailed and needs to be log-transformed to look more normal.

# The double-sided assignment pipe %<>%, from magrittr
# you can use this by library()-ing magrittr
# in addition to library(tidyverse)
# avocado %<>% ... does the same thing as:
# avocado <- avocado %>% ...
# so it gives us a shortcut to pipe an object into some calls,
# and then re-save the new result back into the old variable name
# let's save total_volume_log10 into avocado because we'll need it
# many times in the future
avocado %<>%
  mutate(total_volume_log10 = log10(total_volume))

avocado %>%
  group_by(region, type) %>%
  summarize(mean_total_volume = mean(total_volume_log10),
            sd_total_volume = sd(total_volume_log10))

## `summarise()` has grouped output by 'region'. You can override using the `.groups` argument.

## # A tibble: 92 x 4
## # Groups:   region [46]
##    region              type         mean_total_volume sd_total_volume
##    <chr>               <chr>                    <dbl>           <dbl>
##  1 Albany              conventional              4.95          0.140 
##  2 Albany              organic                   3.29          0.205 
##  3 Atlanta             conventional              5.70          0.0899
##  4 Atlanta             organic                   4.00          0.236 
##  5 BaltimoreWashington conventional              5.88          0.0722
##  6 BaltimoreWashington organic                   4.31          0.225 
##  7 Boise               conventional              4.91          0.0987
##  8 Boise               organic                   3.34          0.209 
##  9 Boston              conventional              5.74          0.0811
## 10 Boston              organic                   4.03          0.359 
## # … with 82 more rows

Above, we group_by() both region and type to get two separate means and SDs of weekly avocados sold in each region, one value each for conventional and organic avocados.

We can pivot_wider() our observation-level data, to get total_volume into two columns, one for the conventional avocados and one for the organic avocados, for each week. Thus, instead of two rows per week per metro region, one for each avocado type, we’ll have one row per week per metro region, with the relevant values for each avocado type in their own columns.

avocado %>%
  # Here, selecting the id_cols and the ONE values_from column I want to pivot out
  select(year, date, region, type, total_volume_log10) %>%
  pivot_wider(names_from = type,
              values_from = total_volume_log10,
              names_prefix = "total_volume_")

## # A tibble: 7,774 x 5
##     year date       region total_volume_conventional total_volume_organic
##    <int> <date>     <chr>                      <dbl>                <dbl>
##  1  2015 2015-12-27 Albany                      4.81                 3.00
##  2  2015 2015-12-20 Albany                      4.74                 3.07
##  3  2015 2015-12-13 Albany                      5.07                 3.00
##  4  2015 2015-12-06 Albany                      4.90                 3.06
##  5  2015 2015-11-29 Albany                      4.71                 2.92
##  6  2015 2015-11-22 Albany                      4.75                 2.93
##  7  2015 2015-11-15 Albany                      4.92                 3.08
##  8  2015 2015-11-08 Albany                      5.04                 3.12
##  9  2015 2015-11-01 Albany                      5.00                 3.01
## 10  2015 2015-10-25 Albany                      4.87                 3.07
## # … with 7,764 more rows

Instructor’s note: Above, I’ve select()-ed to remove all the other outcome columns other than total_volume before pivoting my dataframe. If you wanted to pivot all outcome columns in the data by type, you could, as of tidyr >=1.0.0. The least verbose way to do this would be to feed a negative column name vector into the values_from argument to exclude ID columns from the set of columns to pivot, like so:

avocado %>%
  pivot_wider(names_from = type, values_from = -c(year, date, region))

## # A tibble: 7,774 x 11
##    date        year region avg_price_convent… avg_price_organ… total_volume_con…
##    <date>     <int> <chr>               <dbl>            <dbl>             <dbl>
##  1 2015-12-27  2015 Albany               1.33             1.83            64237.
##  2 2015-12-20  2015 Albany               1.35             1.89            54877.
##  3 2015-12-13  2015 Albany               0.93             1.85           118220.
##  4 2015-12-06  2015 Albany               1.08             1.84            78992.
##  5 2015-11-29  2015 Albany               1.28             1.94            51040.
##  6 2015-11-22  2015 Albany               1.26             1.94            55980.
##  7 2015-11-15  2015 Albany               0.99             1.89            83454.
##  8 2015-11-08  2015 Albany               0.98             1.88           109428.
##  9 2015-11-01  2015 Albany               1.02             1.88            99811.
## 10 2015-10-25  2015 Albany               1.07             1.83            74339.
## # … with 7,764 more rows, and 5 more variables: total_volume_organic <dbl>,
## #   type_conventional <chr>, type_organic <chr>,
## #   total_volume_log10_conventional <dbl>, total_volume_log10_organic <dbl>

Then, by default, since there are multiple values_from columns, each of them has the relevant level of names_from (here, avocado type) appended to the end of the column name with the names_sep separator.

Once the data is pivoted, we can now use mutate() to subtract the organic total dollars spent from the conventional total dollars spent. Then, we can generate similar mean/SD summaries as we did before. Finally, we can make the little tibble print preview a bit more informative by using arrange() to sort the data such that the metro areas with highest conventional - organic sales volume differences get sorted to the top of the tibble.

avocado %>%
  select(year, date, region, type, total_volume_log10) %>%
  pivot_wider(names_from = type,
              values_from = total_volume_log10,
              names_prefix = "total_volume_") %>%
  mutate(total_volume_diff = total_volume_conventional - total_volume_organic) %>%
  group_by(region) %>%
  summarize(mean_volume_diff = mean(total_volume_diff),
            sd_volume_diff = sd(total_volume_diff)) %>%
  arrange(desc(mean_volume_diff))

## # A tibble: 46 x 3
##    region            mean_volume_diff sd_volume_diff
##    <chr>                        <dbl>          <dbl>
##  1 MiamiFtLauderdale             2.25          0.285
##  2 Tampa                         2.01          0.228
##  3 PhoenixTucson                 1.99          0.134
##  4 NewOrleansMobile              1.91          0.247
##  5 GrandRapids                   1.84          0.390
##  6 Orlando                       1.84          0.170
##  7 DallasFtWorth                 1.80          0.177
##  8 Houston                       1.76          0.235
##  9 Sacramento                    1.76          0.149
## 10 Jacksonville                  1.74          0.171
## # … with 36 more rows

Huh, cities in Florida appear to be over-represented in the top several metro areas in buying more conventional than organic avocados. I wonder what’s going on there? We can’t rule out that this is driven by Florida having waaaay cheaper conventional avocados, and relatively more expensive organic avocados, which could be driving Floridian shoppers to choose conventional. We probably want to take this into account in later analyses, so that any effects of conventional vs. organic we see on the volume of avocados sold aren’t driven by a huge price premium for organic avocados. (Stay tuned!)

We can get a fair amount of mileage out of these initial visualizations and summaries, and there’s nothing wrong with using these summaries as a first pass. I do so all the time! However, when you want (or need) to be able to calculate more complex statistics on the weekly data for each metro area, you’ll need more complex tools.

3.5.1 Intro to nest(), unnest(), and tibble list-columns

A df is composed of columns, each of which is a vector of identical length. But nobody ever said the columns of a df all have to be atomic vectors! A df column can be a list vector, where each element of said list-column can contain objects of different data types and lengths. While each element of list-column can contain objects of various and sundry types/lengths, list-columns gain a lot of power when they contain objects of consistent type. (When this is true, you can apply a single function to every row of a list-column, and get a single column of consistent output–this is where we’re headed.)

In this way, a list-column allows you to store multiple observations per unit wrapped up in such a way that the main df still has one row per unit, but you retain the observation-level data because you haven’t actually collapsed across any variables.

First, we’ll create a new variable by folding our observation-level data into a new shape, using nest() from the tidyr package. This is going to package the long-form avocado data into nested form, which is somewhere between long data (one row per observation) and wide data (one row per grouping unit).

avo_by_region <- avocado %>%
  nest(sales = -region)

A nested dataframe like the one we just created has two main column types: key columns and list-columns of sub-dataframes.

Key columns are the columns you might group_by() before, say, generating summary statistics with summarize(). They act as a label for the data contained in that row’s list-columns. List-columns do not contain a single atomic element in each row, but instead _a whole sub-data frame containing all the observation-level data corresponding to that key level.

Instructor’s note: A deeper exploration of the uses for and behaviors of list-columns is beyond the scope of the current tutorial. This tutorial assumes learners have never intentionally worked with non-dataframe list objects, and demonstrates only the case of list-columns containing data nested from an original “normal” dataframe. R for Data Science’s section on dataframe list-columns is a useful reference for more information on list-columns.

Now, for list-columns, instead of actually printing the content of the vector, tibble output gives us a blurb about the object contained inside of each list element.

avo_by_region

## # A tibble: 46 x 2
##    region              sales             
##    <chr>               <list>            
##  1 Albany              <tibble [338 × 6]>
##  2 Atlanta             <tibble [338 × 6]>
##  3 BaltimoreWashington <tibble [338 × 6]>
##  4 Boise               <tibble [338 × 6]>
##  5 Boston              <tibble [338 × 6]>
##  6 BuffaloRochester    <tibble [338 × 6]>
##  7 Charlotte           <tibble [338 × 6]>
##  8 Chicago             <tibble [338 × 6]>
##  9 CincinnatiDayton    <tibble [338 × 6]>
## 10 Columbus            <tibble [338 × 6]>
## # … with 36 more rows

We can see here that each row of sales, the list-column generated by nesting the data, contains a sub-tibble with 338 rows of the other 6 columns of data corresponding to each level of region. The data is still there, safe and sound!

nest() creates list-columns with the general argument format list_col = cols_to_put_into_list_col. Importantly, the columns you wish to put into the list-column can be specified using select() syntax. We can harness this to use the - to exclude specific key columns, and put all other non-key columns into the list-column. For example, if you wanted to hold out multiple key columns, for example, both region and type to split up the sub-data by conventional and organic:

avocado %>%
  nest(sales = -c(region, type))

## # A tibble: 92 x 3
##    type         region              sales             
##    <chr>        <chr>               <list>            
##  1 conventional Albany              <tibble [169 × 5]>
##  2 conventional Atlanta             <tibble [169 × 5]>
##  3 conventional BaltimoreWashington <tibble [169 × 5]>
##  4 conventional Boise               <tibble [169 × 5]>
##  5 conventional Boston              <tibble [169 × 5]>
##  6 conventional BuffaloRochester    <tibble [169 × 5]>
##  7 conventional Charlotte           <tibble [169 × 5]>
##  8 conventional Chicago             <tibble [169 × 5]>
##  9 conventional CincinnatiDayton    <tibble [169 × 5]>
## 10 conventional Columbus            <tibble [169 × 5]>
## # … with 82 more rows

In this way, you can flexibly choose which columns of data to fold into your list-column, and which key-columns you want to fold by.

When you want to unnest the data and go back to original long form, you might not be surprised to hear that the function we’ll need is unnest().

avo_by_region %>%
  unnest(sales)

## # A tibble: 15,545 x 7
##    region date       avg_price total_volume type          year total_volume_log…
##    <chr>  <date>         <dbl>        <dbl> <chr>        <int>             <dbl>
##  1 Albany 2015-12-27      1.33       64237. conventional  2015              4.81
##  2 Albany 2015-12-20      1.35       54877. conventional  2015              4.74
##  3 Albany 2015-12-13      0.93      118220. conventional  2015              5.07
##  4 Albany 2015-12-06      1.08       78992. conventional  2015              4.90
##  5 Albany 2015-11-29      1.28       51040. conventional  2015              4.71
##  6 Albany 2015-11-22      1.26       55980. conventional  2015              4.75
##  7 Albany 2015-11-15      0.99       83454. conventional  2015              4.92
##  8 Albany 2015-11-08      0.98      109428. conventional  2015              5.04
##  9 Albany 2015-11-01      1.02       99811. conventional  2015              5.00
## 10 Albany 2015-10-25      1.07       74339. conventional  2015              4.87
## # … with 15,535 more rows

unnest() takes as its argument the list-column you wish to unnest, and then thy will be done, it takes the data back into long form for you. Notice that the key column, in this case region, is now once again repeated for every row in the data belonging to that key level.

To inspect the contents of a single element of a list-column, we have to index into the list-column. The tidyverse function to index single columns, equivalent to dollar-sign $ indexing dataframe columns in base R, is pull(). pull() is pipe-safe, meaning it takes a tibble/dataframe as its first argument, so we can pipe said dataframe into the function. Then, the second argument that we specify inside of the parentheses is the column name we want to index, “naked”, without quotation marks, as with other tidyverse functions.

avo_by_region %>%
  pull(region)

##  [1] "Albany"              "Atlanta"             "BaltimoreWashington"
##  [4] "Boise"               "Boston"              "BuffaloRochester"   
##  [7] "Charlotte"           "Chicago"             "CincinnatiDayton"   
## [10] "Columbus"            "DallasFtWorth"       "Denver"             
## [13] "Detroit"             "GrandRapids"         "GreatLakes"         
## [16] "HarrisburgScranton"  "HartfordSpringfield" "Houston"            
## [19] "Indianapolis"        "Jacksonville"        "LasVegas"           
## [22] "LosAngeles"          "Louisville"          "MiamiFtLauderdale"  
## [25] "Nashville"           "NewOrleansMobile"    "NewYork"            
## [28] "NorthernNewEngland"  "Orlando"             "Philadelphia"       
## [31] "PhoenixTucson"       "Pittsburgh"          "Portland"           
## [34] "RaleighGreensboro"   "RichmondNorfolk"     "Roanoke"            
## [37] "Sacramento"          "SanDiego"            "SanFrancisco"       
## [40] "Seattle"             "SouthCarolina"       "Spokane"            
## [43] "StLouis"             "Syracuse"            "Tampa"              
## [46] "WestTexNewMexico"

We can see below that base R $ indexing yields the same result:

avo_by_region$region

##  [1] "Albany"              "Atlanta"             "BaltimoreWashington"
##  [4] "Boise"               "Boston"              "BuffaloRochester"   
##  [7] "Charlotte"           "Chicago"             "CincinnatiDayton"   
## [10] "Columbus"            "DallasFtWorth"       "Denver"             
## [13] "Detroit"             "GrandRapids"         "GreatLakes"         
## [16] "HarrisburgScranton"  "HartfordSpringfield" "Houston"            
## [19] "Indianapolis"        "Jacksonville"        "LasVegas"           
## [22] "LosAngeles"          "Louisville"          "MiamiFtLauderdale"  
## [25] "Nashville"           "NewOrleansMobile"    "NewYork"            
## [28] "NorthernNewEngland"  "Orlando"             "Philadelphia"       
## [31] "PhoenixTucson"       "Pittsburgh"          "Portland"           
## [34] "RaleighGreensboro"   "RichmondNorfolk"     "Roanoke"            
## [37] "Sacramento"          "SanDiego"            "SanFrancisco"       
## [40] "Seattle"             "SouthCarolina"       "Spokane"            
## [43] "StLouis"             "Syracuse"            "Tampa"              
## [46] "WestTexNewMexico"

There are definitely still instances where $ indexing columns is the cleanest way to index. However, in this tutorial we’ll be using pipe-safe tidyverse indexing functions moving forward so that we can stick with the %>% pipe for consistency.

Once we’ve used pull() to index a single column, we can use another tidyverse pipe-safe helper function, pluck(), this time for indexing elements of a vector. As a rule of thumb, pull() does what $ indexing does, and pluck() does what [] indexing does. Thus, the below code indexes the sales column, and then specifically the first element in that column vector.

avo_by_region %>%
  pull(sales) %>%
  pluck(1)

## # A tibble: 338 x 6
##    date       avg_price total_volume type          year total_volume_log10
##    <date>         <dbl>        <dbl> <chr>        <int>              <dbl>
##  1 2015-12-27      1.33       64237. conventional  2015               4.81
##  2 2015-12-20      1.35       54877. conventional  2015               4.74
##  3 2015-12-13      0.93      118220. conventional  2015               5.07
##  4 2015-12-06      1.08       78992. conventional  2015               4.90
##  5 2015-11-29      1.28       51040. conventional  2015               4.71
##  6 2015-11-22      1.26       55980. conventional  2015               4.75
##  7 2015-11-15      0.99       83454. conventional  2015               4.92
##  8 2015-11-08      0.98      109428. conventional  2015               5.04
##  9 2015-11-01      1.02       99811. conventional  2015               5.00
## 10 2015-10-25      1.07       74339. conventional  2015               4.87
## # … with 328 more rows

We can also use the dplyr tools we already know to subset avo_by_region() to pluck() the data belonging to a particular region.

avo_by_region %>%
  # This should give us a one-row dataframe
  # so we know the first element is the one we want
  filter(region == "NewYork") %>%
  pull(sales) %>%
  pluck(1)

## # A tibble: 338 x 6
##    date       avg_price total_volume type          year total_volume_log10
##    <date>         <dbl>        <dbl> <chr>        <int>              <dbl>
##  1 2015-12-27      1.17     1129876. conventional  2015               6.05
##  2 2015-12-20      1.23     1139348. conventional  2015               6.06
##  3 2015-12-13      1.12     1254805. conventional  2015               6.10
##  4 2015-12-06      1.2      1068972. conventional  2015               6.03
##  5 2015-11-29      1.16      999170. conventional  2015               6.00
##  6 2015-11-22      1.14     1111803. conventional  2015               6.05
##  7 2015-11-15      1.04     1357393. conventional  2015               6.13
##  8 2015-11-08      1.13     1406262. conventional  2015               6.15
##  9 2015-11-01      1.06     2180520. conventional  2015               6.34
## 10 2015-10-25      1.23     1048046. conventional  2015               6.02
## # … with 328 more rows

Notice that once you’ve got only one row in a nested dataframe, pluck()-ing to unlist the data in the list-column is equivalent to using unnest() to unfold all rows of the list-column.

avo_by_region %>%
  filter(region == "NewYork") %>%
  unnest(sales)

## # A tibble: 338 x 7
##    region  date       avg_price total_volume type         year total_volume_log…
##    <chr>   <date>         <dbl>        <dbl> <chr>       <int>             <dbl>
##  1 NewYork 2015-12-27      1.17     1129876. convention…  2015              6.05
##  2 NewYork 2015-12-20      1.23     1139348. convention…  2015              6.06
##  3 NewYork 2015-12-13      1.12     1254805. convention…  2015              6.10
##  4 NewYork 2015-12-06      1.2      1068972. convention…  2015              6.03
##  5 NewYork 2015-11-29      1.16      999170. convention…  2015              6.00
##  6 NewYork 2015-11-22      1.14     1111803. convention…  2015              6.05
##  7 NewYork 2015-11-15      1.04     1357393. convention…  2015              6.13
##  8 NewYork 2015-11-08      1.13     1406262. convention…  2015              6.15
##  9 NewYork 2015-11-01      1.06     2180520. convention…  2015              6.34
## 10 NewYork 2015-10-25      1.23     1048046. convention…  2015              6.02
## # … with 328 more rows

Instructor’s note: To use exclusively base R to index into list-columns, you have to be careful when you want to use single-bracket indexing (e.g. to conditional-index based in values in another column of the df) versus double-bracket indexing (to actually expose the value of the list-column for printing to console or otherwise manipulation). For a more in-depth exploration of base R indexing into lists, please visit the relevant R for Data Science chapter.

3.6.1 Vectorizing: it’s like a for loop, but not

As you’ve likely experienced, much data processing requires writing code to perform a particular action, and then wrapping it in other code to repeat that action many times over many objects. In R, this usually means doing the same operation to every element in a vector or every row in a df column.

In most computing languages, the most sensible way to do one thing a bunch of times is to use a for loop, which (I know you know this but I’ll spell it out just to be extra) repeats the code inside the loop once for every element of a vector that’s specified at the beginning of the loop code.

To use a for loop to apply some function to every element of a df column, you might write something like:

for (i in 1:nrow(df)) {
  df$new_col[i] = my_function(df$old_col[i])
}

and that would do the trick. In this case, the vector you iterate along contains the row indices of your df, which means that on every iteration of your loop you’re calling my_function() on the next row of data until you get to the end of your df.

This is a perfectly functional way to write code! Today, I’ll argue that it’s prone to typos that might cause you big headaches. Specifically, these typos may cause your code to fail silently (ominous organ riff), or do something unintended without throwing an error. This means you wouldn’t find out that something was wrong unless you visually inspected the output, which you can’t always do after every command (I getcha). Thus, there’s a risk of inducing mistakes in the data that you don’t catch until it’s too late.

For example, once I wrote a for loop just like the above, to perform some functions on every row of a df column, but I wrote df$old_col[1] instead of df$old_col[i]. I accidentally produced the same value in every element of df$new_col, because even though the for loop was iterating as usual, on every iteration of the loop it was calling the same value. I didn’t discover my typo until weeks later. No bueno!

Enter… vectorizing.

It turns out that R has built-in optimized functionality for repeating functions along every element of a vector, just waiting for you to harness!

Generally, any vectorized function:

• takes a vector as input
• does the same thing to every element of that vector
• returns a vector of the same length as output

Plenty of functions are already vectorized, and you likely already use them as such, including:

• math operators (R assumes element-wise math as the default, and element-wise == vectorized here)
  • arithmetic operators: +, -, and the like
  • round(), for example
• string manipulation functions
  • paste() always returns a vector equal to the longest input vector
  • pattern matching functions like grep(), grepl(), and such
• statistical distribution functions
  • the d*(), p*(), and q*() statistical distribution functions (e.g. qnorm()) are all vectorized along their main input arguments
• other stuff too!
  • base::ifelse() is vectorized, as it determines true/false for each pair of elements in the two vector arguments element-wise
  • and MANY more

In general, it’s worth always checking to see if the function you plan to call on a vector is vectorized, and then you can use it on your vector or dataframe column without needing any additional looping.

3.6.2 Tidy vectorizing

If you use mutate() for your everyday column-wise data manipulation needs, using vectorized functions is smooth. In fact, mutate() is built to encourage you to use vectorized functions for column manipulation.

If we go back to the Hass avocado sale data, we can see this at work. For example, let’s take another look at the code we used earlier to create a new column for the weekly total number of avocados sold, but log-10 transformed. Earlier, we saw order-of-magnitude differences in avocado sales between different regions, so we log-transformed the data to get it to be more normally distributed.

avocado %>%
  mutate(total_volume_log10 = log10(total_volume)) %>%
  select(region, date, total_volume, total_volume_log10)

## # A tibble: 15,545 x 4
##    region date       total_volume total_volume_log10
##    <chr>  <date>            <dbl>              <dbl>
##  1 Albany 2015-12-27       64237.               4.81
##  2 Albany 2015-12-20       54877.               4.74
##  3 Albany 2015-12-13      118220.               5.07
##  4 Albany 2015-12-06       78992.               4.90
##  5 Albany 2015-11-29       51040.               4.71
##  6 Albany 2015-11-22       55980.               4.75
##  7 Albany 2015-11-15       83454.               4.92
##  8 Albany 2015-11-08      109428.               5.04
##  9 Albany 2015-11-01       99811.               5.00
## 10 Albany 2015-10-25       74339.               4.87
## # … with 15,535 more rows

We can write this in one call, without needing to wrap it in a loop, because log10() is vectorized, and returns a vector the same length as the input vector. Neat! You can use any vectorized functions “naked” (without any extra looping code) to create new dataframe columns inside of mutate() with no problem.

But what about functions that aren’t already vectorized? Can we use R magic to make them vectorized, if we so wish? Why yes, we can!

avo_by_region %>%
  # Here, I'm using slice() to get only the first 10 rows of the tibble
  # to shorten the eventual output. Otherwise it would output the nrows
  # for every single region in the avocado data
  slice(1:10) %>%
  pull(sales) %>%
  map(~nrow(.x))

## [[1]]
## [1] 338
## 
## [[2]]
## [1] 338
## 
## [[3]]
## [1] 338
## 
## [[4]]
## [1] 338
## 
## [[5]]
## [1] 338
## 
## [[6]]
## [1] 338
## 
## [[7]]
## [1] 338
## 
## [[8]]
## [1] 338
## 
## [[9]]
## [1] 338
## 
## [[10]]
## [1] 338

avo_by_region %>%
  slice(1:10) %>%
  pull(sales) %>%
  map_int(~nrow(.x))

##  [1] 338 338 338 338 338 338 338 338 338 338

avo_by_region_type <- avocado %>%
  # Notice the -c() select() syntax to exclude TWO key columns
  # from the list-column
  nest(sales = -c(region, type))

avo_by_region_type

## # A tibble: 92 x 3
##    type         region              sales             
##    <chr>        <chr>               <list>            
##  1 conventional Albany              <tibble [169 × 5]>
##  2 conventional Atlanta             <tibble [169 × 5]>
##  3 conventional BaltimoreWashington <tibble [169 × 5]>
##  4 conventional Boise               <tibble [169 × 5]>
##  5 conventional Boston              <tibble [169 × 5]>
##  6 conventional BuffaloRochester    <tibble [169 × 5]>
##  7 conventional Charlotte           <tibble [169 × 5]>
##  8 conventional Chicago             <tibble [169 × 5]>
##  9 conventional CincinnatiDayton    <tibble [169 × 5]>
## 10 conventional Columbus            <tibble [169 × 5]>
## # … with 82 more rows

avo_by_region_type %>%
  pull(sales) %>%
  pluck(1)

## # A tibble: 169 x 5
##    date       avg_price total_volume  year total_volume_log10
##    <date>         <dbl>        <dbl> <int>              <dbl>
##  1 2015-12-27      1.33       64237.  2015               4.81
##  2 2015-12-20      1.35       54877.  2015               4.74
##  3 2015-12-13      0.93      118220.  2015               5.07
##  4 2015-12-06      1.08       78992.  2015               4.90
##  5 2015-11-29      1.28       51040.  2015               4.71
##  6 2015-11-22      1.26       55980.  2015               4.75
##  7 2015-11-15      0.99       83454.  2015               4.92
##  8 2015-11-08      0.98      109428.  2015               5.04
##  9 2015-11-01      1.02       99811.  2015               5.00
## 10 2015-10-25      1.07       74339.  2015               4.87
## # … with 159 more rows

avo_by_region_type %>%
  mutate(overall_avg_price = map_dbl(sales, ~.x %>%
                                       pull(avg_price) %>%
                                       mean()
                                     )
         )

## # A tibble: 92 x 4
##    type         region              sales              overall_avg_price
##    <chr>        <chr>               <list>                         <dbl>
##  1 conventional Albany              <tibble [169 × 5]>              1.35
##  2 conventional Atlanta             <tibble [169 × 5]>              1.07
##  3 conventional BaltimoreWashington <tibble [169 × 5]>              1.34
##  4 conventional Boise               <tibble [169 × 5]>              1.08
##  5 conventional Boston              <tibble [169 × 5]>              1.30
##  6 conventional BuffaloRochester    <tibble [169 × 5]>              1.38
##  7 conventional Charlotte           <tibble [169 × 5]>              1.28
##  8 conventional Chicago             <tibble [169 × 5]>              1.37
##  9 conventional CincinnatiDayton    <tibble [169 × 5]>              1.02
## 10 conventional Columbus            <tibble [169 × 5]>              1.07
## # … with 82 more rows

avo_by_region_type %>%
  # I'm indenting the close parentheses and such this way
  # to use the indentation levels to visually cue
  # which chunk of code lives in which level
  # because we are pretty deep here!
  mutate(overall_avg_price = map_dbl(sales,
                                     ~.x %>%
                                       pull(avg_price) %>%
                                       mean()
                                     ),
         avg_total_volume = map_dbl(sales,
                                    ~.x %>%
                                      pull(total_volume) %>%
                                      mean()
                                    )
         )

## # A tibble: 92 x 5
##    type        region           sales          overall_avg_pri… avg_total_volume
##    <chr>       <chr>            <list>                    <dbl>            <dbl>
##  1 convention… Albany           <tibble [169 …             1.35           92903.
##  2 convention… Atlanta          <tibble [169 …             1.07          512789.
##  3 convention… BaltimoreWashin… <tibble [169 …             1.34          773642.
##  4 convention… Boise            <tibble [169 …             1.08           82843.
##  5 convention… Boston           <tibble [169 …             1.30          561541.
##  6 convention… BuffaloRochester <tibble [169 …             1.38          130410.
##  7 convention… Charlotte        <tibble [169 …             1.28          203331.
##  8 convention… Chicago          <tibble [169 …             1.37          759816.
##  9 convention… CincinnatiDayton <tibble [169 …             1.02          247835.
## 10 convention… Columbus         <tibble [169 …             1.07          169268.
## # … with 82 more rows

3.6.3 A caveat: When not to vectorize

Ultimately, I love vectorizing functions like map()–they’ve helped me keep a lot more of my data processing inside of tidyverse functions, and often saved me from copying and pasting code. There are, of course, cases when writing your code to work inside of vectorizing helpers isn’t necessarily optimal. Here are a couple that I’ve run into:

• If you need to recursively reference earlier vector elements while looping
  • Vectorizer functions run every iteration of your vectorizing loop independently, which means that the nth iteration cannot access the output of the (n-1)th iteration or any previous ones. If you need to be able to access previous loop content in later iterations, a for loop is the way to go.
• If you really need a progress bar to print to console, for iterating over looooong vectors
  • True console-based progress bars with must recursively access the console to create a progress bar that stays on one line from 0% to 100%, so they work best inside for loops
  • The good folks at RStudio who maintain purrr are still working on building automatic progress bars for map(), but as of publishing they’re not implemented
  • Check out the progress package, which contains helper functions that work inside for loops to make your own progress bars

3.7.1 Fitting many models

We can use dataframe list-columns, manipulated with map() inside mutate(), to fit this particular model to the avocado sales data from every metropolitan area. Since mutate() is amenable to defining a flexible number of columns, we can first do some light data preprocessing on the observation-level data contained in the sales column, and then in the next command, within mutate(), create a new list-column using map() to fit an lm() model to each row’s data in sales.

# Again, the double-sided pipe in magrittr
# to overwrite the old value of avo_by_region with the new value
# with all the new columns in it
avo_by_region %<>%
  mutate(sales = map(sales,
                     # Yes that's right folks, we are about to use
                     # mutate() inside map() inside mutate()
                     ~.x %>%
                       # Here, let's roughly re-center avg_price
                       # so the "baseline", or 0 point, is set at $1.00/avocado
                       # we already calculate total_volume_log10 earlier
                       # so we should be good there
                       mutate(avg_price_shift1 = avg_price - 1)
                     ),
         # Now that we've preprocesed our columns
         # we can fit our multiple regression
         # using lm() inside of map()
         model = map(sales,
                     ~lm(total_volume_log10 ~ type + avg_price_shift1, data = .x)
                     )
         )

Instructor’s note: You might be tempted to mean-center and perhaps z-score the avg_price column, so that it’s in units of standard deviations of avocado prices. However, you need to be careful how your data is grouped/nested when you center or z-score data, to know whether you’re centering your data within grouping unit or across grouping units. That is, are you centering/scaling by a mean/SD calculated separately for each grouping unit, or for all observations across units? In some cases, centering/scaling data within unit can be dangerous, because you might accidentally be scaling away the effect of interest. Scaling data across units retains between-unit differences of interest, so it’s fine in this case. You would need to unnest your observation-level data to make sure your operations are not implicitly grouped before doing any mean-centering/SD-scaling operations.

Now that we’ve created the model column to hold each of our lm objects, let’s quickly see what the whole dataframe looks like:

avo_by_region

## # A tibble: 46 x 3
##    region              sales              model 
##    <chr>               <list>             <list>
##  1 Albany              <tibble [338 × 7]> <lm>  
##  2 Atlanta             <tibble [338 × 7]> <lm>  
##  3 BaltimoreWashington <tibble [338 × 7]> <lm>  
##  4 Boise               <tibble [338 × 7]> <lm>  
##  5 Boston              <tibble [338 × 7]> <lm>  
##  6 BuffaloRochester    <tibble [338 × 7]> <lm>  
##  7 Charlotte           <tibble [338 × 7]> <lm>  
##  8 Chicago             <tibble [338 × 7]> <lm>  
##  9 CincinnatiDayton    <tibble [338 × 7]> <lm>  
## 10 Columbus            <tibble [338 × 7]> <lm>  
## # … with 36 more rows

Aha! We see that in the list-column model, each of the top rows contains an lm object. Each of these lm objects individually can be operated on like any linear regression object stored in its own variable.

avo_by_region %>%
  pull(model) %>%
  pluck(1)

## 
## Call:
## lm(formula = total_volume_log10 ~ type + avg_price_shift1, data = .x)
## 
## Coefficients:
##      (Intercept)       typeorganic  avg_price_shift1  
##           5.0298           -1.5590           -0.2388

We can even %>% this further into summary() to generate the whole regression output table:

avo_by_region %>%
  pull(model) %>%
  pluck(1) %>%
  summary()

## 
## Call:
## lm(formula = total_volume_log10 ~ type + avg_price_shift1, data = .x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.42005 -0.10962  0.00827  0.10662  0.51076 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       5.02979    0.02068 243.168  < 2e-16 ***
## typeorganic      -1.55898    0.02687 -58.023  < 2e-16 ***
## avg_price_shift1 -0.23883    0.04608  -5.183 3.78e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1693 on 335 degrees of freedom
## Multiple R-squared:  0.9605, Adjusted R-squared:  0.9603 
## F-statistic:  4077 on 2 and 335 DF,  p-value: < 2.2e-16

While this is not in spirit a multiple regression tutorial, we are doing a multiple regression, so now’s a handy time to review interpreting multiple regression coefficients:

• Intercept: This tells us the estimated log-10-transformed volume of conventional avocados sold in a hypothetical week where conventional avocados cost $1.00 each
  • By default the type is dummy-coded so the baseline level is the first alphabetical level, and C comes before O, so the intercept assumes the baseline level is conventional avocados. Fine by me!
  • This is why we created avg_price_shift1 earlier. If we hadn’t shifted the 0 point of avg_price, the intercept would be estimated for a week where conventional avocados cost $0.00 each. I’m assuming avocados are never free (otherwise we millenials would be able to afford mortgages), so shifting the baseline can make the regression coefficient estimates easier to interpret.
• typeorganic: This tells us the estimated log-10-transformed difference in volume of avocados sold, between conventional and organic avocados, in a hypothetical week where both avocado types cost $1.00 each
• avg_price_shift1: This tells us the expected difference in volume of avocados sold between a week where avocados cost $1.00 each and a week where avocados are one price unit, or one dollar, more expensive ($2.00 each)

Importantly, because this is a multiple regression, the typeorganic term tells us the estimate difference in volume of avocados sold with price “held constant”, or if both conventional and organic avocados cost the same in a hypothetical week. If we can estimate the difference volume sold when the two avocado types cost the same, and people are STILL buying more conventional avocados, there must be an extenuating reason why people don’t buy organic avocados even if they’re not more expensive.

The console output that you get from calling summary() on a single lm object is great for visual inspection, but there’s a whole lot of output here that is not the most efficient to parse programmatically with code. As such, we won’t be using map() inside mutate() to call summary() on every single model object in the model column. There are other functions we can use to give us cleaner, more ready-to-use output, that we’ll explore next.

3.7.2 Extracting statistics from models

In order to programmatically extract coefficients, standard errors, and other statistics of interest from each model object in the model column of avo_by_region, our nested avocado sale data, we are going to use the tidy() helper function in the broom R package. broom is designed for generating tidyverse-safe versions of statistical analysis objects to get you from models to plots & interpretation with as little hacking as possible.

First, to demonstrate the output of broom::tidy(), we can call it just on the first model element of the model column of avo_by_region just as we called summary() on it above.

avo_by_region %>%
  pull(model) %>%
  pluck(1) %>%
  # The only thing that's different is here,
  # the last command we call after we've plucked a single lm object:
  broom::tidy()

## # A tibble: 3 x 5
##   term             estimate std.error statistic   p.value
##   <chr>               <dbl>     <dbl>     <dbl>     <dbl>
## 1 (Intercept)         5.03     0.0207    243.   0        
## 2 typeorganic        -1.56     0.0269    -58.0  7.89e-177
## 3 avg_price_shift1   -0.239    0.0461     -5.18 3.78e-  7

Note that I am using the :: to call tidy() directly out of the broom package namespace without library()-ing the package. I prefer to use this syntax when I’m only using a single function from a package in a particular analysis document, to avoid loading too many function names into my session. (sometimes functions in different packages have overlapping names so you want to be careful!)

Calling tidy() on a single lm object gives us a tibble, with a column for each statistic of interest in the regression table, and a row for each term from the multiple regression. This has the same numbers as the core table output we got from calling summary() on the same object, but now with all the added benefits of a tibble. We can use select() to choose specific columns of statistics, filter() to choose specific rows corresponding to particular coefficients, mutate() to calculate new statistics from the ones we already have, and all that good stuff!

Another huge benefit of tidy() returning a tibble of model output for each model is that if we use tidy() inside map() inside mutate() to create a list-column of coefficients tibbles, we can create a list-column of tibble dataframes, which has specific tidyverse abilities.

avo_by_region %<>%
  mutate(coefs = map(model, ~broom::tidy(.x)))

The biggest tidyverse ability is that the resulting coefs list-column of coefficients dataframes can then be pulled out into long form with unnest().

# Storing this long dataframe into its own object
# so we can feed it into ggplot2 calls more quickly later!
avo_coefs <- avo_by_region %>%
  # See the note below for why we're using select()
  # to retain only these two columns
  select(region, coefs) %>%
  unnest(coefs)

avo_coefs

## # A tibble: 138 x 6
##    region              term             estimate std.error statistic   p.value
##    <chr>               <chr>               <dbl>     <dbl>     <dbl>     <dbl>
##  1 Albany              (Intercept)         5.03     0.0207    243.   0        
##  2 Albany              typeorganic        -1.56     0.0269    -58.0  7.89e-177
##  3 Albany              avg_price_shift1   -0.239    0.0461     -5.18 3.78e-  7
##  4 Atlanta             (Intercept)         5.72     0.0128    446.   0        
##  5 Atlanta             typeorganic        -1.57     0.0243    -64.8  1.67e-191
##  6 Atlanta             avg_price_shift1   -0.240    0.0305     -7.88 4.64e- 14
##  7 BaltimoreWashington (Intercept)         5.94     0.0182    327.   0        
##  8 BaltimoreWashington typeorganic        -1.51     0.0229    -65.9  6.99e-194
##  9 BaltimoreWashington avg_price_shift1   -0.165    0.0383     -4.31 2.14e-  5
## 10 Boise               (Intercept)         4.93     0.0107    460.   0        
## # … with 128 more rows

And here we are, a long dataframe with a row for each coefficient from the multiple regression, and a chunk of rows for the coefficients from the model fit to a particular metro area. Usefully, this long dataframe can then be fed into a ggplot2 plotting call, to visualize all the model outputs from each region together. That’s what we’ll get to next!

Notice that as of tidyr >=1.0.0, unnest()’s default behavior is to retain any other list-columns present in the df that aren’t getting unnested. If you unnested a list-column of dataframes that each had multiple rows, e.g. the coefs column, where each sub-dataframe has a row for each coefficient, the other list-column(s) getting retained would have their values repeated for each of the new rows created by unnesting the other column. Sometimes you might want that, but sometimes you might not want to increase your df’s memory footprint so much. In those cases, you should use select() to remove list-columns that you don’t want repeated before unnesting.

3.7.3 Visualizing many statistics

Now that we have our model coefficients in a long form dataframe, with a chunk of rows for each metropolitan area, let’s plot the model estimates for each region!

Here, we’ll focus on plotting the values of the typeorganic term, as this will index each region’s relative buying preference (or anti-preference?) for organic relative to conventional avocados. This estimate will conveniently also be adjusted for average avocado price across regions, so it should not be affected by general differences in avocado prices (and cost of living) from region to region.

avo_coefplot <- avo_coefs %>%
  filter(term == "typeorganic") %>%
  # Here, using fct_reorder() from the forcats pkg of the tidyverse
  # which provides a bunch of nice factor manipulating functions
  # super DUPER handy for reordering factor levels on a plot
  # this reorders the regions in order from highest to lowest estimate
  ggplot(aes(x = fct_reorder(region, desc(estimate)), y = estimate)) +
  # I like to put error bars before, aka UNDERNEATH, the points
  geom_errorbar(aes(ymin = estimate - std.error,
                    ymax = estimate + std.error),
                width = 0) +
  geom_point() +
  # let's annotate
  # to add guides to remind us how to interpret the log-10 scale
  # that stuff is brain-bendy!
  geom_hline(yintercept = c(-1, -2), linetype = 3) +
  annotate("text", x = 8, y = -1.05, hjust = 0,
           label = "-1 on log10 = 10% as many", color = "black") +
  annotate("text", x = 8, y = -1.95, hjust = 0,
           label = "-2 on log10 = 1% as many", color = "black") +
  # since the region names are long and would otherwise overlap on the plot,
  # tilting them to 45 degrees will make them all stick off the plot on an angle
  # the new helper function guide_axis() takes care of specifics
  # when called within scale_x_discrete() as below (since we have discrete labels)
  # (as of ggplot2 v3.3.0)
  scale_x_discrete(guide = guide_axis(angle = 45)) +
  labs(x = "Metro area",
       y = "Log-10 diff in volume of organic avos sold",
       title = "People really do not buy organic avocados",
       subtitle = "especially not people in Florida") +
  theme_bw()

avo_coefplot

On this plot, metro areas on the left of the plot are areas that buy a greater percentage of organic avocados relative to conventional, and metro areas on the right of the plot are areas that buy a lower percentage of avocados relative to conventional.

Be mindful of the log scale when reading these numbers off the plot–thinking on the log scale can be tricky! Since we’re using log-10 scaling here, a log-10 difference of -1 corresponds to a real difference of 10^-1, or a 10% reduction in number. Similarly, a log-10 difference of -2 corresponds to a real difference of 10^-2, or a 1% reduction in number. (A log-10 difference of -1.5 is not as intuitive, we’re all better off just calculating 10^1.5)

So in Seattle, we estimate that when the price premium for organic avocados is taken into account, people buy approximately 10% as many organic avocados as they do conventional avocados. Meanwhile, in Miami & Fort Lauderdale, we estimate that people buy less than 1% as many organic avocados as they do conventional avocados. A log-10 difference of 1 is a whole order of magnitude different in raw avocado units!

Instructor’s note: The following plot takes a bit of fiddling to set up, and requires chaining together many pieces of data wrangling logic. If learners have never made a model-estimate-plus-raw-data plot like this before, this may be a bit beyond scope to live-code. In that case, I recommend taking a couple minutes for learners to brainstorm possible strategies for extending the basic coefficient plot shown above, and then showing a fully rendered copy of the plot below as an example extension plot, qualitatively describing the added layers and annotations.

avo_raw_diffs <- avocado %>%
  # first, prepare to calculate the raw difference in weekly volume
  # of organic - conventional avocados bought
  # (must be that direction bc conventional is baseline in the regression)
  select(region, date, type, total_volume_log10) %>%
  pivot_wider(names_from = type,
              values_from = total_volume_log10) %>%
  mutate(diff_raw = organic - conventional) %>%
  # WestTextNewMexico has some weeks with no organic avocado sales
  # so we need to filter out some NA rows of diff_raw
  filter(!is.na(diff_raw)) %>%
  # for each region, calculate the mean and SE of the mean
  # of this difference
  group_by(region) %>%
  summarize(diff_mean = mean(diff_raw),
            diff_se = sd(diff_raw)/sqrt(length(diff_raw))) %>% 
  # in order for the plot to work
  # (for the regions to be ordered by the model estimate)
  # the model estimates need to be bound on to the raw diffs
  left_join(avo_coefs %>% filter(term == "typeorganic"),
            by = "region")

avo_coefplot +
  # error bars and points for the raw data
  # notice that we use the data argument to feed in different data
  geom_errorbar(aes(y = diff_mean,
                    ymin = diff_mean - diff_se,
                    ymax = diff_mean + diff_se),
                data = avo_raw_diffs,
                color = "olivedrab4",
                width = 0) +
  geom_point(aes(y = diff_mean),
             data = avo_raw_diffs,
             color = "olivedrab4") +
  # hand-adding a legend
  annotate("text", x = 45, y = -1.2, hjust = 1,
           label = "black = model-estimated", color = "black") +
  annotate("text", x = 45, y = -1.3, hjust = 1,
           label = "green = raw", color = "olivedrab4")

I think we’ve gotten pretty close to an answer (not THE answer, but AN answer) to our earlier research question: Do avocado buyers prefer buying conventional or organic avocados? By how much? How do these patterns differ from region to region?

Avocado buyers across the US appear to prefer buying conventional avocados over organic avocados; that is, they buy many more conventional avocados than organic avocados in a given week. This holds even when you adjust for the fact that organic avocados are usually more expensive. The magnitude of this buying preference against organic avocados varies across the country–the Pacific Northwest appears to avoid organic avocados much less than the rest of the country, and Florida appears to avoid organic avocados much more than the rest of the country.

That concludes the core analysis of this tutorial! If you like, continue following along for a rapid-fire demonstration of another use case of map() inside mutate() with nested dataframes to repeat analysis across many datasets containing similarly structured data.