OpenCaseStudies - Health Expenditure
Motivation
Health policy in the Staes is complicated, and several forms of healthcare coverage existed in the United States of America, including both federal goverment-led healthcare policy, and private insurance company. Before making any inference about the relationship between health condition and health policy, it is important for us to have a general idea about healthcare economics in the States. Thus, We are interested in getting sense of the health expenditure, including healthcare coverage and healthcare spending, across States. More specifically, the quesiton are
- Is there a relationship between healthcare coverage and healthcare spending in the United States?
- How does the spending distribution change across geographic regions in the United States?
- Does the relationship between healthcare coverage and healthcare spending in the United States change from 2013 to 2014?
In this case study, we’ll walk you through collecting data, importing data, cleaning data, wrangling data, and visulizing the data, using well-established and commonly used packages, including library(datasets), library(tidyr), library(dplyr), ggplot2, and ggrepel.
What is the data?
Image source from US Department of Health and Human Services
Healthcare data
We will be using the data from the Henry J Kaiser Family Foundation (KFF).
- Health Insurance Coverage of the Total Population - Includes years 2013-2016
- Health Care Expenditures by State of Residence (in millions) - Includes years 1991-2014
We have downloaded, re-named and saved these files in the GitHub repository under the data/KFF/ directory.
Now, before we dig into the data analysis, we need to introduce a set of R packages that we will use to analyze the data.
Data Import
Introduction to “Tidy data”
The tidyverse is “an opinionated collection of R packages designed for data science. All packages share an underlying philosophy and common APIs.”
Another way of putting it is that it’s a set of packages that are useful specifically for data manipulation, exploration and visualization with a common philosphy.
What is this common philosphy?
The common philosphy is called “tidy” data. It is a standard way of mapping the meaning of a dataset to its structure.
In tidy data:
- Each variable forms a column.
- Each observation forms a row.
- Each type of observational unit forms a table.
Below, we are interested in transformating the table on the right to the the table on the left, which is considered “tidy”.
Working with tidy data is useful because it creates a structured way of organizing data values within a data set. This makes the data analysis process more efficient and simplifies the development of data analysis tools that work together. In this way, you can focus on the problem you are investigating, rather than the uninteresting logistics of data.
1. What is in the tidyverse?
We can install and load the set of R packages using install.packages("tidyverse") function.
When we load the tidyverse package using library(tidyverse), there are six core R packages that load:
- readr, for data import.
- tidyr, for data tidying.
- dplyr, for data wrangling.
- ggplot2, for data visualisation.
- purrr, for functional programming.
- tibble, for tibbles, a modern re-imagining of data frames.
Here, we load in the tidyverse.
library(tidyverse)These packages are highlighted in bold here:
Because these packages all share the “tidy” philosphy, the data analysis workflow is easier as you move from package to package.
Here, we will focus on the readr, tidyr and dplyr R packages to import data, to transform data to the “tidy” format, and to wrangle data.
Next, we will give a brief description of the features in each of these packages.
There are several base R functions that allow you read in data into R, which you may be familiar with such as read.table(), read.csv(), and read.delim(). Instead of using these, we will use the functions in the readr R package. The main reasons for this are
- Compared to equivalent base R functions, the functions in
readrare around 10x faser. - You can specify the column types (e.g character, integer, double, logical, date, time, etc)
- All parsing problems are recordered in a data frame.
Read data using the readr R package
library(readr)The main functions in readr are:
readr functions |
Description |
|---|---|
read_delim() |
reads in a flat file data with a given character to separate fields |
read_csv() |
reads in a CSV file |
read_tsv() |
reads in a file with values separated by tabs |
read_lines() |
reads only a certain number of lines from the file |
read_file() |
reads a complete file into a string |
write_csv() |
writes data frame to CSV |
A useful cheatsheet for the functions in the readr package can be found on RStudio’s website:
1. Read in data
Read in health healthcare coverage data
Let’s try reading in some data. We will begin by reading in the healthcare-coverage.csv data.
If we want to see what the header of the file looks like, we can use the read_lines() function to peak at the first few lines.
read_lines(file = "./data/KFF/healthcare-coverage.csv", n_max = 10) [1] "\"Title: Health Insurance Coverage of the Total Population | The Henry J. Kaiser Family Foundation\""
[2] "\"Timeframe: 2013 - 2016\""
[3] "\"Location\",\"2013__Employer\",\"2013__Non-Group\",\"2013__Medicaid\",\"2013__Medicare\",\"2013__Other Public\",\"2013__Uninsured\",\"2013__Total\",\"2014__Employer\",\"2014__Non-Group\",\"2014__Medicaid\",\"2014__Medicare\",\"2014__Other Public\",\"2014__Uninsured\",\"2014__Total\",\"2015__Employer\",\"2015__Non-Group\",\"2015__Medicaid\",\"2015__Medicare\",\"2015__Other Public\",\"2015__Uninsured\",\"2015__Total\",\"2016__Employer\",\"2016__Non-Group\",\"2016__Medicaid\",\"2016__Medicare\",\"2016__Other Public\",\"2016__Uninsured\",\"2016__Total\""
[4] "\"United States\",\"155696900\",\"13816000\",\"54919100\",\"40876300\",\"6295400\",\"41795100\",\"313401200\",\"154347500\",\"19313000\",\"61650400\",\"41896500\",\"5985000\",\"32967500\",\"316159900\",\"155965800\",\"21816500\",\"62384500\",\"43308400\",\"6422300\",\"28965900\",\"318868500\",\"157381500\",\"21884400\",\"62303400\",\"44550200\",\"6192200\",\"28051900\",\"320372000\""
[5] "\"Alabama\",\"2126500\",\"174200\",\"869700\",\"783000\",\"85600\",\"724800\",\"4763900\",\"2202800\",\"288900\",\"891900\",\"718400\",\"143900\",\"522200\",\"4768000\",\"2218000\",\"291500\",\"911400\",\"719100\",\"174600\",\"519400\",\"4833900\",\"2263800\",\"262400\",\"997000\",\"761200\",\"128800\",\"420800\",\"4834100\""
[6] "\"Alaska\",\"364900\",\"24000\",\"95000\",\"55200\",\"60600\",\"102200\",\"702000\",\"345300\",\"26800\",\"130100\",\"55300\",\"37300\",\"100800\",\"695700\",\"355700\",\"22300\",\"128100\",\"60900\",\"47700\",\"90500\",\"705300\",\"324400\",\"20300\",\"145400\",\"68200\",\"55600\",\"96900\",\"710800\""
[7] "\"Arizona\",\"2883800\",\"170800\",\"1346100\",\"842000\",\"N/A\",\"1223000\",\"6603100\",\"2835200\",\"333500\",\"1639400\",\"911100\",\"N/A\",\"827100\",\"6657200\",\"2766500\",\"278400\",\"1711500\",\"949000\",\"189300\",\"844800\",\"6739500\",\"3010700\",\"377000\",\"1468400\",\"1028000\",\"172500\",\"833700\",\"6890200\""
[8] "\"Arkansas\",\"1128800\",\"155600\",\"600800\",\"515200\",\"67600\",\"436800\",\"2904800\",\"1176500\",\"231700\",\"639200\",\"479400\",\"82000\",\"287200\",\"2896000\",\"1293700\",\"200200\",\"641400\",\"484500\",\"63700\",\"268400\",\"2953000\",\"1290900\",\"252900\",\"618600\",\"490000\",\"67500\",\"225500\",\"2945300\""
[9] "\"California\",\"17747300\",\"1986400\",\"8344800\",\"3828500\",\"675400\",\"5594100\",\"38176400\",\"17703700\",\"2778800\",\"9618800\",\"4049000\",\"634400\",\"3916700\",\"38701300\",\"17718300\",\"3444200\",\"10138100\",\"4080100\",\"752700\",\"2980600\",\"39113900\",\"18116200\",\"3195400\",\"9853800\",\"4436000\",\"556100\",\"3030800\",\"39188300\""
[10] "\"Colorado\",\"2852500\",\"426300\",\"697300\",\"549700\",\"118100\",\"654000\",\"5297800\",\"2489400\",\"397900\",\"1053700\",\"619500\",\"214000\",\"602900\",\"5377400\",\"2706000\",\"346900\",\"1036600\",\"708000\",\"148000\",\"475700\",\"5421300\",\"2872600\",\"370000\",\"855800\",\"692400\",\"190100\",\"528400\",\"5509200\"" It looks like the first two lines are descriptive and are not useful. We will tell R to skip reading these in using the skip argument in read_csv(). The third line looks like it contains the column names and starting on the fourth line is where the data starts.
coverage <- read_csv("./data/KFF/healthcare-coverage.csv",
skip = 2, col_names = TRUE)
head(coverage)# A tibble: 6 x 29
Location `2013__Employer` `2013__Non-Grou~ `2013__Medicaid`
<chr> <int> <int> <int>
1 United ~ 155696900 13816000 54919100
2 Alabama 2126500 174200 869700
3 Alaska 364900 24000 95000
4 Arizona 2883800 170800 1346100
5 Arkansas 1128800 155600 600800
6 Califor~ 17747300 1986400 8344800
# ... with 25 more variables: `2013__Medicare` <int>, `2013__Other
# Public` <chr>, `2013__Uninsured` <int>, `2013__Total` <int>,
# `2014__Employer` <int>, `2014__Non-Group` <int>,
# `2014__Medicaid` <int>, `2014__Medicare` <int>, `2014__Other
# Public` <chr>, `2014__Uninsured` <int>, `2014__Total` <int>,
# `2015__Employer` <int>, `2015__Non-Group` <int>,
# `2015__Medicaid` <int>, `2015__Medicare` <int>, `2015__Other
# Public` <chr>, `2015__Uninsured` <int>, `2015__Total` <int>,
# `2016__Employer` <int>, `2016__Non-Group` <dbl>,
# `2016__Medicaid` <int>, `2016__Medicare` <int>, `2016__Other
# Public` <chr>, `2016__Uninsured` <int>, `2016__Total` <int>tail(coverage)# A tibble: 6 x 29
Location `2013__Employer` `2013__Non-Grou~ `2013__Medicaid`
<chr> <int> <int> <int>
1 <NA> NA NA NA
2 *Uninsu~ NA NA NA
3 <NA> NA NA NA
4 For exa~ NA NA NA
5 <NA> NA NA NA
6 *N/A*: ~ NA NA NA
# ... with 25 more variables: `2013__Medicare` <int>, `2013__Other
# Public` <chr>, `2013__Uninsured` <int>, `2013__Total` <int>,
# `2014__Employer` <int>, `2014__Non-Group` <int>,
# `2014__Medicaid` <int>, `2014__Medicare` <int>, `2014__Other
# Public` <chr>, `2014__Uninsured` <int>, `2014__Total` <int>,
# `2015__Employer` <int>, `2015__Non-Group` <int>,
# `2015__Medicaid` <int>, `2015__Medicare` <int>, `2015__Other
# Public` <chr>, `2015__Uninsured` <int>, `2015__Total` <int>,
# `2016__Employer` <int>, `2016__Non-Group` <dbl>,
# `2016__Medicaid` <int>, `2016__Medicare` <int>, `2016__Other
# Public` <chr>, `2016__Uninsured` <int>, `2016__Total` <int>It looks like we now have the right header, but there are a bunch of NAs in the end of the data frame because most of it isn’t useful data.
Let’s take a closer look at the last 30 lines
tail(coverage, n=30)# A tibble: 30 x 29
Location `2013__Employer` `2013__Non-Grou~ `2013__Medicaid`
<chr> <int> <int> <int>
1 Washing~ 3541600 309000 1026800
2 West Vi~ 841300 42600 382500
3 Wiscons~ 3154500 225300 907600
4 Wyoming 305900 19500 74200
5 Notes NA NA NA
6 The maj~ NA NA NA
7 <NA> NA NA NA
8 "In thi~ NA NA NA
9 <NA> NA NA NA
10 Data ex~ NA NA NA
# ... with 20 more rows, and 25 more variables: `2013__Medicare` <int>,
# `2013__Other Public` <chr>, `2013__Uninsured` <int>,
# `2013__Total` <int>, `2014__Employer` <int>, `2014__Non-Group` <int>,
# `2014__Medicaid` <int>, `2014__Medicare` <int>, `2014__Other
# Public` <chr>, `2014__Uninsured` <int>, `2014__Total` <int>,
# `2015__Employer` <int>, `2015__Non-Group` <int>,
# `2015__Medicaid` <int>, `2015__Medicare` <int>, `2015__Other
# Public` <chr>, `2015__Uninsured` <int>, `2015__Total` <int>,
# `2016__Employer` <int>, `2016__Non-Group` <dbl>,
# `2016__Medicaid` <int>, `2016__Medicare` <int>, `2016__Other
# Public` <chr>, `2016__Uninsured` <int>, `2016__Total` <int>It looks like there is a line with a string Notes in it and everything below that line should not be read in. We can use the n_max argument here.
coverage <- read_csv("./data/KFF/healthcare-coverage.csv",
skip = 2, col_names = TRUE)
coverage <- read_csv("./data/KFF/healthcare-coverage.csv",
skip = 2, col_names = TRUE,
n_max = which(coverage$Location == "Notes")-1)
tail(coverage)# A tibble: 6 x 29
Location `2013__Employer` `2013__Non-Grou~ `2013__Medicaid`
<chr> <int> <int> <int>
1 Vermont 317700 26200 123400
2 Virginia 4661600 364800 773200
3 Washing~ 3541600 309000 1026800
4 West Vi~ 841300 42600 382500
5 Wiscons~ 3154500 225300 907600
6 Wyoming 305900 19500 74200
# ... with 25 more variables: `2013__Medicare` <int>, `2013__Other
# Public` <chr>, `2013__Uninsured` <int>, `2013__Total` <int>,
# `2014__Employer` <int>, `2014__Non-Group` <int>,
# `2014__Medicaid` <int>, `2014__Medicare` <int>, `2014__Other
# Public` <chr>, `2014__Uninsured` <int>, `2014__Total` <int>,
# `2015__Employer` <int>, `2015__Non-Group` <int>,
# `2015__Medicaid` <int>, `2015__Medicare` <int>, `2015__Other
# Public` <chr>, `2015__Uninsured` <int>, `2015__Total` <int>,
# `2016__Employer` <int>, `2016__Non-Group` <dbl>,
# `2016__Medicaid` <int>, `2016__Medicare` <int>, `2016__Other
# Public` <chr>, `2016__Uninsured` <int>, `2016__Total` <int>That’s better!
Read in healthcare spending data
Now because we are also going to want to use in healthcare-spending.csv, let’s read it in now.
spending <- read_csv("./data/KFF/healthcare-spending.csv",
skip = 2, col_names = TRUE)
spending <- read_csv("./data/KFF/healthcare-spending.csv",
skip = 2, col_names = TRUE,
n_max = which(spending$Location == "Notes")-1)
tail(spending)# A tibble: 6 x 25
Location `1991__Total He~ `1992__Total He~ `1993__Total He~
<chr> <int> <int> <int>
1 Vermont 1330 1421 1522
2 Virginia 14829 15599 16634
3 Washing~ 12674 13859 14523
4 West Vi~ 4672 5159 5550
5 Wiscons~ 12694 13669 14636
6 Wyoming 1023 1067 1171
# ... with 21 more variables: `1994__Total Health Spending` <int>,
# `1995__Total Health Spending` <int>, `1996__Total Health
# Spending` <int>, `1997__Total Health Spending` <int>, `1998__Total
# Health Spending` <int>, `1999__Total Health Spending` <int>,
# `2000__Total Health Spending` <int>, `2001__Total Health
# Spending` <int>, `2002__Total Health Spending` <int>, `2003__Total
# Health Spending` <int>, `2004__Total Health Spending` <int>,
# `2005__Total Health Spending` <int>, `2006__Total Health
# Spending` <int>, `2007__Total Health Spending` <int>, `2008__Total
# Health Spending` <int>, `2009__Total Health Spending` <int>,
# `2010__Total Health Spending` <int>, `2011__Total Health
# Spending` <int>, `2012__Total Health Spending` <int>, `2013__Total
# Health Spending` <int>, `2014__Total Health Spending` <int>2. Take a glimpse() at your data
One last thing in this section. One way to look at our data would be to use head() or tail(), as we just saw. Another one you might have heard of is the str() function. One you might not have heard of is the glimpse() function. It’s used for a special type of object in R called a tibble. Let’s read the help file to learn more.
?tibble::tibbleIt’s kind of like print() where it shows you columns running down the page. Let’s try it out. If we look at our data, say the coverage data frame, we see that it is not “tidy”:
glimpse(coverage)Observations: 52
Variables: 29
$ Location <chr> "United States", "Alabama", "Alaska", "Ar...
$ `2013__Employer` <int> 155696900, 2126500, 364900, 2883800, 1128...
$ `2013__Non-Group` <int> 13816000, 174200, 24000, 170800, 155600, ...
$ `2013__Medicaid` <int> 54919100, 869700, 95000, 1346100, 600800,...
$ `2013__Medicare` <int> 40876300, 783000, 55200, 842000, 515200, ...
$ `2013__Other Public` <chr> "6295400", "85600", "60600", "N/A", "6760...
$ `2013__Uninsured` <int> 41795100, 724800, 102200, 1223000, 436800...
$ `2013__Total` <int> 313401200, 4763900, 702000, 6603100, 2904...
$ `2014__Employer` <int> 154347500, 2202800, 345300, 2835200, 1176...
$ `2014__Non-Group` <int> 19313000, 288900, 26800, 333500, 231700, ...
$ `2014__Medicaid` <int> 61650400, 891900, 130100, 1639400, 639200...
$ `2014__Medicare` <int> 41896500, 718400, 55300, 911100, 479400, ...
$ `2014__Other Public` <chr> "5985000", "143900", "37300", "N/A", "820...
$ `2014__Uninsured` <int> 32967500, 522200, 100800, 827100, 287200,...
$ `2014__Total` <int> 316159900, 4768000, 695700, 6657200, 2896...
$ `2015__Employer` <int> 155965800, 2218000, 355700, 2766500, 1293...
$ `2015__Non-Group` <int> 21816500, 291500, 22300, 278400, 200200, ...
$ `2015__Medicaid` <int> 62384500, 911400, 128100, 1711500, 641400...
$ `2015__Medicare` <int> 43308400, 719100, 60900, 949000, 484500, ...
$ `2015__Other Public` <chr> "6422300", "174600", "47700", "189300", "...
$ `2015__Uninsured` <int> 28965900, 519400, 90500, 844800, 268400, ...
$ `2015__Total` <int> 318868500, 4833900, 705300, 6739500, 2953...
$ `2016__Employer` <int> 157381500, 2263800, 324400, 3010700, 1290...
$ `2016__Non-Group` <dbl> 21884400, 262400, 20300, 377000, 252900, ...
$ `2016__Medicaid` <int> 62303400, 997000, 145400, 1468400, 618600...
$ `2016__Medicare` <int> 44550200, 761200, 68200, 1028000, 490000,...
$ `2016__Other Public` <chr> "6192200", "128800", "55600", "172500", "...
$ `2016__Uninsured` <int> 28051900, 420800, 96900, 833700, 225500, ...
$ `2016__Total` <int> 320372000, 4834100, 710800, 6890200, 2945...Read the State information using the datasets R package
Since our goal is to get sense of the health expenditure, including healthcare coverage and healthcare spending, across States, it would be nice add some information about each state. Namely, the state abbreviation and state region (i.e. north, south, etc).
For this we use the state dataset in the datasets R package.
Before we begin, let’s look at what states are there:
unique(coverage$Location) [1] "United States" "Alabama" "Alaska"
[4] "Arizona" "Arkansas" "California"
[7] "Colorado" "Connecticut" "Delaware"
[10] "District of Columbia" "Florida" "Georgia"
[13] "Hawaii" "Idaho" "Illinois"
[16] "Indiana" "Iowa" "Kansas"
[19] "Kentucky" "Louisiana" "Maine"
[22] "Maryland" "Massachusetts" "Michigan"
[25] "Minnesota" "Mississippi" "Missouri"
[28] "Montana" "Nebraska" "Nevada"
[31] "New Hampshire" "New Jersey" "New Mexico"
[34] "New York" "North Carolina" "North Dakota"
[37] "Ohio" "Oklahoma" "Oregon"
[40] "Pennsylvania" "Rhode Island" "South Carolina"
[43] "South Dakota" "Tennessee" "Texas"
[46] "Utah" "Vermont" "Virginia"
[49] "Washington" "West Virginia" "Wisconsin"
[52] "Wyoming" We see there are more than 50 states because “United States” and “District of Columbia” are both included.
Let’s look what states are inside the state dataset.
library(datasets)
data(state)
unique(state.name) [1] "Alabama" "Alaska" "Arizona" "Arkansas"
[5] "California" "Colorado" "Connecticut" "Delaware"
[9] "Florida" "Georgia" "Hawaii" "Idaho"
[13] "Illinois" "Indiana" "Iowa" "Kansas"
[17] "Kentucky" "Louisiana" "Maine" "Maryland"
[21] "Massachusetts" "Michigan" "Minnesota" "Mississippi"
[25] "Missouri" "Montana" "Nebraska" "Nevada"
[29] "New Hampshire" "New Jersey" "New Mexico" "New York"
[33] "North Carolina" "North Dakota" "Ohio" "Oklahoma"
[37] "Oregon" "Pennsylvania" "Rhode Island" "South Carolina"
[41] "South Dakota" "Tennessee" "Texas" "Utah"
[45] "Vermont" "Virginia" "Washington" "West Virginia"
[49] "Wisconsin" "Wyoming" Ah, ok. So let’s start by dealing with DC as a special case.
state.abb <- c(state.abb, "DC")
state.region <- as.factor(c(as.character(state.region), "South"))
state.name <- c(state.name, "District of Columbia")We will deal with the “United States” in the next section.
Data Wrangling
What is “Tidy Data”?
Glance at “Tidy Data”
A subset of the data analysis process can be thought about in the following way:
where each of these steps need their own tools and software to complete.
After we import the data into R, if we are going to take advantage of the “tidyverse”, this means we need to transform the data into a form that is “tidy”. If you recall, in tidy data:
- Each variable forms a column.
- Each observation forms a row.
- Each type of observational unit forms a table.
For example, consider the following dataset:
Here:
- each row represents one company (row names are companies)
- each column represent one time point
- the stock prices are defined for each row/column pair
Alternatively, a data set can be structured in the following way:
- each row represents one time point (but no row names)
- the first column defines the time variable and the last three columns contain the stock prices for three companies
In both cases, the data is the same, but the structure is different. This can be frustrating to deal with as an analyst because the meaning of the values (rows and columns) in the two data sets are different. Providing a standardized way of organizing values within a data set would alleviate a major portion of this frustration.
For motivation, a tidy version of the stock data we looked at above looks like this: (we’ll learn how the functions work in just a moment)
In this “tidy” data set, we have three columns representing three variables (time, company name and stock price). Every row represents contains one stock price from a particular time and for a specific company.
If we consider our coverage dataframe, we see it is also not in a tidy format. Each row contains information about the coverage level by Location across years and types of coverage.
coverage[1:5, 1:5]# A tibble: 5 x 5
Location `2013__Employer` `2013__Non-Grou~ `2013__Medicaid`
<chr> <int> <int> <int>
1 United ~ 155696900 13816000 54919100
2 Alabama 2126500 174200 869700
3 Alaska 364900 24000 95000
4 Arizona 2883800 170800 1346100
5 Arkansas 1128800 155600 600800
# ... with 1 more variable: `2013__Medicare` <int>Now, let’s use the tidyr R package to transform our data into a tidy format.
The tidyr R package
1. What is the tidyr R package ?
tidyr is an R package that transforms data sets to a tidy format.
This package is installed and loaded when you load the tidyverse using library(tidyverse). However, you can also just load the library by itself.
library(tidyr)The main functions in tidyr are:
tidyr functions |
Description |
|---|---|
gather() |
takes multiple columns, and gathers them into key-value pairs, making “wide” data longer |
separate() |
turns a single character column into multiple columns, making “long” data wider |
spread() |
spread rows into multiple columns, transforming “long” data into “wide” format |
We’ll explore what it means to go between a “wide” and “long” data format using gather() , separate(), and spread().
A tidyr cheatsheet for the functions in the tidyr package can be found on RStudio’s website:
2. Convert data from wide format to long format using gather()
Let’s start by looking at the gather() help file
?gatherThis function gathers multiple columns and collapses them into new key-value pairs. This transform data from wide format into a long format.
- The
keyis the name of the new column that you are creating which contains the values of the column headings that you are gathering - The
valueis the name of the new column that will contain the values themselves - The third argument defines the columns to gather
For example, here we create a column titled year_type and coverage. We also want to keep the Location column as it is because it also contains observational level data.
coverage <- gather(coverage, "year_type", "tot_coverage", -Location)
coverage# A tibble: 1,456 x 3
Location year_type tot_coverage
<chr> <chr> <chr>
1 United States 2013__Employer 155696900
2 Alabama 2013__Employer 2126500
3 Alaska 2013__Employer 364900
4 Arizona 2013__Employer 2883800
5 Arkansas 2013__Employer 1128800
6 California 2013__Employer 17747300
7 Colorado 2013__Employer 2852500
8 Connecticut 2013__Employer 2030500
9 Delaware 2013__Employer 473700
10 District of Columbia 2013__Employer 324300
# ... with 1,446 more rowsNow we see each row contains one observation. Namely, a Location, a year_type and coverage. It would be nice to separate out the information in the year_type column into two columns. We can implement same techniques to the healthcare spending dataset.
Convert healthcare spending data to a long format (tidy format)
Let’s do the same for the spending data. In this case I will use year and spending for the key and value. We also want to keep Location like before.
spending <- gather(spending, "year", "tot_spending", -Location)
spending# A tibble: 1,248 x 3
Location year tot_spending
<chr> <chr> <int>
1 United States 1991__Total Health Spending 675896
2 Alabama 1991__Total Health Spending 10393
3 Alaska 1991__Total Health Spending 1458
4 Arizona 1991__Total Health Spending 9269
5 Arkansas 1991__Total Health Spending 5632
6 California 1991__Total Health Spending 81438
7 Colorado 1991__Total Health Spending 8460
8 Connecticut 1991__Total Health Spending 10950
9 Delaware 1991__Total Health Spending 1938
10 District of Columbia 1991__Total Health Spending 2800
# ... with 1,238 more rowsWe will explore how to do that in the Data Wrangling section below. For now let’s learn more about the tidyr package.
3. Convert data from long format to wide format using spread()
In contrast to gathering multiple columns into key-value pairs, we can spread a key-value pair across multiple columns.
The function spread() does just that. It transforms data from a long format into a wide format.
- The
keyis the name of the column in your data set that contains the values of the column headings that you are spreading across multiple columns - The
valueis the name of the column that contains the values for the multiple columns
spread(coverage, year_type, tot_coverage)# A tibble: 52 x 29
Location `2013__Employer` `2013__Medicaid` `2013__Medicare`
<chr> <chr> <chr> <chr>
1 Alabama 2126500 869700 783000
2 Alaska 364900 95000 55200
3 Arizona 2883800 1346100 842000
4 Arkansas 1128800 600800 515200
5 Califor~ 17747300 8344800 3828500
6 Colorado 2852500 697300 549700
7 Connect~ 2030500 532000 475300
8 Delaware 473700 192700 141300
9 Distric~ 324300 174900 59900
10 Florida 8023400 3190900 3108800
# ... with 42 more rows, and 25 more variables: `2013__Non-Group` <chr>,
# `2013__Other Public` <chr>, `2013__Total` <chr>,
# `2013__Uninsured` <chr>, `2014__Employer` <chr>,
# `2014__Medicaid` <chr>, `2014__Medicare` <chr>,
# `2014__Non-Group` <chr>, `2014__Other Public` <chr>,
# `2014__Total` <chr>, `2014__Uninsured` <chr>, `2015__Employer` <chr>,
# `2015__Medicaid` <chr>, `2015__Medicare` <chr>,
# `2015__Non-Group` <chr>, `2015__Other Public` <chr>,
# `2015__Total` <chr>, `2015__Uninsured` <chr>, `2016__Employer` <chr>,
# `2016__Medicaid` <chr>, `2016__Medicare` <chr>,
# `2016__Non-Group` <chr>, `2016__Other Public` <chr>,
# `2016__Total` <chr>, `2016__Uninsured` <chr>In the real world, analyzing data rarely involves data that can be easily imported and ready for analysis. According to Wikipedia:
Data munging or data wrangling is loosely the process of manually converting or mapping data from one “raw” form into another format that allows for more convenient consumption of the data with the help of semi-automated tools.
As you will see in class, one of the most time-consuming aspects of the data analysis process is “data wrangling”. This is also is a trendy term for cleaning up a messy data set.
R provides incredibly powerful and flexible language for data wrangling. However, the syntax is somewhat hard to get used to. We will therefore introducing a package that makes the syntax much more like the English language. This package is dplyr.
The dplyr R package
1. What is the dplyr R package ?
dplyr is a powerful R-package to transform and summarize tabular data with rows and columns.
The package contains a set of functions (or “verbs”) to perform common data manipulation operations such as filtering for rows, selecting specific columns, re-ordering rows, adding new columns and summarizing data.
In addition, dplyr contains a useful function to perform another common task which is the is the “split-apply-combine” concept. We will discuss that in a little bit.
2. Compare dplyr R package compare with base functions R
If you are familiar with R, you are probably familiar with base R functions such as split(), subset(), apply(), sapply(), lapply(), tapply() and aggregate(). Compared to base functions in R, the functions in dplyr are easier to work with, are more consistent in the syntax and are targeted for data analysis around data frames instead of just vectors.
The important dplyr verbs to remember are:
dplyr verbs |
Description |
|---|---|
select() |
select columns |
filter() |
filter rows |
arrange() |
re-order or arrange rows |
mutate() |
create new columns |
summarize() |
summarize values |
group_by() |
allows for group operations in the “split-apply-combine” concept |
3. Pipe operator: %>%
Before we go any futher, let’s introduce the pipe operator: %>%. In our stocks example, we briefly saw this symbol. It is called the pipe operator. dplyr imports this operator from another package (magrittr) see help file here. This operator allows you to pipe the output from one function to the input of another function. Instead of nesting functions (reading from the inside to the outside), the idea of of piping is to read the functions from left to right.
Now in stocks example, we pipe the stocks data frame to the function that will gather multiple columns into key-value pairs.
dplyr verbs in action: separate(), unite(), …
First, let’s separate the year_type column in the coverage dataset to two columns: year and health coverage type.
To do this, we will use the separate() function in the tidyr package.
Note:
separate()= separate one column into multiple columnsunite()= unite multiple columns into one
Learn separate() and unite() in the spending dataset
coverage %>%
separate(year_type, sep="__",
into=c("year", "type"))# A tibble: 1,456 x 4
Location year type tot_coverage
<chr> <chr> <chr> <chr>
1 United States 2013 Employer 155696900
2 Alabama 2013 Employer 2126500
3 Alaska 2013 Employer 364900
4 Arizona 2013 Employer 2883800
5 Arkansas 2013 Employer 1128800
6 California 2013 Employer 17747300
7 Colorado 2013 Employer 2852500
8 Connecticut 2013 Employer 2030500
9 Delaware 2013 Employer 473700
10 District of Columbia 2013 Employer 324300
# ... with 1,446 more rowsWe see that we now have two columns, except the year column was converted to a character. If we look at the help file ?separate, we see we can use the convert=TRUE argument to convert the character to an integer.
coverage <-
coverage %>%
separate(year_type, sep="__",
into=c("year", "type"),
convert = TRUE)
coverage# A tibble: 1,456 x 4
Location year type tot_coverage
<chr> <int> <chr> <chr>
1 United States 2013 Employer 155696900
2 Alabama 2013 Employer 2126500
3 Alaska 2013 Employer 364900
4 Arizona 2013 Employer 2883800
5 Arkansas 2013 Employer 1128800
6 California 2013 Employer 17747300
7 Colorado 2013 Employer 2852500
8 Connecticut 2013 Employer 2030500
9 Delaware 2013 Employer 473700
10 District of Columbia 2013 Employer 324300
# ... with 1,446 more rowsNext, we see that the tot_coverage column is also a character. Gah!
Let’s fix that. We can use the mutate_at() function to do this. We are asking R to take tot_coverage column and convert it to an integer and then replace the old column with the new converted column
coverage <-
coverage %>%
mutate_at("tot_coverage", as.integer)
# Add the abbreviation of States
coverage$abb <- state.abb[match(coverage$Location, state.name)]
coverage$region <- state.region[match(coverage$Location, state.name)]
coverage# A tibble: 1,456 x 6
Location year type tot_coverage abb region
<chr> <int> <chr> <int> <chr> <fct>
1 United States 2013 Employer 155696900 <NA> <NA>
2 Alabama 2013 Employer 2126500 AL South
3 Alaska 2013 Employer 364900 AK West
4 Arizona 2013 Employer 2883800 AZ West
5 Arkansas 2013 Employer 1128800 AR South
6 California 2013 Employer 17747300 CA West
7 Colorado 2013 Employer 2852500 CO West
8 Connecticut 2013 Employer 2030500 CT Northeast
9 Delaware 2013 Employer 473700 DE South
10 District of Columbia 2013 Employer 324300 DC South
# ... with 1,446 more rowsThe coverage data looks good now. We see that there are different years and different types of healthcare coverage.
Also, you may want to link the coverage data with our location information.
# Add the abbreviation of States
coverage$abb <- state.abb[match(coverage$Location, state.name)]
coverage$region <- state.region[match(coverage$Location, state.name)]
coverage# A tibble: 1,456 x 6
Location year type tot_coverage abb region
<chr> <int> <chr> <int> <chr> <fct>
1 United States 2013 Employer 155696900 <NA> <NA>
2 Alabama 2013 Employer 2126500 AL South
3 Alaska 2013 Employer 364900 AK West
4 Arizona 2013 Employer 2883800 AZ West
5 Arkansas 2013 Employer 1128800 AR South
6 California 2013 Employer 17747300 CA West
7 Colorado 2013 Employer 2852500 CO West
8 Connecticut 2013 Employer 2030500 CT Northeast
9 Delaware 2013 Employer 473700 DE South
10 District of Columbia 2013 Employer 324300 DC South
# ... with 1,446 more rows(*) What is the range of years and types of healthcare in the coverage dataset?
table(coverage$type, coverage$year)
2013 2014 2015 2016
Employer 52 52 52 52
Medicaid 52 52 52 52
Medicare 52 52 52 52
Non-Group 52 52 52 52
Other Public 52 52 52 52
Total 52 52 52 52
Uninsured 52 52 52 52Implement separate() and unite() in the spending dataset
Next, we will look at the spending data. We see the year column has information that we do not want. We only care about the year.
spending# A tibble: 1,248 x 3
Location year tot_spending
<chr> <chr> <int>
1 United States 1991__Total Health Spending 675896
2 Alabama 1991__Total Health Spending 10393
3 Alaska 1991__Total Health Spending 1458
4 Arizona 1991__Total Health Spending 9269
5 Arkansas 1991__Total Health Spending 5632
6 California 1991__Total Health Spending 81438
7 Colorado 1991__Total Health Spending 8460
8 Connecticut 1991__Total Health Spending 10950
9 Delaware 1991__Total Health Spending 1938
10 District of Columbia 1991__Total Health Spending 2800
# ... with 1,238 more rowsLet’s use the separate() function with convert=TRUE to separate the year column into columns. Then, we introduce another dplyr action verb: select().
The two most basic functions are select() and filter() which selects columns and filters rows, respectively.
4. Select columns using select()
In the separate() function, we create two new columns called year and name. Then, we ask to return all the columns, except name. To select all the columns except a specific column, use the “-” (subtraction) operator (also known as negative indexing).
spending <-
spending %>%
separate(year, sep="__", into=c("year", "name"), convert = TRUE) %>%
select(-name)
spending# A tibble: 1,248 x 3
Location year tot_spending
<chr> <int> <int>
1 United States 1991 675896
2 Alabama 1991 10393
3 Alaska 1991 1458
4 Arizona 1991 9269
5 Arkansas 1991 5632
6 California 1991 81438
7 Colorado 1991 8460
8 Connecticut 1991 10950
9 Delaware 1991 1938
10 District of Columbia 1991 2800
# ... with 1,238 more rowsThe function select() is much more powerful though. To select a range of columns by name, use the “:” (colon) operator
coverage %>%
select(year:type)# A tibble: 1,456 x 2
year type
<int> <chr>
1 2013 Employer
2 2013 Employer
3 2013 Employer
4 2013 Employer
5 2013 Employer
6 2013 Employer
7 2013 Employer
8 2013 Employer
9 2013 Employer
10 2013 Employer
# ... with 1,446 more rowsTo select all columns that start with the character string “t”, use the function starts_with()
coverage %>%
select(starts_with("t"))# A tibble: 1,456 x 2
type tot_coverage
<chr> <int>
1 Employer 155696900
2 Employer 2126500
3 Employer 364900
4 Employer 2883800
5 Employer 1128800
6 Employer 17747300
7 Employer 2852500
8 Employer 2030500
9 Employer 473700
10 Employer 324300
# ... with 1,446 more rowsSome additional options to select columns based on a specific criteria include
ends_with()= Select columns that end with a character stringcontains()= Select columns that contain a character stringmatches()= Select columns that match a regular expressionone_of()= Select columns names that are from a group of names
5. Select rows using filter()
Let’s say we want to know how many peopled had health insurance coverage in Maryland?
First, we can filter the rows for years in 2007.
coverage %>%
filter(Location == "Maryland")# A tibble: 28 x 6
Location year type tot_coverage abb region
<chr> <int> <chr> <int> <chr> <fct>
1 Maryland 2013 Employer 3172400 MD South
2 Maryland 2013 Non-Group 320800 MD South
3 Maryland 2013 Medicaid 889800 MD South
4 Maryland 2013 Medicare 751500 MD South
5 Maryland 2013 Other Public 124400 MD South
6 Maryland 2013 Uninsured 682000 MD South
7 Maryland 2013 Total 5940900 MD South
8 Maryland 2014 Employer 3558800 MD South
9 Maryland 2014 Non-Group 361700 MD South
10 Maryland 2014 Medicaid 807900 MD South
# ... with 18 more rowsNote: you can use the boolean operators (e.g. >, <, >=, <=, !=, %in%) to create logical tests.
For example, if we wanted only years after 2014, we can add a second criteria:
coverage %>%
filter(Location == "Maryland",
year > 2014)# A tibble: 14 x 6
Location year type tot_coverage abb region
<chr> <int> <chr> <int> <chr> <fct>
1 Maryland 2015 Employer 3431400 MD South
2 Maryland 2015 Non-Group 371400 MD South
3 Maryland 2015 Medicaid 856800 MD South
4 Maryland 2015 Medicare 705500 MD South
5 Maryland 2015 Other Public 141200 MD South
6 Maryland 2015 Uninsured 394300 MD South
7 Maryland 2015 Total 5900500 MD South
8 Maryland 2016 Employer 3210600 MD South
9 Maryland 2016 Non-Group 443000 MD South
10 Maryland 2016 Medicaid 926300 MD South
11 Maryland 2016 Medicare 827000 MD South
12 Maryland 2016 Other Public 153800 MD South
13 Maryland 2016 Uninsured 372100 MD South
14 Maryland 2016 Total 5932800 MD South (*) Has the number of uninsured has increased or decreased in Maryland between 2013 and 2016?
coverage %>%
filter(Location == "Maryland",
type == "Uninsured")# A tibble: 4 x 6
Location year type tot_coverage abb region
<chr> <int> <chr> <int> <chr> <fct>
1 Maryland 2013 Uninsured 682000 MD South
2 Maryland 2014 Uninsured 343000 MD South
3 Maryland 2015 Uninsured 394300 MD South
4 Maryland 2016 Uninsured 372100 MD South What happened between 2013 and 2014?
6. Arrange or re-order rows using arrange()
Now, let’s say we want to see which states has the To arrange (or re-order) rows by a particular column such as the population, list the name of the column you want to arrange the rows by
coverage %>%
arrange(tot_coverage)# A tibble: 1,456 x 6
Location year type tot_coverage abb region
<chr> <int> <chr> <int> <chr> <fct>
1 Vermont 2013 Other Public 9900 VT Northeast
2 Vermont 2014 Other Public 9900 VT Northeast
3 Rhode Island 2013 Other Public 12100 RI Northeast
4 Wyoming 2014 Other Public 13600 WY West
5 Delaware 2013 Other Public 13800 DE South
6 Vermont 2016 Other Public 14600 VT Northeast
7 New Hampshire 2013 Other Public 15100 NH Northeast
8 Wyoming 2016 Other Public 16400 WY West
9 Vermont 2015 Other Public 16500 VT Northeast
10 North Dakota 2014 Other Public 17300 ND North Central
# ... with 1,446 more rows(*) In 2016, what were the top three states with the largest Employer type of healthcare coverage?
Hint: use the desc() function inside of arrange() to order rows in a descending order.
coverage %>%
filter(Location != "United States", year == 2016, type == "Employer") %>%
arrange(desc(tot_coverage)) %>%
head(n=3)# A tibble: 3 x 6
Location year type tot_coverage abb region
<chr> <int> <chr> <int> <chr> <fct>
1 California 2016 Employer 18116200 CA West
2 Texas 2016 Employer 13607200 TX South
3 New York 2016 Employer 9767500 NY Northeast7. Join two datasets using join()
Here, we’re gioing to demonstrate how to join two datasets using series of join() function, including left_join(), right_join(), inner_join(), …
Up until now, we have been working with three datasets coverage and spending separtely. Next, we will combine these together.
If we want to combine, say, coverage and spending together, we have to decide a few things. Both share a Location column and a year column. However, the range of years is different between datasets.
table(coverage$year)
2013 2014 2015 2016
364 364 364 364 table(spending$year)
1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005
52 52 52 52 52 52 52 52 52 52 52 52 52 52 52
2006 2007 2008 2009 2010 2011 2012 2013 2014
52 52 52 52 52 52 52 52 52 Do we want a dataset with all the years available or only a portion of it? Because there is spending information from 1991-2014, and coverage information from 2013-2016.
dplyr has a list of join functions that are useful to combine datasets. To read more about them, Jenny Bryan has created a nice cheatsheet.
If we look at the help file
?dplyr::joinWe see there are several options for us to pick from. Let’s try one out. We’ll start with left_join() and see what that does.
hc <- left_join(coverage, spending, by = c("Location", "year"))
head(hc)# A tibble: 6 x 7
Location year type tot_coverage abb region tot_spending
<chr> <int> <chr> <int> <chr> <fct> <int>
1 United States 2013 Employer 155696900 <NA> <NA> 2435624
2 Alabama 2013 Employer 2126500 AL South 33788
3 Alaska 2013 Employer 364900 AK West 7684
4 Arizona 2013 Employer 2883800 AZ West 41481
5 Arkansas 2013 Employer 1128800 AR South 20500
6 California 2013 Employer 17747300 CA West 278168tail(hc)# A tibble: 6 x 7
Location year type tot_coverage abb region tot_spending
<chr> <int> <chr> <int> <chr> <fct> <int>
1 Vermont 2016 Total 622500 VT Northeast NA
2 Virginia 2016 Total 8175000 VA South NA
3 Washington 2016 Total 7297300 WA West NA
4 West Virginia 2016 Total 1814100 WV South NA
5 Wisconsin 2016 Total 5766100 WI North Central NA
6 Wyoming 2016 Total 571700 WY West NAWhat did it do? We see that the new hc dataset includes all the years from 2013-2016 (as that is the range of years in coverage), but because the spending dataset only goes to 2014, the tot_spending is reported as NA for years 2015 and 2016.
What about a right_join()?
hc <- right_join(coverage, spending, by = c("Location", "year"))
head(hc)# A tibble: 6 x 7
Location year type tot_coverage abb region tot_spending
<chr> <int> <chr> <int> <chr> <fct> <int>
1 United States 1991 <NA> NA <NA> <NA> 675896
2 Alabama 1991 <NA> NA <NA> <NA> 10393
3 Alaska 1991 <NA> NA <NA> <NA> 1458
4 Arizona 1991 <NA> NA <NA> <NA> 9269
5 Arkansas 1991 <NA> NA <NA> <NA> 5632
6 California 1991 <NA> NA <NA> <NA> 81438tail(hc)# A tibble: 6 x 7
Location year type tot_coverage abb region tot_spending
<chr> <int> <chr> <int> <chr> <fct> <int>
1 Wyoming 2014 Non-Group 31600 WY West 4856
2 Wyoming 2014 Medicaid 54900 WY West 4856
3 Wyoming 2014 Medicare 65600 WY West 4856
4 Wyoming 2014 Other Public 13600 WY West 4856
5 Wyoming 2014 Uninsured 58100 WY West 4856
6 Wyoming 2014 Total 572000 WY West 4856Here, we see every row in the spending dataset is there, but with NAs for the years that there was no coverage data.
There is also a full_join() and inner_join(). If we want the intersection of years from coverage and spending (meaning only 2013 and 2014), we should use inner_join().
hc <- inner_join(coverage, spending, by = c("Location", "year"))
head(hc)# A tibble: 6 x 7
Location year type tot_coverage abb region tot_spending
<chr> <int> <chr> <int> <chr> <fct> <int>
1 United States 2013 Employer 155696900 <NA> <NA> 2435624
2 Alabama 2013 Employer 2126500 AL South 33788
3 Alaska 2013 Employer 364900 AK West 7684
4 Arizona 2013 Employer 2883800 AZ West 41481
5 Arkansas 2013 Employer 1128800 AR South 20500
6 California 2013 Employer 17747300 CA West 278168tail(hc)# A tibble: 6 x 7
Location year type tot_coverage abb region tot_spending
<chr> <int> <chr> <int> <chr> <fct> <int>
1 Vermont 2014 Total 617000 VT Northeast 6389
2 Virginia 2014 Total 8258800 VA South 62847
3 Washington 2014 Total 7085000 WA West 55819
4 West Virginia 2014 Total 1825500 WV South 17491
5 Wisconsin 2014 Total 5747200 WI North Central 50109
6 Wyoming 2014 Total 572000 WY West 4856Yes, that’s what we want!
Next, if we are only interested in looking at US states, we can remove the rows corresponding to the Location == "United States"
hc <- hc %>%
filter(Location != "United States")Another problem is that inside our hc dataset, we have seen there are multiple types of healthcare coverage
table(hc$type)
Employer Medicaid Medicare Non-Group Other Public
102 102 102 102 102
Total Uninsured
102 102 The total type is not really a formal type of healthcare coverage. It really represents just the total number of people in the state. This is useful information and we can include it as a column called tot_pop. How can we do this?
Well, one way would be to use the join functions again in dplyr.
pop <- hc %>%
filter(type == "Total") %>%
select(Location, year, tot_coverage)
pop# A tibble: 102 x 3
Location year tot_coverage
<chr> <int> <int>
1 Alabama 2013 4763900
2 Alaska 2013 702000
3 Arizona 2013 6603100
4 Arkansas 2013 2904800
5 California 2013 38176400
6 Colorado 2013 5297800
7 Connecticut 2013 3578900
8 Delaware 2013 909300
9 District of Columbia 2013 652100
10 Florida 2013 19429000
# ... with 92 more rowshc <- hc %>%
filter(type != "Total") %>%
left_join(pop, by = c("Location", "year")) %>%
rename(tot_coverage = tot_coverage.x, tot_pop = tot_coverage.y)
hc# A tibble: 612 x 8
Location year type tot_coverage abb region tot_spending tot_pop
<chr> <int> <chr> <int> <chr> <fct> <int> <int>
1 Alabama 2013 Emplo~ 2126500 AL South 33788 4.76e6
2 Alaska 2013 Emplo~ 364900 AK West 7684 7.02e5
3 Arizona 2013 Emplo~ 2883800 AZ West 41481 6.60e6
4 Arkansas 2013 Emplo~ 1128800 AR South 20500 2.90e6
5 California 2013 Emplo~ 17747300 CA West 278168 3.82e7
6 Colorado 2013 Emplo~ 2852500 CO West 34090 5.30e6
7 Connecticut 2013 Emplo~ 2030500 CT North~ 34223 3.58e6
8 Delaware 2013 Emplo~ 473700 DE South 9038 9.09e5
9 District o~ 2013 Emplo~ 324300 DC South 7443 6.52e5
10 Florida 2013 Emplo~ 8023400 FL South 150547 1.94e7
# ... with 602 more rowsWe can check to make sure that the total is no longer listed as a type of healthcare coverage.
table(hc$type)
Employer Medicaid Medicare Non-Group Other Public
102 102 102 102 102
Uninsured
102 We are now ready to try answering our first question that we asked:
- Is there a relationship between healthcare coverage and healthcare spending in the United States?
Let’s pick out the type==Employer and year==2013.
hc.employer.2013 <- hc %>%
filter(type == "Employer", year == "2013")
plot(hc.employer.2013$tot_spending,
hc.employer.2013$tot_coverage, log = "xy",
xlab = "spending", ylab = "coverage")
We see there is a strong relationship. However, we also see that healthcare coverage and spending is also strongly related to population size
par(mfrow=c(1,2))
plot(hc.employer.2013$tot_pop,
hc.employer.2013$tot_coverage, log = "xy",
xlab = "population size", ylab = "coverage")
plot(hc.employer.2013$tot_pop,
hc.employer.2013$tot_spending, log = "xy",
xlab = "population size", ylab = "spending")
This means we need to take into account the population size of each state when we are comparing the heathcare coverage and spending.
8. Add columns using mutate()
Instead of the absolute number of people who are covered (tot_coverage), we will calculate the proportion of people who are coverage in each state, year and type.
For this, we will use the mutate() function in dplyr.
hc <- hc %>%
mutate(prop_coverage = tot_coverage/tot_pop)
hc# A tibble: 612 x 9
Location year type tot_coverage abb region tot_spending tot_pop
<chr> <int> <chr> <int> <chr> <fct> <int> <int>
1 Alabama 2013 Empl~ 2126500 AL South 33788 4.76e6
2 Alaska 2013 Empl~ 364900 AK West 7684 7.02e5
3 Arizona 2013 Empl~ 2883800 AZ West 41481 6.60e6
4 Arkansas 2013 Empl~ 1128800 AR South 20500 2.90e6
5 Califor~ 2013 Empl~ 17747300 CA West 278168 3.82e7
6 Colorado 2013 Empl~ 2852500 CO West 34090 5.30e6
7 Connect~ 2013 Empl~ 2030500 CT North~ 34223 3.58e6
8 Delaware 2013 Empl~ 473700 DE South 9038 9.09e5
9 Distric~ 2013 Empl~ 324300 DC South 7443 6.52e5
10 Florida 2013 Empl~ 8023400 FL South 150547 1.94e7
# ... with 602 more rows, and 1 more variable: prop_coverage <dbl>We need to add another column to our dataset. We will add the spending per capita (or spending per person) in dollars and name this column spending_capita.
How we will do this?
The tot_spending column is reported in millions (1e6). Therefore, to calculate spending_capita we will need to adjust for this scaling factor to report it on the original scale (just dollars) and then divide by tot_pop.
hc <- hc %>%
mutate(spending_capita = (tot_spending*1e6) / tot_pop)
hc %>% select(prop_coverage, spending_capita)# A tibble: 612 x 2
prop_coverage spending_capita
<dbl> <dbl>
1 0.446 7093.
2 0.520 10946.
3 0.437 6282.
4 0.389 7057.
5 0.465 7286.
6 0.538 6435.
7 0.567 9562.
8 0.521 9940.
9 0.497 11414.
10 0.413 7749.
# ... with 602 more rowsNow we are ready to go back to our first question
- Is there a relationship between healthcare coverage and healthcare spending in the United States?
hc.employer.2013 <- hc %>%
filter(type == "Employer", year == "2013")
plot(hc.employer.2013$spending_capita,
hc.employer.2013$prop_coverage, log = "xy",
xlab = "spending per capita",
ylab = "proportion of Employer coverage")
Yes, it looks like there is a relationship for Employer healthcare coverage in 2013.
We will continue to explore the other types of coverages later on. For now, we get back to to learning more action verbs in dplyr.
Our second question that we were interested in was:
- Which US states spend the most and which spend the least on healthcare? How does the spending distribution change across geographic regions in the United States?
To answer these questions, we need to learn how to calculate summary statistics in our data.
9. Create summaries of columns using summarise()
The summarise() function in dplyr will create summary statistics for a given column in the data frame such as finding the max, min, average. For example, to compute the average spending per capita, we can apply the mean() function to the column spending_captia and call the summary value avg_spending_capita.
hc %>%
summarise(avg_spending_capita = mean(spending_capita))# A tibble: 1 x 1
avg_spending_capita
<dbl>
1 8246.There are many other summary statistics you could consider such sd(), min(), median(), mean(), sum(), n() (returns the length of vector), first() (returns first value in vector), last() (returns last value in vector) and n_distinct() (number of distinct values in vector).
Also note, this is the average across all states, and all years. This is not very informative.
If you recall, our question asked about which states spent the most, so we want an average spending per capita for each state.
For this, we need to introduc another function in dplyr called group_by().
10. Group operations using group_by()
The group_by() verb is and incredibly powerful function in dplyr. As we mentioned before it’s related to concept of “split-apply-combine”.
In our example above, we want to split the data frame by some variable (e.g. Location), apply a function to the individual data frames (mean) and then combine the output back into a summary data frame.
Let’s see how that would look
hc %>%
group_by(Location) %>%
summarise(avg_spending_capita = mean(spending_capita))# A tibble: 51 x 2
Location avg_spending_capita
<chr> <dbl>
1 Alabama 7244.
2 Alaska 11331.
3 Arizona 6397.
4 Arkansas 7324.
5 California 7416.
6 Colorado 6602.
7 Connecticut 9730.
8 Delaware 10127.
9 District of Columbia 11698.
10 Florida 7945.
# ... with 41 more rowsThat’s better. Here we are averaging across the years 2013 and 2014.
(*) What are the top 3 states that have the largestaverage spending per capita? What about the top 3 states with the smallest average spending per capita?
# smallest
hc %>%
group_by(Location) %>%
summarise(avg_spending_capita = mean(spending_capita)) %>%
arrange(avg_spending_capita) %>%
head(n=3)# A tibble: 3 x 2
Location avg_spending_capita
<chr> <dbl>
1 Utah 5842.
2 Arizona 6397.
3 Georgia 6513.# largest
hc %>%
group_by(Location) %>%
summarise(avg_spending_capita = mean(spending_capita)) %>%
arrange(desc(avg_spending_capita)) %>%
head(n=3)# A tibble: 3 x 2
Location avg_spending_capita
<chr> <dbl>
1 District of Columbia 11698.
2 Alaska 11331.
3 Massachusetts 10535.(*) How does the spending distribution change across geographic regions in the United States?
Hint: Calculate the mean and standard deviation of spending per capita for each geographic region in the US.
hc %>%
group_by(region) %>%
summarise(avg_spending_capita = mean(spending_capita),
sd_spending_capita = sd(spending_capita))# A tibble: 4 x 3
region avg_spending_capita sd_spending_capita
<fct> <dbl> <dbl>
1 North Central 8404. 541.
2 Northeast 9592. 546.
3 South 7994. 1273.
4 West 7498. 1329.Another way to visualize distributions is to use boxplots.
Create four boxplots represening the spending per capita distribution for each of the four regions using the boxplot() function in R.
boxplot(hc$spending_capita ~ hc$region)
Now that we have our data in a tidy format, next, we will learn about how to do this using the ggplot2 R package in the tidyverse.
Data Visualization
As you have already seen, there are many functions available in base R that can create plots (e.g. plot(), boxplot()). Others include: hist(), qqplot(), etc. These functions are great because they come with a basic installation of R and can be quite powerful when you need a quick visualization of something when you are exploring data.
We are choosing to introduce ggplot2 because, in our opinion, it’s one of the simplest ways for beginners to create relatively complicated plots that are intuitive and aesthically pleasing.
The ggplot2 R package
The reasons ggplot2 is generally intuitive for beginners is the use of grammar of graphics or the gg in ggplot2. The idea is that you can construct many sentences by learning just a few nouns, adjectives, and verbs. There are specific “words” that we will need to learn and once we do, you will be able to create (or “write”) hundreds of different plots.
The critical part to making graphics using ggplot2 is the data needs to be in a tidy format. Given that we have just spend the last two lectures learning about how to work with tidy data, we are primed to take advantage of all that ggplot2 has to offer!
We will show how it’s easy to pipe tidy data (output) as input to other functions that creates plots. This all works because we are working within the tidyverse.
ggplot2 cheatsheet
The cheatsheet looks like the following: 
1. What is the ggplot() function?
As explained by Hadley Wickham:
the grammar tells us that a statistical graphic is a mapping from data to aesthetic attributes (colour, shape, size) of geometric objects (points, lines, bars). The plot may also contain statistical transformations of the data and is drawn on a specific coordinates system.
ggplot2 Terminology
- ggplot - the main function where you specify the data set and variables to plot (this is where we define the
xandyvariable names) - geoms - geometric objects
- e.g.
geom_point(),geom_bar(),geom_line(),geom_histogram()
- e.g.
- aes - aesthetics
- shape, transparency, color, fill, linetype
- scales - define how your data will be plotted
- continuous, discrete, log, etc
There are three ways to initialize a ggplot() object.
An empty ggplot object
library(ggplot2)
p <- ggplot() A ggplot object associated with a dataset
p <- hc %>%
filter(year==2014) %>%
ggplot()or a ggplot object with a dataset and x and y defined
p <- hc %>%
filter(year==2014) %>%
ggplot(aes(x = spending_capita, y = prop_coverage))p2. Create scatter plots using geom_point()
The function aes() is an aesthetic mapping function inside the ggplot() object. We use this function to specify plot attributes (e.g. x and y variable names) that will not change as we add more layers.
Anything that goes in the ggplot() object becomes a global setting. From there, we use the geom objects to add more layers to the base ggplot() object. These will define what we are interested in illustrating using the data.
If you recall, our first question that we were interested in was
- Is there a relationship between healthcare coverage and healthcare spending in the United States?
Before, we were using base R to create something like this:
hc.employer.2013 <- hc %>%
filter(type == "Employer", year == "2013")
plot(hc.employer.2013$spending_capita,
hc.employer.2013$prop_coverage,
xlab = "spending per capita",
ylab = "coverage proportion")
Let’s re-create this plot with ggplot2 using the geom_point() geometry.
p <- hc %>%
filter(type == "Employer", year == "2013") %>%
ggplot(aes(x = spending_capita, y = prop_coverage))
p + geom_point() +
xlab("spending per capita") +
ylab("coverage proportion")
We used the xlab() and ylab() functions in ggplot2 to specify the x-axis and y-axis labels.
Note, we do not have to assign (<-) the plot to anything:
hc %>%
filter(type == "Employer", year == "2013") %>%
ggplot(aes(x = spending_capita, y = prop_coverage)) +
geom_point() +
xlab("spending per capita") +
ylab("coverage proportion")
It’s also simple to fit a linear regression model and plot it on top of scatter plot using the geom_smooth() (or stat_smooth()) functions.
hc %>%
filter(type == "Employer", year == "2013") %>%
ggplot(aes(x = spending_capita, y = prop_coverage)) +
geom_point() +
xlab("spending per capita") +
ylab("coverage proportion") +
geom_smooth(method = "lm", col = "red")
The standard error bounds are computed and included in the plot.
It would be nice to know which state is represented by which state. For this, we will introduce another geom called geom_text().
3. Add layers of text using geom_text()
In our dataset, we have information about the abbreviation for each state. We could add the abbreviations for each state next to the point on the plot to assess which states have a higher or lower coverage for a given amount of money they spend per capita.
hc %>%
filter(type == "Employer", year == "2013") %>%
ggplot(aes(x = spending_capita, y = prop_coverage)) +
geom_point() +
xlab("spending per capita") +
ylab("coverage proportion") +
geom_smooth(method = "lm", col = "red") +
geom_text(aes(label=abb))
That is cool, but it would be even better if we could nudge the text over a bit. Let’s look at the help file for geom_text():
?ggplot2::geom_textWe see there is an argument called nudge_x and nudge_y. We can use these to nudge the text over a bit so the text is not directly on top of the points.
hc %>%
filter(type == "Employer", year == "2013") %>%
ggplot(aes(x = spending_capita, y = prop_coverage)) +
geom_point() +
xlab("spending per capita") +
ylab("coverage proportion") +
geom_smooth(method = "lm", col = "red") +
geom_text(aes(label=abb), nudge_x = 150)
(*) Color each point (or state) by what region they are from.
## add your code here
hc %>%
filter(type == "Employer", year == "2013") %>%
ggplot(aes(x = spending_capita, y = prop_coverage,
color = region)) +
geom_point() +
xlab("spending per capita") +
ylab("coverage proportion") +
geom_smooth(method = "lm", col = "red") +
geom_text(aes(label=abb), nudge_x = 150)
() Try to explore the package ggrepel, and check if you can improve the quality of visualization using the geom_text_repel* from ggrepel instead of x_nudge from ggplot2.
## add your code here
library(ggrepel)
hc %>%
filter(type == "Employer", year == "2013") %>%
ggplot(aes(x = spending_capita, y = prop_coverage,
color = region)) +
geom_point() +
xlab("spending per capita") +
ylab("coverage proportion") +
geom_smooth(method = "lm", col = "red") +
geom_text_repel(aes(label=abb)) 
4. Facet across a variable using facet_wrap
Ok, getting back to our original question:
- Is there a relationship between healthcare coverage and healthcare spending in the United States?
We saw there was a positive relationship, but this was only for one type of healthcare coverage (Employer) and one year. What about the other types?
For this, we will introduce facets. The idea of faceting is to stratify the data by some variable and make the same plot for each strata.
For example, if we wanted to facet by the type variable, we will add a layer to our ggplot() object using the facet_grid() or facet_wrap() functions. The function expects the row and column variables to be separated by a ~.
hc %>%
filter(year == "2013") %>%
ggplot(aes(x = spending_capita, y = prop_coverage,
color = region)) +
geom_point() +
xlab("spending per capita") +
ylab("coverage proportion") +
geom_smooth(method = "lm", col = "red") +
geom_text_repel(aes(label=abb)) +
facet_wrap(~type)
We see that the proportion of people covered have different scales in the y-axis. Let’s read the help file to see if there is some way to not restrict the y-axis to be the same.
?ggplot2::facet_gridYes, we see there is an argument called scales that can be free_y, (free columns), free_x (free rows), and free (both). Let’s try free_y and look at a different year (year=="2014"):
hc %>%
filter(year == "2014") %>%
ggplot(aes(x = spending_capita, y = prop_coverage,
color = region)) +
geom_point() +
xlab("spending per capita") +
ylab("coverage proportion") +
geom_smooth(method = "lm", col = "red") +
geom_text_repel(aes(label=abb)) +
facet_wrap(~type, scales="free_y")
Given we know Other Public refers to the military or Veterans Adminstration, we can see states like HI, VA, NV have a larger proportion of military or VA Other Public type coverage. While a state like AK has a similar proportion of Other Public coverage, it has a much larger spending per capita.
We also see a negative relationship with the Uninsured type. The more states spend, the less uninsured people in the state.
5. Create boxplots using geom_boxplot()
Next, let’s revisit the second question.
- Which US states spend the most and which spend the least on healthcare? How does the spending distribution change across geographic regions in the United States?
Let’s try making a boxplot with ggplot2. If you recall, the way to do this in base R was:
boxplot(hc$spending_capita ~ hc$region)
Now, we introduce the geom_boxplot() function. Note, we needed to tell ggplot2 what needs to be along the x and y axis in aes().
hc %>%
ggplot(aes(x = region, y = spending_capita)) +
geom_boxplot()
6. Facet by two variables using facet_grid
- Does the relationship between healthcare coverage and healthcare spending in the United States change from 2013 to 2014?
Let’s try faceting by both year and type. Note that we can facet by rows putting a column name before the ~ and facet by columns putting a column name after the ~. We are also using facet_grid() instead of facet_wrap().
p <- hc %>%
ggplot(aes(x = spending_capita, y = prop_coverage,
color = region)) +
geom_point() +
xlab("spending per capita") +
ylab("coverage proportion") +
geom_smooth(method = "lm", col = "red") +
geom_text_repel(aes(label=abb))
p + facet_grid(year~type, scales="free")
p + facet_grid(year~type)
Summary
The total healthcare expenditure is associated with the population. To make a fair comparison, we create “healthcare expenditure per capita.” Further, the exploratory analysis via data visualization showed higher spending in healthcare per capita is positively associated with higher employer coverage proportion and is negatively associated with the porportion of uninsured population across the States.