Data fundamentals: Describe, wrangle, tidy

Drug treatment dataset

To elaborate further, we can use the example given in Wickham (2014), a drug treatment dataset in which two different treatments were administered to participants.

The data could be represented as:

An alternative organisation could be:

Both present the same information unambiguously – Table 3 is simply Table 2 transposed. However, neither is tidy as the observations are spread across both the rows and columns. This means that we need to apply different procedures to extract, perform computations on, and visually represent, these data.

Much better would be to organise the table into a tidy form. To do this we need to identify the variables:

1. person: a categorical nominal variable which takes three values: John Smith, Jane Doe, Mary Johnson.
2. treatment: a categorical nominal variable which takes values: a and b.
3. result: a measurement ratio (I think) variable with six recorded values (including the missing value): -, 16, 3, 2, 11,

Each observation is then a test result returned for each combination of person and treatment.

So, a tidy organisation for this dataset would be:

Gapminder population dataset

In chapter 12 of Wickham and Grolemund (2017), the benefits of this layout, particularly for working with R, are demonstrated using the gapminder dataset. I recommend reading this short chapter in full. We will be applying similar approaches in the technique part of this class (which follows shortly) and also the homework. To consolidate our conceptual understanding of tidy data let’s quickly look at the gapminder data, as it is a dataset we’re probably more likely to encounter.

First, a tidy version of the data:

So the variables:

1. country: a categorical nominal variable.
2. year: a date (cyclic ratio) variable.
3. cases: a ratio (count) variable.
4. population: a ratio (count) variable.

Each observation is therefore a recorded count of cases and population for a country in a year.

An alternative organisation of this dataset that appears in Wickham and Grolemund (2017) is below. This is untidy as the observations are spread across two rows. This makes operations that we might want to perform on the cases and population variables – for example computing exposure rates – somewhat tedious.

Imagine that the gapminder dataset instead reported values of cases separately by gender. A type of representation I’ve often seen in social science domains, probably as it is helpful for data entry, is where observations are spread across the columns. This too creates problems for performing aggregate functions, but also for specifying visualization designs (in ggplot2) as we will discover in the next session.

Transform with dplyr

dplyr is one of the most important packages for supporting modern data analysis workflows. The package provides a grammar of data manipulation, with access to functions that can be variously combined to support most data processing and transformation activity. Once you become familiar with dplyr functions (or verbs) you will find yourself generating analysis templates to re-use whenever you work on a new dataset.

All dplyr functions work in the same way:

1. Start with a data frame.
2. Pass some arguments to the function which control what you do to the data frame.
3. Return the updated data frame.

So every dplyr function expects a data frame and will always return a data frame.

Use pipes %>% with dplyr

dplyr is most effective when its functions are chained together – you will see this shortly as we explore the New York bikeshare data. This chaining of functions can be achieved using the pipe operator (%>%). Pipes are used for passing information in a program. They take the output of a set of code (a dplyr specification) and make it the input of the next set (another dplyr specification).

Pipes can be easily applied to dplyr functions, and the functions of all packages that form the tidyverse. I mentioned in Session 1 that ggplot2 provides a framework for specifying a layered grammar of graphics (more on this in Session 3). Together with the pipe operator (%>%), dplyr supports a layered grammar of data manipulation.

count() rows

This might sound a little abstract so let’s use and combine some dplyr functions to generate some statistical summaries on the New York bikeshare data.

First we’ll count the number of trips made in Jun 2020 by gender. dplyr has a convenience function for counting, so we could run the code below, also in the 02-template.Rmd for this session. I’ve commented the code block to convey what each line achieves.

ny_trips %>%  # Take the ny_trips data frame.
  count(gender, sort=TRUE) # Run the count function over the data frame and set the sort parameter to TRUE.
##    gender       n
## 1    male 1044621
## 2  female  586361
## 3 unknown  251291

There are a few things happening in the count() function. It takes the gender variable from ny_trips, organises or groups the rows in the data frame according to its values (female | male | unknown), counts the rows and then orders the summarised output descending on the counts.

summarise() over rows

Often you will want to do more than simply counting and you may also want to be more explicit in the way the data frame is grouped for computation. We’ll demonstrate this here with a more involved analysis of the usage data and using some key aggregate functions (Table 10).

A common workflow is to combine group_by() and summarise(), and in this case arrange() to replicate the count() example.

