Module road map

Overview

Artwork by @allison_horst

To wrangle data is to clean and prepare data for analysis. Our goal is to get the data into a tidy shape that can be readily plotted and analyzed. I like to think about data analysis in terms of painting a room. In painting, an excellent result comes from careful and thorough preparation. About 90% of painting is preparation – cleaning the walls, taping off trim, protecting the floor, getting the right tools, etc. Data analysis is similar – and data wrangling is the preparation. For most projects, most of your time will be spent preparing the data for analysis. If you skip this step or do a shoddy job – your data analysis will be poor.

The dplyr functions to wrangle data

dplyr is a R package that is part of the tidyverse — it’s the workhorse for wrangling your data.

In this Module, we will explore the following wrangling functions:

1. Choose certain cases (rows) in your dataframe with filter()

2. Choose certain variables (columns) in your dataframe with select()

3. Collapse many values down to summaries with summarize()

4. Create or modify variables with mutate()

5. Give variables a new name with rename()

6. Order your cases (rows) with arrange()

7. Merge two dataframes together with join()

Introduction to the data

In this Module we’re going to work with a very large data repository compiled by The Washington Post (WAPO) that tracked the path of every opioid pain pill, from manufacturer to pharmacy, in the United States between 2006 and 2014. The data scientists at WAPO created a package (called arcos) that makes the data publicly available.

The opioid epidemic in America has been devastating for individuals and communities. To provide context to the data exploration in this Module, please watch the following video from Last Week Tonight with John Oliver. Please note that the video contains some lewd language, you can skip the video without it affecting your grade if you like. You can also peruse this National Public Radio article and this Nature paper that also provide context.

We will work with the county level data provided by WAPO in this Module. I already downloaded it from the arcos package and saved it as a dataframe in our RStudio Cloud project. Let’s import the the dataframe.

opioid_counties <- here("data", "opioid_counties.Rds") %>% 
  read_rds()

Each county in the US has up to nine rows of data in this dataframe, one for each year from 2006 to 2014. The dataframe includes the following variables:

• fips is a 5 digit code that identifies each county

• county

• state is the state abbreviation

• year

• number_pills is the number of pills purchased – specifically the number of Oxycodone & Hydrocodone pills purchased by pharmacies in the county

• population is the total number of people living in the county

Let’s take a glimpse at the data.

opioid_counties %>% 
  glimpse()

## Rows: 27,758
## Columns: 6
## $ fips         <chr> "45001", "45001", "45001", "45001", "45001", "45001", "45…
## $ county       <chr> "ABBEVILLE", "ABBEVILLE", "ABBEVILLE", "ABBEVILLE", "ABBE…
## $ state        <chr> "SC", "SC", "SC", "SC", "SC", "SC", "SC", "SC", "SC", "LA…
## $ year         <int> 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 200…
## $ number_pills <dbl> 363620, 402940, 424590, 467230, 539280, 566560, 589010, 5…
## $ population   <int> 25821, 25745, 25699, 25347, 25643, 25515, 25387, 25233, 2…

Let’s also view the first set of rows. The head() function will display the first set of rows (we specify to see 18 rows in this case by including n = 18). If you don’t specify a number of rows, then the default number of rows (6) will be printed.

opioid_counties %>% 
  head(n = 18)

The filter() function

The filter() dplyr verb is used to choose specific rows of data that meet certain conditions. R provides the standard set of selection operators:

• >

• >=

• <

• <=

• != (not equal to)

• == (equal to)

Boolean operators are also permitted:

• & is “and”

• | is “or”

• ! is “not”

opioid_counties2014 <- opioid_counties %>% 
  filter(year == 2014)

opioid_counties2014 %>% 
  head()

opioid_counties2014_CO <- opioid_counties %>% 
  filter(state == "CO" & year == 2014)

opioid_counties2014_CO %>% 
  head()

The select() function

The select() dplyr verb is used to select only certain columns (i.e., variables) in your dataframe. Let’s choose just state, county, and number_pills.

opioid_counties_subset <- opioid_counties %>% 
  select(state, county, number_pills)

opioid_counties_subset %>% 
  head()

You can also use select() to reorder variables. For example, if you simply wanted to move state, county, and number_pills to the front of your dataframe, but still keep the other variables, you could modify the select() function above as follows:

select(state, county, number_pills, everything())

The everything() function puts all the other variables after the three specifically mentioned.

The summarize function

The summarize() dplyr verb is used to summarize your data – for example, to create summary statistics.

opioid_counties %>% 
  summarize(mean_number_pills = mean(number_pills, na.rm = TRUE), 
            sd_number_pills = sd(number_pills, na.rm = TRUE))

opioid_counties %>% 
  group_by(year) %>% 
  summarize(mean_number_pills = mean(number_pills, na.rm = TRUE)) 

sum_pills_year <- opioid_counties %>% 
  group_by(year) %>% 
  summarize(total_pills = sum(number_pills, na.rm = TRUE))

sum_pills_year

sum_pills_year %>% 
  ggplot(aes(x = year)) +
  geom_line(aes(y = total_pills)) +
  theme_bw() +
  labs(title = "Change in total number of opioid pills purchased in the US from 2006 to 2014?",
       y = "Number of opioid pills puchased",
       x = "Year") 

The mutate() function

The mutate() dplyr verb is used to create new variables.

opioid_counties <- opioid_counties %>% 
  mutate(pills_person = number_pills/population) 

opioid_counties %>% 
  head()

The rename() function

The rename() dplyr verb is used to rename a variable. Let’s imagine that we want to change the name of number_pills to opioid_pills.

opioid_counties_rename <- opioid_counties %>% 
  rename(opioid_pills = number_pills)

opioid_counties_rename %>% 
  head()

The arrange() function

The arrange() dplyr verb will sort a dataframe by a variable or set of variables. Let’s filter to choose only CO counties in 2010 and sort by number of pills purchased per person.

opioid_counties %>% 
  filter(year == 2010 & state == "CO") %>% 
  arrange(pills_person)

If instead you wanted to sort from high to low, then use the desc() function:

opioid_counties %>% 
  filter(year == 2010 & state == "CO") %>% 
  arrange(desc(pills_person))

We can also create a chart of the same data – a bar chart with county on the y-axis and pill_person on the x-axis will do the trick. If we’d like to order the counties by the pills per capita, then we need to use the fct_reorder function (see code below). Rather than y = county, we need y = fct_reorder(county, pills_person). fct_reorder in this instance needs two arguments, the variable that you want to put on the y-axis (county), and the ordering that you desire (we want to order by pills_person).

opioid_counties %>% 
  filter(year == 2010 & state == "CO") %>% 
ggplot(aes(x = pills_person, y = fct_reorder(county, pills_person))) +
  geom_col(fill = "#E24E42") +
  labs(title = "Number of pills purchased per person in 2010 by county in CO",
       x = "Number of pills per person",
       y = "County",
       caption = "Source: The Washington Post, ARCOS") +
  theme_minimal()

The join functions

The join dplyr verbs are used to join/merge two dataframes together. Let’s read in a related dataframe, called adi_counties.Rds, that provides the Area Disadvantage Index (ADI) for each county. The ADI is based on a measure created by the Health Resources & Services Administration (HRSA). It is based on US Census records, and serves as a measure of concentrated socio-economic disadvantage in the county. In the dataframe, there is one row of data per county in the US, and then a series of variables that describe the county. In this example, we’ll use the following variables:

• fips is the county ID.

• ADI is the area disadvantage index – a higher score means more disadvantage

• WhiteNonHispanicPct2010 is the percent of county residents who identify as non-Hispanic White.

For this example, we will dichotomize WhiteNonHispanicPct2010, naming it mincounty. This binary indicator will be used to compare counties where fewer than 50% of the residents identify as non-Hispanic White to all others. We’ll do this using the case_when() function in dplyr. Each case_when() statement has two parts – the part before the tilde (~) and the part after the tilde. The part before the tilde is called the left hand side (LHS) and the part after the tilde is called the right hand side (RHS). case_when() uses the LHS part of the statement as the condition to test for and the RHS as the value to return. case_when() works in sequence through the conditions, finding when the LHS is true for each case (i.e., row of data), and then returning the value associated with that condition.