ny_trips %>% # Take the ny_trips data frame.
  group_by(gender) %>% # Group by gender.
  summarise(count=n()) %>% # Count the number of observations per group.
  arrange(desc(count)) # Arrange the grouped and summarised (collapsed) rows according to count.
## # A tibble: 3 x 2
##  gender    count
##  <chr>     <int>
## 1 male    1044621
## 2 female   586361
## 3 unknown  251291

In ny_trips there is a variable measuring trip duration in seconds (trip_duration) and distinguishing casual users from those formally registered to use the scheme (user_type - Customer vs. Subscriber). It may be instructive to calculate some summary statistics to see how trip duration varies between these groups.

The code below uses group_by(), summarise() and arrange() in exactly the same way, but with the addition of other aggregate functions profiles the trip_duration variable according to central tendency and by user_type.

ny_trips %>% # Take the ny_trips data frame.
  group_by(user_type) %>% # Group by user type.
  summarise( # Summarise over the grouped rows, generate a new variable for each type of summary.
    count=n(),
    avg_duration=mean(trip_duration/60),
    median_duration=median(trip_duration/60),
    sd_duration=sd(trip_duration/60),
    min_duration=min(trip_duration/60),
    max_duration=max(trip_duration/60)
    ) %>%
  arrange(desc(count)) # Arrange on the count variable.

## # A tibble: 2 x 6
##  user_type    count avg_duration median_duration sd_duration min_duration max_duration
##  <chr>        <int>        <dbl>           <dbl>       <dbl>        <dbl>        <dbl>
## 1 Subscriber 1306688         20.2            14.4        110.         1.02         33090
## 2 Customer    575585         43.6            23.2        393.         1.02         46982

Clearly there are some outlier trips that may need to be examined. Bikeshare schemes are built to incentivise short journeys of <30 minutes, but the maximum trip duration recorded above is clearly erroneous – 32 days. Ignoring these sorts of outliers by calculating the trip durations at the 95th percentiles is instructive. The max trip duration at the 95th percentile for Subscribers was almost 27 minutes and for Customers was 1 hours 26 mins. It makes sense that more casual users may have longer trip durations, as they are more likely to be tourists or occasional cyclists using the scheme for non-utility trips.

Returning to the breakdown of usage by gender, an interesting question is whether or not the male-female split in bikehare is similar to that of the cycling population of New York City as a whole. This might tell us something about whether the bikeshare scheme could be representative of wider cycling. This could be achieved with the code below. A couple of new additions: we use filter(), to remove observations where the gender of the cyclist is unknown. We also use mutate() for the first time, which allows us to modify or create new variables.

ny_trips %>% # Take the ny_trips data frame.
  filter(gender != "unknown") %>% # Filter out rows with the value "unknown" on gender.
  group_by(gender) %>% # Group by gender.
  summarise(count=n()) %>% # Count the number of observations per group.
  mutate(prop=count/sum(count)) %>% # Add a new column called `prop`, divide the value in the row for the variable count by the sum of the count variable across all rows.
  arrange(desc(count)) # Arrange on the count variable.

## # A tibble: 2 x 3
##  gender   count  prop
##  <chr>    <int> <dbl>
## 1 male   1044621 0.640
## 2 female  586361 0.360

As I’ve commented each line you hopefully get a sense of what is happening in the code above. I mentioned that dplyr functions read like verbs. This is a very deliberate design decision. With the code laid out as above – each dplyr verb occupying a single line, separated by a pipe (%>%) – you can generally understand the code with a cursory glance. There are obvious benefits to this. Once you are familiar with dplyr it becomes very easy to read, write and share code.

Manipulate dates with lubridate

Let’s continue this investigation of usage by gender, and whether bikeshare might be representative of regular cycling, by profiling how usage varies over time. To do this we will need to work with date-time variables. The lubridate package provides various convenience functions for this.

In the code block below we extract the day of week and hour of day from the start_time variable using lubridate’s day accessor functions. Documentation on these can be accessed in the usual way (?<function-name>), but reading down the code it should be clear to you how this works. Next we count the number of trips made by hour of day, day of week and gender. The summarised data frame will be re-used several times in our analysis, so we store it as an object with a suitable name (ny_temporal) using the assignment operator.