In the example below, mincounty is a new variable that is created using case_when(). Each row in the dataframe will be evaluated, if WhiteNonHispanicPct2010 < 50, then “minority county” will be recorded for mincounty. The second case_when() statement is only evaluated if the first case_when() statement is false. In the second case_when() statement, if WhiteNonHispanicPct2010 is greater than or equal to 50, then “white county” will be recorded for mincounty. Any rows that don’t match any of the conditions will be assigned NA (NA is R’s missing data indicator) for mincounty.

adi <- here("data", "adi_counties.Rds") %>% 
  read_rds() %>% 
  mutate(mincounty = case_when(WhiteNonHispanicPct2010 < 50 ~ "minority county", 
                               WhiteNonHispanicPct2010 >= 50 ~ "white county")) %>%
  select(fips, ADI, WhiteNonHispanicPct2010, mincounty)

adi %>% 
  head()

Now, let’s merge the adi dataframe with the 2010 opioid data, since fips is in both dataframes, and it uniquely identifies each county, we can join the dataframes together using fips as the key.

opioid_counties2010 <-
  opioid_counties %>% 
  filter(year == 2010) %>% 
  left_join(adi, by = "fips")

opioid_counties2010 %>% 
  head()

Advanced tip: There are many different types of joins available, and the best choice depends on the qualities of the dataframes to be joined and the analyst’s needs. You can explore the types here.

With these data compiled, let’s take the opportunity to determine if counties with more socio-economic disadvantage had more pills coming into their communities. This is an important public health question. The socio-economic vitality of a community plays an important role in promoting health, well-being, and prosperity for residents. On the flip side, communities characterized by socio-economic strain pose threats to the health, well-being, and prosperity of residents. As such, public health scientists recognize socio-economic advantage vs. disadvantage, of both an individual and his/her community, as a key social determinant of health. HealthyPeople.gov defines social determinants of health and describes their importance:

Social determinants of health are conditions in the environments in which people are born, live, learn, work, play, worship, and age that affect a wide range of health, functioning, and quality-of-life outcomes and risks. Conditions (e.g., social, economic, and physical) in these various environments and settings (e.g., school, church, workplace, and neighborhood) have been referred to as “place.” In addition to the more material attributes of “place,” the patterns of social engagement and sense of security and well-being are also affected by where people live. Resources that enhance quality of life can have a significant influence on population health outcomes. Examples of these resources include safe and affordable housing, access to education, public safety, availability of healthy foods, local emergency/health services, and environments free of life-threatening toxins. Understanding the relationship between how population groups experience “place” and the impact of “place” on health is fundamental to the social determinants of health—including both social and physical determinants.

You probably know from the news cycle over the past decade that the opioid epidemic hit poor White communities with few resources especially hard. If you’d like to read more about the opioid epidemic and the role of neighborhood disadvantage, particularly in White communities, please read this interesting report by the Brookings Institute entitled “How can Policy Address the Opioid Crisis and Despair in America”.

Let’s use our data to take a look at the relationship between ADI and number of pills per person purchased. Given our knowledge that the opioid pandemic was especially wide spread in predominantly non-Hispanic White communities, we will facet our scatterplot by the mincounty indicator created earlier.

opioid_counties2010 %>% 
  select(ADI, pills_person, mincounty) %>% 
  drop_na() %>% 
  ggplot(mapping = aes(x = ADI, y = pills_person)) +
  geom_jitter(color = "dark gray") + 
  geom_smooth(method = "loess", se = FALSE, formula = y~x) +
  facet_wrap(~mincounty) +
  theme_bw() +
  labs(title = "Is county-level socio-economic disadvantage associated with more opioid pills per capita?",
       x = "Area Disadvantage Index",
       y = "Number of pills purchased per person")

Importing and exporting data in R

Before we close out this Module, let’s practice importing (reading in) and exporting (writing out) data files.

adi <- here("data", "adi_counties.Rds") %>% 
  read_rds()

adi <- here("data", "adi_counties.csv") %>% 
  read_csv()

opioid_counties2010 %>% 
  saveRDS(here("data", "opioid_counties2010.Rds"))

opioid_counties2010 %>% 
  write_csv(here("data", "opioid_counties2010.csv"))

Resources

Go further, learn more:

• dplyr tidyverse site (including cheat sheets)

• R for Data Science