# Create a hod dow summary by gender and assign it the name "ny_temporal".
ny_temporal <- ny_trips %>%  # Take the ny_trips data frame.
  mutate(
    day=wday(start_time, label=TRUE), # Create a new column identify dow.
    hour=hour(start_time)) %>% # Create a new column identify hod.
  group_by(gender, day, hour) %>% # Group by day, hour, gender.
  summarise(count=n()) %>% # Count the grouped rows.
  ungroup()

In Figure 2 below these derived data are plotted. The template contains ggplot2 code for creating the graphic. Don’t obsess too much on it – more on this next session. The plot demonstrates a familiar weekday-weekend pattern of usage. Trip frequencies peak in the morning and evening rush hours during weekdays and mid/late-morning and afternoon during weekends. This is consistent with typical travel behaviour. Notice though that the weekday afternoon peak is much larger than the morning peak. There are several speculative explanations for this and re-running the plot on Subscriber users only may be instructive. A secondary observation is that whilst men and women share this overall pattern of usage, the relative number of trips taken each day of week varies. Men make many more trips at peak times during the start of the week than they do later in the week. The same pattern does not appear for women. This is certainly something to follow up on, for example by collecting data over a longer period of time.

Figure 2: Line charts generated with ggplot2. Plot data computed using dplyr and lubridate.

Relate tables with join()

Trip distance is not recorded directly in the ny_trips table, but may be important for profiling usage behaviour. Calculating trip distance is eminently achievable as the ny_trips table contains the origin and destination station of every trip and the ny_stations table contains coordinates corresponding to those stations. To relate the two tables, we need to specify a join between them.

A sensible approach is to:

1. Select all uniquely cycled trip pairs (origin-destination pairs) that appear in the ny_trips table.
2. Bring in the corresponding coordinate pairs representing the origin and destination stations by joining on the ny_stations table.
3. Calculate the distance between the coordinate pairs representing the origin and destination.

The code below is one way of achieving this.

od_pairs <- ny_trips %>% # Take the ny_trips data frame.
select(start_station_id, end_station_id) %>% unique() %>% # Select trip origin and destination (OD) station columns and extract unique OD pairs.
  left_join(ny_stations %>% select(stn_id, longitude, latitude), by=c("start_station_id"="stn_id")) %>% # Select lat, lon columns from ny_stations and join on the origin column.
  rename(o_lon=longitude, o_lat=latitude) %>% # Rename new lat, lon columns -- associate with origin station.
  left_join(ny_stations %>% select(stn_id, longitude, latitude), by=c("end_station_id"="stn_id")) %>% # Select lat, lon columns from ny_stations and join on the destination column.
  rename(d_lon=longitude, d_lat=latitude) %>%  # Rename new lat, lon columns -- associate with destination station.
  rowwise() %>% # For computing distance calculation one row-at-a-time.
  mutate(dist=geosphere::distHaversine(c(o_lat, o_lon), c(d_lat, d_lon))/1000) %>% # Calculate distance and express in kms.
  ungroup()

The code block above introduces some new functions: select() to pick or drop variables, rename() to rename variables and a convenience function for calculating straight line distance from polar coordinates (distHaversine()). The key function to emphasise is the left_join(). If you’ve worked with relational databases and SQL, dplyr’s join functions will be familiar. In a left_join, all the values from the main table are retained, the one on the left – ny_trips, and variables from the table on the right (ny_stations) are added. We specify explicitly the variable on which the tables should be joined with the by= parameter, station_id in this case. If there is a station_id in ny_trips that doesn’t exist in ny_stations then NA is returned.

Other join functions provided by dplyr are in the table below. Rather than discussing each, I recommend consulting Chapter 13 of Wickham and Grolemund (2017).

Figure 3: Histograms generated with ggplot2. Plot data computed using dplyr and lubridate

From the newly created distance variable, we can calculate the average (mean) trip distance for the 1.9m trips – 1.6km. This might seem very short, but remember that the distance calculation is problematic in that these are straight-line distances between pairs of docking stations. Really we should be calculating network distances derived from the cycle network in New York. A separate reason – discovered when generating a histogram on the dist variable – is that there are a large number of trips (124,403) that start and end at the same docking station. Initially these might seem to be unsuccessful hires – people failing to undock a bike for example. We could investigate this further by paying attention to the docking stations at which same origin-destination trips occur, as in the code block below.

ny_trips %>%
  filter(start_station_id==end_station_id) %>%
  group_by(start_station_id) %>% summarise(count=n()) %>%
  left_join(ny_stations %>%  select(stn_id, name), by=c("start_station_id"="stn_id")) %>%
  arrange(desc(count))

## # A tibble: 958 x 3
##    start_station_id count name
##    <chr>            <int> <chr>
##  1 ny3423            2017 West Drive & Prospect Park West
##  2 ny3881            1263 12 Ave & W 125 St
##  3 ny514             1024 12 Ave & W 40 St
##  4 ny3349             978 Grand Army Plaza & Plaza St West
##  5 ny3992             964 W 169 St & Fort Washington Ave
##  6 ny3374             860 Central Park North & Adam Clayton Powell Blvd
##  7 ny3782             837 Brooklyn Bridge Park - Pier 2
##  8 ny3599             829 Franklin Ave & Empire Blvd
##  9 ny3521             793 Lenox Ave & W 111 St
## 10 ny2006             782 Central Park S & 6 Ave
## # … with 948 more rows

All of the top 10 docking stations are either in parks, near parks or located along river promenades. This coupled with the fact that these trips occur in much greater relative number for casual than regular users (Customer vs Subscriber) is further evidence that these are valid trips.

Write functions of your own

Through most of the course we will be making use of functions written by others – mainly those developed for packages that form the tidyverse and therefore that follow a consistent syntax. However, there may be times where you need to abstract over some of your code to make functions of your own. Chapter 19 of Wickham and Grolemund (2017) presents some helpful guidelines around the circumstances under which the data scientist typically tends to write functions. Most often this is when you find yourself copy-pasting the same chunks of code with minimal adaptation.

Functions have three key characteristics:

1. They are (usually) named – the name should be expressive and communicate what the function does (we talk about dplyr verbs).
2. They have brackets <function()> usually containing arguments – inputs which determine what the function does and returns.
3. Immediately followed by <function()> are {} used to contain the body – in this is code that performs a distinct task, described by the function’s name.

Effective functions are short, perform single discrete operations and are intuitive.

You will recall that in the ny_trips table there is a variable called birth_year. From this we can derive cyclists’ approximate age. Below I have written a function get_age() for doing this. The function expects two arguments: yob – a year of birth as type chr; yref – a reference year. In the body, lubridate’s as.period function is used to calculate the time in years that elapsed, the value that the function returns. Once defined, and loaded into the session by being executed, it can be used (as below).

# Function for calculating time elapsed between two dates in years (age).
get_age <- function(yob, yref) {
    period <- lubridate::as.period(lubridate::interval(yob, yref),unit = "year")
    return(period$year)
}

ny_trips <- ny_trips %>% # Take the ny_trips data frame.
  mutate(
    age=get_age(as.POSIXct(birth_year, format="%Y"), as.POSIXct("2020", format="%Y")) # Calculate age from birth_date.
    )

We can use the two new derived variables – distance travelled and age – in our analysis. In Figure 4, we explore how approximate travel speeds vary by age, gender and trip distance. The code used to generate the summary data and plot is in the template file. Again the average “speed” calculation should be treated very cautiously as it is based on straight line distances and it is very difficult to select out “utility” from “leisure” trips. I have tried to do the latter by selecting trips that occur only on weekdays and that are made by Subscriber cyclists. Additionally, due to the heavy subsetting data become a little volatile for certain age groups and so I’ve aggregated the age variable into 5-year bands. Collecting more data is probably a good idea.

There are nevertheless some interesting patterns. Men tend to cycle at faster speeds than do women, although this gap narrows with the older age groups. The effect of age on speed cycled is more apparent for the longer trips. This trend is reasonably strong, although the volatility in the older age groups for trips >4.5km suggests we probably need more data and a more involved analysis to establish this. For example, it may be that the comparatively rare occurrence of trips in the 65-70 age group is made by only a small subset of cyclists. A larger dataset would result in a regression to the mean effect and negate any noise caused by outlier individuals. Certainly Figure 1 is an interesting data graphic – and the type of exploratory analysis demonstrated here, using dplyr functions, is most definitely consistent with that identified in the previous session when introducing Social Data Science.

Figure 4: Line charts generated with ggplot2. Plot data computed using dplyr and lubridate