Section 3.11 Case Studies
So far, we’ve introduced the case studies and read the raw data into R.
Let’s load the raw data that we previously saved using the
here package.
library(here)
load(here::here("data","raw_data", "case_study_1.rda"))
load(here::here("data", "raw_data", "case_study_2.rda"))
#This loads all the data objects that we previously saved in our raw_data directory. Recall that this directory is located within a directory called data that is located within the directory where our project is located.
## here() starts at /Users/carriewright/Documents/GitHub/Coursera/tidyversecourse
Now, we will work to get the data into two tidy formatted datasets that will include the information needed to answer our questions of interest.
Subsection 3.11.1 Case Study #1: Health Expenditures
We’ve already read in the datasets we’ll use for this health expenditures case study, but they’re not yet cleaned and wrangled. So, we’ll do that here!
As a reminder, we’re ultimately interested in answering the following questions with these data:
-
Is there a relationship between health care coverage and health care spending in the United States?
-
How does the spending distribution change across geographic regions in the United States?
-
Does the relationship between health care coverage and health care spending in the United States change from 2013 to 2014?
This means that we’ll need all the data from the variables necessary to answer this question in our tidy dataset.
Subsubsection 3.11.1.1 health care Coverage Data
Let’s remind ourselves before we get to wrangling what data we have when it comes to health care coverage.
coverage
## # A tibble: 52 × 29 ## Location `2013__Employer` `2013__Non-Grou… `2013__Medicaid` `2013__Medicare` ## <chr> <dbl> <dbl> <dbl> <dbl> ## 1 United S… 155696900 13816000 54919100 40876300 ## 2 Alabama 2126500 174200 869700 783000 ## 3 Alaska 364900 24000 95000 55200 ## 4 Arizona 2883800 170800 1346100 842000 ## 5 Arkansas 1128800 155600 600800 515200 ## 6 Californ… 17747300 1986400 8344800 3828500 ## 7 Colorado 2852500 426300 697300 549700 ## 8 Connecti… 2030500 126800 532000 475300 ## 9 Delaware 473700 25100 192700 141300 ## 10 District… 324300 30400 174900 59900 ## # … with 42 more rows, and 24 more variables: 2013__Other Public <chr>, ## # 2013__Uninsured <dbl>, 2013__Total <dbl>, 2014__Employer <dbl>, ## # 2014__Non-Group <dbl>, 2014__Medicaid <dbl>, 2014__Medicare <dbl>, ## # 2014__Other Public <chr>, 2014__Uninsured <dbl>, 2014__Total <dbl>, ## # 2015__Employer <dbl>, 2015__Non-Group <dbl>, 2015__Medicaid <dbl>, ## # 2015__Medicare <dbl>, 2015__Other Public <chr>, 2015__Uninsured <dbl>, ## # 2015__Total <dbl>, 2016__Employer <dbl>, 2016__Non-Group <dbl>, …
At a glance, we see that state-level information is stored in rows (with the exception of the first row, which stores country-level information) with columns corresponding to the amount of money spent on each type of health care, by year.
States Data.
To work with these data, we’ll also want to be able to switch between full state names and two letter abbreviations. There’s data in R available to you for just this purpose!
library(datasets)
data(state)
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"
Before going any further, let’s add some information about Washington, D.C, the nation’s capital, which is not a state, but a territory.
state.abb <- c(state.abb, "DC")
state.region <- as.factor(c(as.character(state.region), "South"))
state.name <- c(state.name, "District of Columbia")
state_data <- tibble(Location = state.name,
abb = state.abb,
region = state.region)
state_data
## # A tibble: 51 × 3 ## Location abb region ## <chr> <chr> <fct> ## 1 Alabama AL South ## 2 Alaska AK West ## 3 Arizona AZ West ## 4 Arkansas AR South ## 5 California CA West ## 6 Colorado CO West ## 7 Connecticut CT Northeast ## 8 Delaware DE South ## 9 Florida FL South ## 10 Georgia GA South ## # … with 41 more rows
If we focus in on the columns within this dataframe, we see that we have a number of different types of health care (i.e. employer, medicare, medicaid, etc.) for each year between 2013 and 2016:
names(coverage)
## [1] "Location" "2013__Employer" "2013__Non-Group" ## [4] "2013__Medicaid" "2013__Medicare" "2013__Other Public" ## [7] "2013__Uninsured" "2013__Total" "2014__Employer" ## [10] "2014__Non-Group" "2014__Medicaid" "2014__Medicare" ## [13] "2014__Other Public" "2014__Uninsured" "2014__Total" ## [16] "2015__Employer" "2015__Non-Group" "2015__Medicaid" ## [19] "2015__Medicare" "2015__Other Public" "2015__Uninsured" ## [22] "2015__Total" "2016__Employer" "2016__Non-Group" ## [25] "2016__Medicaid" "2016__Medicare" "2016__Other Public" ## [28] "2016__Uninsured" "2016__Total"
While a lot of information in here will be helpful, it’s not in a tidy format. This is because, each variable is not in a separate column. For example, each column includes year, the type of coverage and the amount spent by state. We’ll want to use each piece of information separately downstream as we start to visualize and analyze these data. So, let’s work to get these pieces of information separated out now.
To accomplish this, the first thing we’ll have to do is reshape the data, using the
pivot_longer() function from the tidyr package. As a reminder, this function gathers multiple columns and collapses them into new name-value pairs. This transform data from wide format into a long format, where:
-
The first argument defines the columns to gather
-
The
names_toargument is the name of the new column that you are creating which contains the values of the column headings that you are gathering -
The
values_toargument is the name of the new column that will contain the values themselves; you can indicate the name of this column with thevalues_toargument.
Here, we create a column titled
year_type and tot_coverage, storing this newly formatted dataframe back into the variable name coverage. We also want to keep the Location column as it is because it also contains observational level data.
coverage <- coverage %>%
mutate(across(starts_with("20"),
as.integer)) %>% ## Convert all year-based columns to integer
pivot_longer(-Location, ## Use all columns BUT 'Location'
names_to = "year_type",
values_to = "tot_coverage")
coverage
## Warning in mask$eval_all_mutate(quo): NAs introduced by coercion ## Warning in mask$eval_all_mutate(quo): NAs introduced by coercion ## Warning in mask$eval_all_mutate(quo): NAs introduced by coercion ## Warning in mask$eval_all_mutate(quo): NAs introduced by coercion ## # A tibble: 1,456 × 3 ## Location year_type tot_coverage ## <chr> <chr> <int> ## 1 United States 2013__Employer 155696900 ## 2 United States 2013__Non-Group 13816000 ## 3 United States 2013__Medicaid 54919100 ## 4 United States 2013__Medicare 40876300 ## 5 United States 2013__Other Public 6295400 ## 6 United States 2013__Uninsured 41795100 ## 7 United States 2013__Total 313401200 ## 8 United States 2014__Employer 154347500 ## 9 United States 2014__Non-Group 19313000 ## 10 United States 2014__Medicaid 61650400 ## # … with 1,446 more rows
Great! We still have
Location stored in a single column, but we’ve separated out year_type and tot_coverage into their own columns, storing all of the information in a long data format.
Unfortunately, the
year_type column still contains two pieces of information. We’ll want to separate these out to ensure that the data are in a properly tidy format. To do this, we’ll use the separate() function, which allows us to separate out the information stored in a single column into two columns. We’ll also 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 × 4 ## Location year type tot_coverage ## <chr> <int> <chr> <int> ## 1 United States 2013 Employer 155696900 ## 2 United States 2013 Non-Group 13816000 ## 3 United States 2013 Medicaid 54919100 ## 4 United States 2013 Medicare 40876300 ## 5 United States 2013 Other Public 6295400 ## 6 United States 2013 Uninsured 41795100 ## 7 United States 2013 Total 313401200 ## 8 United States 2014 Employer 154347500 ## 9 United States 2014 Non-Group 19313000 ## 10 United States 2014 Medicaid 61650400 ## # … with 1,446 more rows
Perfect! We now have the four columns we wanted, each storing a separate piece of information, and the year column is an integer, as you would want it to be!
Let’s go one step further and add in the state-level abbreviations and region for each row. We’ll utilize our
state datasets that we read in previously to accomplish this! Because we formatted the state data as a tibble, we can simply join it with our coverage dataset to get the state and region information.
coverage <- coverage %>%
left_join(state_data, by = "Location")
coverage
## # A tibble: 1,456 × 6 ## Location year type tot_coverage abb region ## <chr> <int> <chr> <int> <chr> <fct> ## 1 United States 2013 Employer 155696900 <NA> <NA> ## 2 United States 2013 Non-Group 13816000 <NA> <NA> ## 3 United States 2013 Medicaid 54919100 <NA> <NA> ## 4 United States 2013 Medicare 40876300 <NA> <NA> ## 5 United States 2013 Other Public 6295400 <NA> <NA> ## 6 United States 2013 Uninsured 41795100 <NA> <NA> ## 7 United States 2013 Total 313401200 <NA> <NA> ## 8 United States 2014 Employer 154347500 <NA> <NA> ## 9 United States 2014 Non-Group 19313000 <NA> <NA> ## 10 United States 2014 Medicaid 61650400 <NA> <NA> ## # … with 1,446 more rows
Perfect! At this point, each row is an observation and each column stores a single piece of information. This dataset is now in good shape!
Subsubsection 3.11.1.2 health care Spending Data
We’ll have to take a similar approach when it comes to tidying the spending data as it has a similar structure to how the coverage data were stored.
spending
## # A tibble: 52 × 25 ## Location `1991__Total He… `1992__Total He… `1993__Total He… `1994__Total He… ## <chr> <dbl> <dbl> <dbl> <dbl> ## 1 United S… 675896 731455 778684 820172 ## 2 Alabama 10393 11284 12028 12742 ## 3 Alaska 1458 1558 1661 1728 ## 4 Arizona 9269 9815 10655 11364 ## 5 Arkansas 5632 6022 6397 6810 ## 6 Californ… 81438 87949 91963 94245 ## 7 Colorado 8460 9215 9803 10382 ## 8 Connecti… 10950 11635 12081 12772 ## 9 Delaware 1938 2111 2285 2489 ## 10 District… 2800 3098 3240 3255 ## # … with 42 more rows, and 20 more variables: ## # 1995__Total Health Spending <dbl>, 1996__Total Health Spending <dbl>, ## # 1997__Total Health Spending <dbl>, 1998__Total Health Spending <dbl>, ## # 1999__Total Health Spending <dbl>, 2000__Total Health Spending <dbl>, ## # 2001__Total Health Spending <dbl>, 2002__Total Health Spending <dbl>, ## # 2003__Total Health Spending <dbl>, 2004__Total Health Spending <dbl>, ## # 2005__Total Health Spending <dbl>, 2006__Total Health Spending <dbl>, …
Here, we reshape the data using
year and tot_spending for the key and value. We also want to keep Location like before. Then, 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. (This process is also referred to as negative indexing.)
# take spending data from wide to long
spending <- spending %>%
pivot_longer(-Location,
names_to = "year",
values_to = "tot_spending")
# separate year and name columns
spending <- spending %>%
separate(year, sep="__",
into = c("year", "name"),
convert = TRUE) %>%
select(-name)
# look at the data
spending
## # A tibble: 1,248 × 3 ## Location year tot_spending ## <chr> <int> <dbl> ## 1 United States 1991 675896 ## 2 United States 1992 731455 ## 3 United States 1993 778684 ## 4 United States 1994 820172 ## 5 United States 1995 869578 ## 6 United States 1996 917540 ## 7 United States 1997 969531 ## 8 United States 1998 1026103 ## 9 United States 1999 1086280 ## 10 United States 2000 1162035 ## # … with 1,238 more rows
Perfect, we have a tidy dataset and the type of information stored in each column is appropriate for the information being stored in the column!
Subsubsection 3.11.1.3 Join the Data
At this point, we have a
coverage dataset and a spending dataset, but ultimately, we want all of this information in a single tidy data frame. To do this, we’ll have to join the datasets together.
We have to decide what type of join we want to do. For our questions, we only want information from years that are found in both the
coverage and the spending datasets. This means that we’ll want to do an inner_join(). This will keep the data from the intersection of years from coverage and spending (meaning only 2013 and 2014). We’ll store this in a new variable: hc.
# inner join to combine data frames
hc <- inner_join(coverage, spending,
by = c("Location", "year"))
hc
## # A tibble: 728 × 7 ## Location year type tot_coverage abb region tot_spending ## <chr> <int> <chr> <int> <chr> <fct> <dbl> ## 1 United States 2013 Employer 155696900 <NA> <NA> 2435624 ## 2 United States 2013 Non-Group 13816000 <NA> <NA> 2435624 ## 3 United States 2013 Medicaid 54919100 <NA> <NA> 2435624 ## 4 United States 2013 Medicare 40876300 <NA> <NA> 2435624 ## 5 United States 2013 Other Public 6295400 <NA> <NA> 2435624 ## 6 United States 2013 Uninsured 41795100 <NA> <NA> 2435624 ## 7 United States 2013 Total 313401200 <NA> <NA> 2435624 ## 8 United States 2014 Employer 154347500 <NA> <NA> 2562824 ## 9 United States 2014 Non-Group 19313000 <NA> <NA> 2562824 ## 10 United States 2014 Medicaid 61650400 <NA> <NA> 2562824 ## # … with 718 more rows
Great, we’ve combined the information in our datasets. But, we’ve got a bit of extraneous information remaining. For example, we want to look only at the state-level. So, let’s filter out the country-level summary row:
# filter to only include state level
hc <- hc %>%
filter(Location != "United States")
Another problem is that inside our
hc dataset, there are multiple types of health care coverage.
table(hc$type)
## ## Employer Medicaid Medicare Non-Group Other Public Total ## 102 102 102 102 102 102 ## Uninsured ## 102
The "Total" type is not really a formal type of health care 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. To accomplish this, we’ll first store this information in a data frame called pop.
pop <- hc %>%
filter(type == "Total") %>%
select(Location, year, tot_coverage)
pop
## # A tibble: 102 × 3 ## Location year tot_coverage ## <chr> <int> <int> ## 1 Alabama 2013 4763900 ## 2 Alabama 2014 4768000 ## 3 Alaska 2013 702000 ## 4 Alaska 2014 695700 ## 5 Arizona 2013 6603100 ## 6 Arizona 2014 6657200 ## 7 Arkansas 2013 2904800 ## 8 Arkansas 2014 2896000 ## 9 California 2013 38176400 ## 10 California 2014 38701300 ## # … with 92 more rows
We can then, using a
left_join to ensure we keep all of the rows in the hc data frame in tact, add this population level information while simultaneously removing the rows where type is "Total" from the dataset. Finally, we’ll rename the columns to be informative of the information stored within:
# ad population level information
hc <- 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 × 8 ## Location year type tot_coverage abb region tot_spending tot_pop ## <chr> <int> <chr> <int> <chr> <fct> <dbl> <int> ## 1 Alabama 2013 Employer 2126500 AL South 33788 4763900 ## 2 Alabama 2013 Non-Group 174200 AL South 33788 4763900 ## 3 Alabama 2013 Medicaid 869700 AL South 33788 4763900 ## 4 Alabama 2013 Medicare 783000 AL South 33788 4763900 ## 5 Alabama 2013 Other Public 85600 AL South 33788 4763900 ## 6 Alabama 2013 Uninsured 724800 AL South 33788 4763900 ## 7 Alabama 2014 Employer 2202800 AL South 35263 4768000 ## 8 Alabama 2014 Non-Group 288900 AL South 35263 4768000 ## 9 Alabama 2014 Medicaid 891900 AL South 35263 4768000 ## 10 Alabama 2014 Medicare 718400 AL South 35263 4768000 ## # … with 602 more rows
From here, instead of only storing 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, storing this information in prop_coverage.
# add proportion covered
hc <- hc %>%
mutate(prop_coverage = tot_coverage/tot_pop)
hc
## # A tibble: 612 × 9 ## Location year type tot_coverage abb region tot_spending tot_pop ## <chr> <int> <chr> <int> <chr> <fct> <dbl> <int> ## 1 Alabama 2013 Employer 2126500 AL South 33788 4763900 ## 2 Alabama 2013 Non-Group 174200 AL South 33788 4763900 ## 3 Alabama 2013 Medicaid 869700 AL South 33788 4763900 ## 4 Alabama 2013 Medicare 783000 AL South 33788 4763900 ## 5 Alabama 2013 Other Public 85600 AL South 33788 4763900 ## 6 Alabama 2013 Uninsured 724800 AL South 33788 4763900 ## 7 Alabama 2014 Employer 2202800 AL South 35263 4768000 ## 8 Alabama 2014 Non-Group 288900 AL South 35263 4768000 ## 9 Alabama 2014 Medicaid 891900 AL South 35263 4768000 ## 10 Alabama 2014 Medicare 718400 AL South 35263 4768000 ## # … with 602 more rows, and 1 more variable: prop_coverage <dbl>
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. We can again use mutate() to accomplish this:
# get spending capita in dollars
hc <- hc %>%
mutate(spending_capita = (tot_spending*1e6) / tot_pop)
hc
## # A tibble: 612 × 10 ## Location year type tot_coverage abb region tot_spending tot_pop ## <chr> <int> <chr> <int> <chr> <fct> <dbl> <int> ## 1 Alabama 2013 Employer 2126500 AL South 33788 4763900 ## 2 Alabama 2013 Non-Group 174200 AL South 33788 4763900 ## 3 Alabama 2013 Medicaid 869700 AL South 33788 4763900 ## 4 Alabama 2013 Medicare 783000 AL South 33788 4763900 ## 5 Alabama 2013 Other Public 85600 AL South 33788 4763900 ## 6 Alabama 2013 Uninsured 724800 AL South 33788 4763900 ## 7 Alabama 2014 Employer 2202800 AL South 35263 4768000 ## 8 Alabama 2014 Non-Group 288900 AL South 35263 4768000 ## 9 Alabama 2014 Medicaid 891900 AL South 35263 4768000 ## 10 Alabama 2014 Medicare 718400 AL South 35263 4768000 ## # … with 602 more rows, and 2 more variables: prop_coverage <dbl>, ## # spending_capita <dbl>
Yes! At this point we have a single tidy data frame storing all the information we’ll need to answer our questions!
Let’s save our new tidy data for case study #1.
save(hc, file = here::here("data", "tidy_data", "case_study_1_tidy.rda"))
Subsection 3.11.2 Case Study #2: Firearms
For our second case study, we’re interested in the following question: At the state-level, what is the relationship between firearm legislation strength and annual rate of fatal police shootings? Time to wrangle all those many datasets we read in previously!
Subsubsection 3.11.2.1 Census Data
Let’s take a look at the raw data to remind ourselves of what information we have:
census
## # A tibble: 236,844 × 19 ## SUMLEV REGION DIVISION STATE NAME SEX ORIGIN RACE AGE CENSUS2010POP ## <chr> <dbl> <dbl> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> ## 1 040 3 6 01 Alabama 0 0 1 0 37991 ## 2 040 3 6 01 Alabama 0 0 1 1 38150 ## 3 040 3 6 01 Alabama 0 0 1 2 39738 ## 4 040 3 6 01 Alabama 0 0 1 3 39827 ## 5 040 3 6 01 Alabama 0 0 1 4 39353 ## 6 040 3 6 01 Alabama 0 0 1 5 39520 ## 7 040 3 6 01 Alabama 0 0 1 6 39813 ## 8 040 3 6 01 Alabama 0 0 1 7 39695 ## 9 040 3 6 01 Alabama 0 0 1 8 40012 ## 10 040 3 6 01 Alabama 0 0 1 9 42073 ## # … with 236,834 more rows, and 9 more variables: ESTIMATESBASE2010 <dbl>, ## # POPESTIMATE2010 <dbl>, POPESTIMATE2011 <dbl>, POPESTIMATE2012 <dbl>, ## # POPESTIMATE2013 <dbl>, POPESTIMATE2014 <dbl>, POPESTIMATE2015 <dbl>, ## # POPESTIMATE2016 <dbl>, POPESTIMATE2017 <dbl>
These data look reasonably tidy to start; however, the information stored in each column is not particularly clear at a glance. For example, what is a RACE of 1? What does that mean?
Well, if we look at the data dictionary in the document sc-est2017-alldata6.pdf, we learn that:
The key for
SEX is as follows:
-
0 = Total
-
1 = Male
-
2 = Female
The key for
ORIGIN is as follows:
-
0 = Total
-
1 = Not Hispanic
-
2 = Hispanic
The key for
RACE is as follows:
-
1 = White Alone
-
2 = Black or African American Alone
-
3 = American Indian and Alaska Native Alone
-
4 = Asian Alone
-
5 = Native Hawaiian and Other Pacific Islander Alone
-
6 = Two or more races
With that information in mind, we can then use the
dplyr package to filter, group_by, and summarize the data in order to calculate the necessary statistics we’ll need to answer our question.
For each state, we add rows in the column
POPESTIMATE2015 since we are looking at the year 2015. Setting the ORIGIN or SEX equal to 0 ensures we don’t add duplicate data, since 0 is the key for both Hispanic and non Hispanic residents and total male and female residents. We group by each state since all data in this study should be at the state level.
We store each of these pieces of information in its own column within the new dataframe we’ve created
census_stats
# summarize by ethnicity
census_stats <- census %>%
filter(ORIGIN == 0, SEX == 0) %>%
group_by(NAME) %>%
summarize(white = sum(POPESTIMATE2015[RACE == 1])/sum(POPESTIMATE2015)*100,
black = sum(POPESTIMATE2015[RACE == 2])/sum(POPESTIMATE2015)*100)
# add hispanic information
census_stats$hispanic <- census %>%
filter(SEX == 0) %>%
group_by(NAME) %>%
summarize(x = sum(POPESTIMATE2015[ORIGIN == 2])/sum(POPESTIMATE2015[ORIGIN == 0])*100) %>%
pull(x)
# add male information
census_stats$male <- census %>%
filter(ORIGIN == 0) %>%
group_by(NAME) %>%
summarize(x = sum(POPESTIMATE2015[SEX == 1])/sum(POPESTIMATE2015[SEX == 0])*100) %>%
pull(x)
# add total population information
census_stats$total_pop <- census %>%
filter(ORIGIN == 0, SEX == 0 ) %>%
group_by(NAME) %>%
summarize(total = sum(POPESTIMATE2015)) %>%
pull(total)
# lowercase state name for consistency
census_stats$NAME <- tolower(census_stats$NAME)
census_stats
## # A tibble: 51 × 6 ## NAME white black hispanic male total_pop ## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> ## 1 alabama 69.5 26.7 4.13 48.5 4850858 ## 2 alaska 66.5 3.67 6.82 52.4 737979 ## 3 arizona 83.5 4.80 30.9 49.7 6802262 ## 4 arkansas 79.6 15.7 7.18 49.1 2975626 ## 5 california 73.0 6.49 38.7 49.7 39032444 ## 6 colorado 87.6 4.47 21.3 50.3 5440445 ## 7 connecticut 80.9 11.6 15.3 48.8 3593862 ## 8 delaware 70.3 22.5 8.96 48.4 944107 ## 9 district of columbia 44.1 48.5 10.7 47.4 672736 ## 10 florida 77.7 16.9 24.7 48.9 20268567 ## # … with 41 more rows
We can approach the age data similarly, where we get the number of people within each state at each age:
# get state-level age information
age_stats <- census %>%
filter(ORIGIN == 0, SEX == 0) %>%
group_by(NAME, AGE) %>%
summarize(sum_ages = sum(POPESTIMATE2015))
age_stats
## `summarise()` has grouped output by 'NAME'. You can override using the `.groups` argument. ## # A tibble: 4,386 × 3 ## # Groups: NAME [51] ## NAME AGE sum_ages ## <chr> <dbl> <dbl> ## 1 Alabama 0 59080 ## 2 Alabama 1 58738 ## 3 Alabama 2 57957 ## 4 Alabama 3 58800 ## 5 Alabama 4 59329 ## 6 Alabama 5 59610 ## 7 Alabama 6 59977 ## 8 Alabama 7 62282 ## 9 Alabama 8 62175 ## 10 Alabama 9 61249 ## # … with 4,376 more rows
This information is in a long format, but it likely makes more sense to store this information in a wide format, where each column is a different state and each row is an age. To do this:
age_stats <- age_stats %>%
pivot_wider(names_from = "NAME",
values_from = "sum_ages")
age_stats
## # A tibble: 86 × 52 ## AGE Alabama Alaska Arizona Arkansas California Colorado Connecticut ## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> ## 1 0 59080 11253 86653 38453 500834 66222 36414 ## 2 1 58738 11109 86758 38005 499070 66528 36559 ## 3 2 57957 11009 86713 37711 499614 66144 36887 ## 4 3 58800 10756 86914 38381 498536 67065 37745 ## 5 4 59329 10895 87624 38443 510026 68443 38962 ## 6 5 59610 10537 87234 38582 498754 69823 39182 ## 7 6 59977 10352 89215 38630 497444 69691 39871 ## 8 7 62282 10431 93236 40141 516916 71415 41438 ## 9 8 62175 10302 93866 40677 518117 72384 42359 ## 10 9 61249 10055 92531 39836 511610 72086 43032 ## # … with 76 more rows, and 44 more variables: Delaware <dbl>, ## # District of Columbia <dbl>, Florida <dbl>, Georgia <dbl>, Hawaii <dbl>, ## # Idaho <dbl>, Illinois <dbl>, Indiana <dbl>, Iowa <dbl>, Kansas <dbl>, ## # Kentucky <dbl>, Louisiana <dbl>, Maine <dbl>, Maryland <dbl>, ## # Massachusetts <dbl>, Michigan <dbl>, Minnesota <dbl>, Mississippi <dbl>, ## # Missouri <dbl>, Montana <dbl>, Nebraska <dbl>, Nevada <dbl>, ## # New Hampshire <dbl>, New Jersey <dbl>, New Mexico <dbl>, New York <dbl>, …
Now that we’ve made the data easier to work with, we need to find a way to get the median. One method is to take the cumulative sum of each column and then divide all the rows by the last row in each respective column, calculating a percentile/quantile for each age. To do this, we first remove the AGE column, as we don’t want to calculate the median for this column. We then apply the
cumsum() function and an anonymous function using purrr’s map_dfc function. This is a special variation of the map() function that returns a dataframe instead of a list by combining the data by column. But, of course, we do still want the AGE information in there, so we add that column back in using mutate() and then reorder the columns so that AGE is at the front again using select().
age_stats %>%
select(-AGE) %>%
map(cumsum) %>%
map(function(x) x/x[nrow(age_stats)]) %>%
glimpse
## List of 51 ## $ Alabama : num [1:86] 0.0122 0.0243 0.0362 0.0484 0.0606 ... ## $ Alaska : num [1:86] 0.0152 0.0303 0.0452 0.0598 0.0746 ... ## $ Arizona : num [1:86] 0.0127 0.0255 0.0382 0.051 0.0639 ... ## $ Arkansas : num [1:86] 0.0129 0.0257 0.0384 0.0513 0.0642 ... ## $ California : num [1:86] 0.0128 0.0256 0.0384 0.0512 0.0643 ... ## $ Colorado : num [1:86] 0.0122 0.0244 0.0366 0.0489 0.0615 ... ## $ Connecticut : num [1:86] 0.0101 0.0203 0.0306 0.0411 0.0519 ... ## $ Delaware : num [1:86] 0.0116 0.0233 0.0348 0.0466 0.0587 ... ## $ District of Columbia: num [1:86] 0.0146 0.0276 0.0407 0.0533 0.0657 ... ## $ Florida : num [1:86] 0.011 0.0219 0.0328 0.0437 0.0547 ... ## $ Georgia : num [1:86] 0.0128 0.0256 0.0384 0.0514 0.0647 ... ## $ Hawaii : num [1:86] 0.0128 0.0258 0.039 0.0519 0.065 ... ## $ Idaho : num [1:86] 0.0139 0.0274 0.0413 0.0549 0.069 ... ## $ Illinois : num [1:86] 0.0123 0.0245 0.0366 0.0488 0.0611 ... ## $ Indiana : num [1:86] 0.0126 0.0253 0.038 0.0507 0.0635 ... ## $ Iowa : num [1:86] 0.0127 0.0254 0.038 0.0506 0.063 ... ## $ Kansas : num [1:86] 0.0133 0.0268 0.0404 0.054 0.0677 ... ## $ Kentucky : num [1:86] 0.0126 0.0251 0.0376 0.05 0.0624 ... ## $ Louisiana : num [1:86] 0.0137 0.0272 0.0405 0.0536 0.0667 ... ## $ Maine : num [1:86] 0.00949 0.01906 0.02881 0.0386 0.04839 ... ## $ Maryland : num [1:86] 0.0123 0.0245 0.0367 0.049 0.0614 ... ## $ Massachusetts : num [1:86] 0.0106 0.0213 0.0319 0.0427 0.0536 ... ## $ Michigan : num [1:86] 0.0114 0.023 0.0345 0.046 0.0577 ... ## $ Minnesota : num [1:86] 0.0127 0.0256 0.0383 0.0511 0.0639 ... ## $ Mississippi : num [1:86] 0.0128 0.0255 0.0382 0.0512 0.0641 ... ## $ Missouri : num [1:86] 0.0123 0.0248 0.0371 0.0494 0.0618 ... ## $ Montana : num [1:86] 0.0123 0.0244 0.0364 0.0484 0.0604 ... ## $ Nebraska : num [1:86] 0.0141 0.0281 0.042 0.0557 0.0696 ... ## $ Nevada : num [1:86] 0.0125 0.0248 0.0374 0.0498 0.0626 ... ## $ New Hampshire : num [1:86] 0.00932 0.01866 0.02852 0.03814 0.04831 ... ## $ New Jersey : num [1:86] 0.0115 0.0232 0.0349 0.0468 0.0589 ... ## $ New Mexico : num [1:86] 0.0124 0.025 0.0376 0.0504 0.0634 ... ## $ New York : num [1:86] 0.0122 0.0241 0.036 0.0478 0.0598 ... ## $ North Carolina : num [1:86] 0.012 0.024 0.0359 0.0479 0.06 ... ## $ North Dakota : num [1:86] 0.015 0.0296 0.0438 0.0576 0.0709 ... ## $ Ohio : num [1:86] 0.012 0.024 0.036 0.048 0.0601 ... ## $ Oklahoma : num [1:86] 0.0136 0.0272 0.0409 0.0546 0.0683 ... ## $ Oregon : num [1:86] 0.0114 0.0229 0.0344 0.046 0.0577 ... ## $ Pennsylvania : num [1:86] 0.0111 0.0222 0.0333 0.0445 0.0558 ... ## $ Rhode Island : num [1:86] 0.0103 0.0205 0.0309 0.0413 0.0518 ... ## $ South Carolina : num [1:86] 0.0118 0.0236 0.0354 0.0474 0.0594 ... ## $ South Dakota : num [1:86] 0.0144 0.029 0.0433 0.0575 0.0714 ... ## $ Tennessee : num [1:86] 0.0123 0.0245 0.0367 0.049 0.0612 ... ## $ Texas : num [1:86] 0.0146 0.0292 0.0435 0.0578 0.0724 ... ## $ Utah : num [1:86] 0.0171 0.034 0.051 0.0676 0.0846 ... ## $ Vermont : num [1:86] 0.0096 0.0195 0.0292 0.039 0.0489 ... ## $ Virginia : num [1:86] 0.0123 0.0245 0.0367 0.0489 0.0612 ... ## $ Washington : num [1:86] 0.0124 0.0248 0.0373 0.0497 0.0623 ... ## $ West Virginia : num [1:86] 0.0108 0.0219 0.0331 0.0443 0.0554 ... ## $ Wisconsin : num [1:86] 0.0116 0.0232 0.0349 0.0467 0.0586 ... ## $ Wyoming : num [1:86] 0.0132 0.0262 0.0392 0.0523 0.0655 ...
We can see that we create a list of vectors for each state.
Now let’s use
map_dfc():
# calculate median age for each state
age_stats <- age_stats %>%
select(-AGE) %>%
map_dfc(cumsum) %>%
map_dfc(function(x) x/x[nrow(age_stats)]) %>%
mutate(AGE = age_stats$AGE) %>%
select(AGE, everything())
glimpse(age_stats)
## Rows: 86 ## Columns: 52 ## $ AGE <dbl> 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 1… ## $ Alabama <dbl> 0.01217929, 0.02428807, 0.03623586, 0.04835742,… ## $ Alaska <dbl> 0.01524840, 0.03030168, 0.04521944, 0.05979438,… ## $ Arizona <dbl> 0.01273885, 0.02549314, 0.03824081, 0.05101803,… ## $ Arkansas <dbl> 0.01292266, 0.02569476, 0.03836806, 0.05126652,… ## $ California <dbl> 0.01283122, 0.02561725, 0.03841722, 0.05118957,… ## $ Colorado <dbl> 0.01217217, 0.02440058, 0.03655841, 0.04888552,… ## $ Connecticut <dbl> 0.01013228, 0.02030490, 0.03056879, 0.04107142,… ## $ Delaware <dbl> 0.01163322, 0.02326537, 0.03477678, 0.04661654,… ## $ `District of Columbia` <dbl> 0.01457035, 0.02759032, 0.04071434, 0.05331661,… ## $ Florida <dbl> 0.01095628, 0.02192577, 0.03277035, 0.04371311,… ## $ Georgia <dbl> 0.01284902, 0.02562823, 0.03837038, 0.05144559,… ## $ Hawaii <dbl> 0.01275029, 0.02580557, 0.03896040, 0.05191402,… ## $ Idaho <dbl> 0.01390752, 0.02744397, 0.04125812, 0.05494312,… ## $ Illinois <dbl> 0.01233645, 0.02447549, 0.03656625, 0.04878623,… ## $ Indiana <dbl> 0.01264485, 0.02525899, 0.03798841, 0.05071888,… ## $ Iowa <dbl> 0.01267030, 0.02538518, 0.03802630, 0.05063696,… ## $ Kansas <dbl> 0.01334268, 0.02681406, 0.04039729, 0.05398568,… ## $ Kentucky <dbl> 0.01257763, 0.02509149, 0.03759246, 0.05002287,… ## $ Louisiana <dbl> 0.01366990, 0.02722998, 0.04047237, 0.05358525,… ## $ Maine <dbl> 0.00949098, 0.01906255, 0.02881486, 0.03860408,… ## $ Maryland <dbl> 0.01233418, 0.02449538, 0.03670807, 0.04897542,… ## $ Massachusetts <dbl> 0.01064012, 0.02127391, 0.03194332, 0.04267411,… ## $ Michigan <dbl> 0.01142560, 0.02298902, 0.03448267, 0.04603339,… ## $ Minnesota <dbl> 0.01271147, 0.02555680, 0.03831076, 0.05107128,… ## $ Mississippi <dbl> 0.01277126, 0.02554419, 0.03822836, 0.05115873,… ## $ Missouri <dbl> 0.01232034, 0.02475678, 0.03708173, 0.04938873,… ## $ Montana <dbl> 0.01225303, 0.02440298, 0.03639442, 0.04838391,… ## $ Nebraska <dbl> 0.01411096, 0.02814111, 0.04197798, 0.05572191,… ## $ Nevada <dbl> 0.01250131, 0.02483232, 0.03739746, 0.04979506,… ## $ `New Hampshire` <dbl> 0.009324624, 0.018656767, 0.028516676, 0.038139… ## $ `New Jersey` <dbl> 0.01153672, 0.02318527, 0.03489888, 0.04677209,… ## $ `New Mexico` <dbl> 0.01241485, 0.02496946, 0.03762395, 0.05035577,… ## $ `New York` <dbl> 0.01217189, 0.02407385, 0.03599014, 0.04784986,… ## $ `North Carolina` <dbl> 0.01200077, 0.02400872, 0.03589836, 0.04789096,… ## $ `North Dakota` <dbl> 0.01498293, 0.02957241, 0.04375387, 0.05755379,… ## $ Ohio <dbl> 0.01195482, 0.02398530, 0.03604946, 0.04804797,… ## $ Oklahoma <dbl> 0.01360302, 0.02720220, 0.04094481, 0.05455859,… ## $ Oregon <dbl> 0.01137846, 0.02290605, 0.03443489, 0.04598165,… ## $ Pennsylvania <dbl> 0.01105806, 0.02218984, 0.03332365, 0.04452838,… ## $ `Rhode Island` <dbl> 0.01029817, 0.02052152, 0.03085757, 0.04130916,… ## $ `South Carolina` <dbl> 0.01182747, 0.02361386, 0.03537920, 0.04736508,… ## $ `South Dakota` <dbl> 0.01441743, 0.02898356, 0.04325579, 0.05748470,… ## $ Tennessee <dbl> 0.01234189, 0.02454161, 0.03672630, 0.04900204,… ## $ Texas <dbl> 0.01460720, 0.02916560, 0.04354217, 0.05781300,… ## $ Utah <dbl> 0.01705642, 0.03397314, 0.05104564, 0.06758513,… ## $ Vermont <dbl> 0.009603574, 0.019524225, 0.029198261, 0.039019… ## $ Virginia <dbl> 0.01229627, 0.02446226, 0.03665131, 0.04888352,… ## $ Washington <dbl> 0.01241119, 0.02476856, 0.03725623, 0.04972054,… ## $ `West Virginia` <dbl> 0.01083507, 0.02189516, 0.03312702, 0.04427952,… ## $ Wisconsin <dbl> 0.01158854, 0.02323645, 0.03492551, 0.04674149,… ## $ Wyoming <dbl> 0.01320760, 0.02620022, 0.03920990, 0.05229295,…
Great, we have a tidy dataframe with a column for each state storing important census information for both ethnicity and age. Now onto the other datasets!
Subsubsection 3.11.2.2 Violent Crime
For crime, we have the following data:
crime
## # A tibble: 510 × 14 ## State Area ...3 Population `Violent\ncrime… `Murder and \nnonne… ## <chr> <chr> <chr> <chr> <dbl> <dbl> ## 1 ALABAMA Metropoli… <NA> 3708033 NA NA ## 2 <NA> <NA> Area ac… 0.97099999… 18122 283 ## 3 <NA> <NA> Estimat… 1 18500 287 ## 4 <NA> Cities ou… <NA> 522241 NA NA ## 5 <NA> <NA> Area ac… 0.97399999… 3178 32 ## 6 <NA> <NA> Estimat… 1 3240 33 ## 7 <NA> Nonmetrop… <NA> 628705 NA NA ## 8 <NA> <NA> Area ac… 0.99399999… 1205 28 ## 9 <NA> <NA> Estimat… 1 1212 28 ## 10 <NA> State Tot… <NA> 4858979 22952 348 ## # … with 500 more rows, and 8 more variables: Rape ## (revised ## definition)2 <dbl>, ## # Rape ## (legacy ## definition)3 <dbl>, Robbery <dbl>, Aggravated ## assault <dbl>, ## # Property ## crime <dbl>, Burglary <dbl>, Larceny- ## theft <dbl>, ## # Motor ## vehicle ## theft <dbl>
If we take a look at what information is stored in each column...
colnames(crime)
## [1] "State" ## [2] "Area" ## [3] "...3" ## [4] "Population" ## [5] "Violent\ncrime1" ## [6] "Murder and \nnonnegligent \nmanslaughter" ## [7] "Rape\n(revised\ndefinition)2" ## [8] "Rape\n(legacy\ndefinition)3" ## [9] "Robbery" ## [10] "Aggravated \nassault" ## [11] "Property \ncrime" ## [12] "Burglary" ## [13] "Larceny-\ntheft" ## [14] "Motor \nvehicle \ntheft"
you see that it’s kind of a mess and there’s a whole bunch of information in there that we’re not necessarily interested in for this analysis.
Because of the messy names here (we’ll clean them up in a bit), we’ll see the column index to select columns instead of the complicated names. Also, we print a specified row of violent crime to observe the
X__1 group we are looking for – Rate per 100,000 inhabitants (per the study.)
violentcrime <- crime %>%
select(c(1,3,5))
violentcrime
## # A tibble: 510 × 3 ## State ...3 `Violent\ncrime1` ## <chr> <chr> <dbl> ## 1 ALABAMA <NA> NA ## 2 <NA> Area actually reporting 18122 ## 3 <NA> Estimated total 18500 ## 4 <NA> <NA> NA ## 5 <NA> Area actually reporting 3178 ## 6 <NA> Estimated total 3240 ## 7 <NA> <NA> NA ## 8 <NA> Area actually reporting 1205 ## 9 <NA> Estimated total 1212 ## 10 <NA> <NA> 22952 ## # … with 500 more rows
Great, so we’re starting to home in on the data we’re interested in but we’re ultimately interested in Rate per 100,000 inhabitants, so we need get all rows where the second column is equal to
Rate per 100,000 inhabitants.
However, as we can see above, the value for State in these rows is
NA, so we need to fill() that value with the state name that is listed in a previous row. Then we can select the rows where the second column is Rate per 100,000 inhabitants. After that, we no longer need the second column, so we’ll remove it.
violentcrime <- violentcrime %>%
fill(State) %>%
filter(.[[2]] == "Rate per 100,000 inhabitants") %>%
rename( violent_crime = `Violent\ncrime1`) %>%
select(-`...3`)
violentcrime
## # A tibble: 52 × 2 ## State violent_crime ## <chr> <dbl> ## 1 ALABAMA 472. ## 2 ALASKA 730. ## 3 ARIZONA 410. ## 4 ARKANSAS 521. ## 5 CALIFORNIA 426. ## 6 COLORADO 321 ## 7 CONNECTICUT 218. ## 8 DELAWARE 499 ## 9 DISTRICT OF COLUMBIA4 1269. ## 10 FLORIDA 462. ## # … with 42 more rows
If we look closely at our data, we’ll notice that some of our state names have 6s at the end of them. This will cause problems later on.
violentcrime$State[20]
## [1] "MAINE6"
So, let’s clean that up now be removing those trailing numeric values and converting the names to lower case:
# lower case and remove numbers from State column
violentcrime <- violentcrime %>%
mutate(State = tolower(gsub('[0-9]+', '', State)))
violentcrime
## # A tibble: 52 × 2 ## State violent_crime ## <chr> <dbl> ## 1 alabama 472. ## 2 alaska 730. ## 3 arizona 410. ## 4 arkansas 521. ## 5 california 426. ## 6 colorado 321 ## 7 connecticut 218. ## 8 delaware 499 ## 9 district of columbia 1269. ## 10 florida 462. ## # … with 42 more rows
We’ve now got ourselves a tidy dataset with violent crime information that’s ready to be joined with our census_stats data!
# join with census data
firearms <- left_join(census_stats, violentcrime,
by = c("NAME" = "State"))
firearms
## # A tibble: 51 × 7 ## NAME white black hispanic male total_pop violent_crime ## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> ## 1 alabama 69.5 26.7 4.13 48.5 4850858 472. ## 2 alaska 66.5 3.67 6.82 52.4 737979 730. ## 3 arizona 83.5 4.80 30.9 49.7 6802262 410. ## 4 arkansas 79.6 15.7 7.18 49.1 2975626 521. ## 5 california 73.0 6.49 38.7 49.7 39032444 426. ## 6 colorado 87.6 4.47 21.3 50.3 5440445 321 ## 7 connecticut 80.9 11.6 15.3 48.8 3593862 218. ## 8 delaware 70.3 22.5 8.96 48.4 944107 499 ## 9 district of columbia 44.1 48.5 10.7 47.4 672736 1269. ## 10 florida 77.7 16.9 24.7 48.9 20268567 462. ## # … with 41 more rows
Subsubsection 3.11.2.3 Brady Scores
The study by AJPH groups the scores using 7 different categories. The study removed all weightings of the different laws in favor of a “1 law 1 point” system, since the weightings were “somewhat arbitrary.”
For the purpose of practice and simplification we will just keep the first line of “total state points” from the Brady Scorecard as they are given. This will be where our analysis differs from the study. We need to transform the data frame so that we have a column of state names and a column of the corresponding total scores.
brady
## # A tibble: 116 × 54 ## `States can recei… `Category Point… `Sub Category P… Points AL AK AR ## <chr> <dbl> <dbl> <dbl> <chr> <chr> <chr> ## 1 TOTAL STATE POINTS NA NA NA -18 -30 -24 ## 2 CATEGORY 1: KEEP… 50 NA NA <NA> <NA> <NA> ## 3 BACKGROUND CHECKS… NA 25 NA AL AK AR ## 4 Background Checks… NA NA 25 <NA> <NA> <NA> ## 5 Background Checks… NA NA 20 <NA> <NA> <NA> ## 6 Background Checks… NA NA 5 <NA> <NA> <NA> ## 7 Verifiy Legal Pur… NA NA 20 <NA> <NA> <NA> ## 8 TOTAL NA NA NA 0 0 0 ## 9 <NA> NA NA NA <NA> <NA> <NA> ## 10 OTHER LAWS TO STO… NA 12 NA AL AK AR ## # … with 106 more rows, and 47 more variables: AZ <chr>, CA <chr>, CO <chr>, ## # CT <chr>, DE <chr>, FL <chr>, GA <chr>, HI <chr>, ID <chr>, IL <chr>, ## # IN <chr>, IA <chr>, KS <chr>, KY <chr>, LA <chr>, MA <chr>, MD <chr>, ## # ME <chr>, MI <chr>, MN <chr>, MO <chr>, MT <chr>, MS <chr>, NC <chr>, ## # ND <chr>, NE <chr>, NH <chr>, NJ <chr>, NM <chr>, NV <chr>, NY <chr>, ## # OK <chr>, OH <chr>, OR <chr>, PA <chr>, RI <chr>, SC <chr>, SD <chr>, ## # TN <chr>, TX <chr>, UT <chr>, VA <chr>, VT <chr>, WA <chr>, WI <chr>, …
This dataset includes a lot of information, but we’re interested in the brady scores for each state. These are stored in the row where the first column is equal to "TOTAL STATE POINTS," so we
filter() to only include that row. We then want to only receive the scores for each state, and not the information in the first few columns, so we specify that using select(). With the information we’re interested in, we then take the data from wide to long using pivot_longer(), renaming the columns as we go. Finally, we specify that the information in the brady_scores column is numeric, not a character.
brady <- brady %>%
rename(Law = `States can receive a maximum of 100 points`) %>%
filter(Law == "TOTAL STATE POINTS") %>%
select((ncol(brady) - 49):ncol(brady)) %>%
pivot_longer(everything(),
names_to = "state",
values_to = "brady_scores") %>%
mutate_at("brady_scores", as.numeric)
brady
## # A tibble: 50 × 2 ## state brady_scores ## <chr> <dbl> ## 1 AL -18 ## 2 AK -30 ## 3 AR -24 ## 4 AZ -39 ## 5 CA 76 ## 6 CO 22 ## 7 CT 73 ## 8 DE 41 ## 9 FL -20.5 ## 10 GA -18 ## # … with 40 more rows
Only problem now is that we have the two letter state code, rather than the full state name we’ve been joining on so far here. We can, however, use the
state datasets we used in the first case study here!
brady <- brady %>%
left_join(rename(state_data, state = abb),
by = "state") %>%
select(Location, brady_scores) %>%
rename(state = Location) %>%
mutate(state = tolower(state))
brady
## # A tibble: 50 × 2 ## state brady_scores ## <chr> <dbl> ## 1 alabama -18 ## 2 alaska -30 ## 3 arkansas -24 ## 4 arizona -39 ## 5 california 76 ## 6 colorado 22 ## 7 connecticut 73 ## 8 delaware 41 ## 9 florida -20.5 ## 10 georgia -18 ## # … with 40 more rows
Now, it’s time to join this information into our growing dataframe
firearms:
firearms <- left_join(firearms, brady, by = c("NAME" = "state"))
firearms
## # A tibble: 51 × 8 ## NAME white black hispanic male total_pop violent_crime brady_scores ## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> ## 1 alabama 69.5 26.7 4.13 48.5 4850858 472. -18 ## 2 alaska 66.5 3.67 6.82 52.4 737979 730. -30 ## 3 arizona 83.5 4.80 30.9 49.7 6802262 410. -39 ## 4 arkansas 79.6 15.7 7.18 49.1 2975626 521. -24 ## 5 california 73.0 6.49 38.7 49.7 39032444 426. 76 ## 6 colorado 87.6 4.47 21.3 50.3 5440445 321 22 ## 7 connecticut 80.9 11.6 15.3 48.8 3593862 218. 73 ## 8 delaware 70.3 22.5 8.96 48.4 944107 499 41 ## 9 district of columbia 44.1 48.5 10.7 47.4 672736 1269. NA ## 10 florida 77.7 16.9 24.7 48.9 20268567 462. -20.5 ## # … with 41 more rows
Subsubsection 3.11.2.4 The Counted Fatal Shootings
We’re making progress, but we have a ways to go still! Let’s get working on incorporating data from The Counted.
As a reminder, we have a datasets here with data from 2015:
counted15
## # A tibble: 1,146 × 6 ## gender raceethnicity state classification lawenforcementage… armed ## <chr> <chr> <chr> <chr> <chr> <chr> ## 1 Male Black GA Death in custo… Chatham County Sh… No ## 2 Male White OR Gunshot Washington County… Firea… ## 3 Male White HI Struck by vehi… Kauai Police Depa… No ## 4 Male Hispanic/Latino KS Gunshot Wichita Police De… No ## 5 Male Asian/Pacific Islander WA Gunshot Mason County Sher… Firea… ## 6 Male White CA Gunshot San Francisco Pol… Non-l… ## 7 Male Hispanic/Latino AZ Gunshot Chandler Police D… Firea… ## 8 Male Hispanic/Latino CO Gunshot Evans Police Depa… Other ## 9 Male White CA Gunshot Stockton Police D… Knife ## 10 Male Black CA Taser Los Angeles Count… No ## # … with 1,136 more rows
The data from each year are in a similar format with each row representing a different individual and the columns being consistent between the two datasets.
Because of this consistent format, we can combine these two datasets using
bind_rows(). By specifying id = "dataset", a column called dataset will store which dataset each row came from originally. We can then use mutate() and ifelse() to conditionally specify the year -- 2015 or 2016 -- from which the data originated. We’ll also be sure to change the two letter state abbreviation to the lower case state name, to allow for each merging.
counted15 <- counted15 %>%
mutate(state = tolower(state.name[match(state, state.abb)]))
At this point, we have a lot of information at the individual level, but we’d like to summarize this at the state level by ethnicity, gender, and armed status. The researchers “calculated descriptive statistics for the proportion of victims that were male, armed, and non-White,” so we’ll do the same. We can accomplish this using
dplyr. The tally() function will be particularly helpful here to count the number of observations in each group. We’re calculating this for each state as well as calculating the annualized rate per 1,000,000 residents. This utilizes the total_pop column from the census_stats data frame we used earlier.
# get overall stats
counted_stats <- counted15 %>%
group_by(state) %>%
filter(classification == "Gunshot") %>%
tally() %>%
rename("gunshot_tally" = "n")
# get summary for subset of population
gunshot_filtered <- counted15 %>%
group_by(state) %>%
filter(classification == "Gunshot",raceethnicity != "white", armed != "No", gender == "Male") %>%
tally() %>%
rename("gunshot_filtered" = "n")
# join data together
counted_stats <- left_join(counted_stats, gunshot_filtered, by = "state") %>%
mutate(total_pop = census_stats$total_pop[match(state, census_stats$NAME)],
gunshot_rate = (gunshot_tally/total_pop)*1000000/2) %>%
select(-total_pop)
counted_stats
## # A tibble: 50 × 4 ## state gunshot_tally gunshot_filtered gunshot_rate ## <chr> <int> <int> <dbl> ## 1 alabama 18 15 1.86 ## 2 alaska 4 4 2.71 ## 3 arizona 43 37 3.16 ## 4 arkansas 5 4 0.840 ## 5 california 196 150 2.51 ## 6 colorado 29 27 2.67 ## 7 connecticut 2 2 0.278 ## 8 delaware 3 2 1.59 ## 9 district of columbia 5 4 3.72 ## 10 florida 64 54 1.58 ## # … with 40 more rows
Time to merge this into the data frame we’ve been compiling:
firearms <- left_join(firearms, counted_stats, by = c("NAME" = "state"))
Subsubsection 3.11.2.5 Unemployment Data
Let’s recall the table we scraped from the web, which is currently storing our unemployment data:
unemployment
## # A tibble: 54 × 3 ## State `2015rate` Rank ## <chr> <chr> <chr> ## 1 "United States" "5.3" "" ## 2 "" "" "" ## 3 "North Dakota" "2.8" "1" ## 4 "Nebraska" "3.0" "2" ## 5 "South Dakota" "3.1" "3" ## 6 "New Hampshire" "3.4" "4" ## 7 "Hawaii" "3.6" "5" ## 8 "Utah" "3.6" "5" ## 9 "Vermont" "3.6" "5" ## 10 "Minnesota" "3.7" "8" ## # … with 44 more rows
Let’s first rename the columns to clean things up. You’ll note that there are more rows in this data frame (due to an empty row, the United States, and a note being in this dataset); however, when we
left_merge() in just a second these will disappear, so we can ignore them for now.
unemployment <- unemployment %>%
rename("state" = "State",
"unemployment_rate" = "2015rate",
"unemployment_rank" = "Rank") %>%
mutate(state = tolower(state)) %>%
arrange(state)
unemployment
## # A tibble: 54 × 3 ## state unemployment_rate unemployment_rank ## <chr> <chr> <chr> ## 1 "" "" "" ## 2 "alabama" "6.1" "42" ## 3 "alaska" "6.5" "47" ## 4 "arizona" "6.1" "42" ## 5 "arkansas" "5.0" "24" ## 6 "california" "6.2" "44" ## 7 "colorado" "3.9" "10" ## 8 "connecticut" "5.7" "35" ## 9 "delaware" "4.9" "22" ## 10 "district of columbia" "6.9" "51" ## # … with 44 more rows
Let’s do that join now. Let’s add unemployment information to our growing data frame!
firearms <- left_join(firearms, unemployment, by = c("NAME" = "state"))
If we take a look at the data we now have in our growing data frame, using
glimpse(), we see that type is correct for most of our variables except unemployment_rate and unemployment_rank. This is due to that "Note" and empty ("") row in the unemployment dataset. So, let’s be sure to get that variable to a numeric now as it should be:
glimpse(firearms)
## Rows: 51 ## Columns: 13 ## $ NAME <chr> "alabama", "alaska", "arizona", "arkansas", "califor… ## $ white <dbl> 69.50197, 66.51368, 83.52295, 79.57623, 72.97555, 87… ## $ black <dbl> 26.7459489, 3.6679906, 4.7978011, 15.6634268, 6.4910… ## $ hispanic <dbl> 4.129434, 6.821197, 30.873010, 7.180439, 38.727129, … ## $ male <dbl> 48.46650, 52.36978, 49.71595, 49.12855, 49.67815, 50… ## $ total_pop <dbl> 4850858, 737979, 6802262, 2975626, 39032444, 5440445… ## $ violent_crime <dbl> 472.4, 730.2, 410.2, 521.3, 426.3, 321.0, 218.5, 499… ## $ brady_scores <dbl> -18.0, -30.0, -39.0, -24.0, 76.0, 22.0, 73.0, 41.0, … ## $ gunshot_tally <int> 18, 4, 43, 5, 196, 29, 2, 3, 5, 64, 29, 2, 7, 22, 19… ## $ gunshot_filtered <int> 15, 4, 37, 4, 150, 27, 2, 2, 4, 54, 25, 2, 7, 21, 15… ## $ gunshot_rate <dbl> 1.8553419, 2.7101042, 3.1607133, 0.8401593, 2.510731… ## $ unemployment_rate <chr> "6.1", "6.5", "6.1", "5.0", "6.2", "3.9", "5.7", "4.… ## $ unemployment_rank <chr> "42", "47", "42", "24", "44", "10", "35", "22", "51"…
# convert type for unemployment columns
firearms <- firearms %>%
mutate_at("unemployment_rate", as.numeric) %>%
mutate_at("unemployment_rank", as.integer)
Subsubsection 3.11.2.6 Population Density: 2015
Population density for 2015 can be calculated from the Census data in combination with the land area data we’ve read in. This is calculated (rather than simply imported) because accurate data for state population in 2015 was not available in a downloadable format nor was it easy to scrape.
From the census data, we can obtain total population counts:
totalPop <- census %>%
filter(ORIGIN == 0, SEX == 0 ) %>%
group_by(NAME) %>%
summarize(total = sum(POPESTIMATE2015)) %>%
mutate(NAME = tolower(NAME))
totalPop
## # A tibble: 51 × 2 ## NAME total ## <chr> <dbl> ## 1 alabama 4850858 ## 2 alaska 737979 ## 3 arizona 6802262 ## 4 arkansas 2975626 ## 5 california 39032444 ## 6 colorado 5440445 ## 7 connecticut 3593862 ## 8 delaware 944107 ## 9 district of columbia 672736 ## 10 florida 20268567 ## # … with 41 more rows
Then, we select
LND110210D by looking at the land table and comparing values on other sites (such as the census or Wikipedia) to find the correct column. This column corresponds to land area in square miles. We’ll convert all state names to lower case for easy merging with our growing data frame in a few steps.
landSqMi <- land %>%
select(Areaname, land_area = LND110210D) %>%
mutate(Areaname = tolower(Areaname))
landSqMi
## # A tibble: 3,198 × 2 ## Areaname land_area ## <chr> <dbl> ## 1 united states 3531905. ## 2 alabama 50645. ## 3 autauga, al 594. ## 4 baldwin, al 1590. ## 5 barbour, al 885. ## 6 bibb, al 623. ## 7 blount, al 645. ## 8 bullock, al 623. ## 9 butler, al 777. ## 10 calhoun, al 606. ## # … with 3,188 more rows
Since
landSqMi gives us area for each town in addition to the states, we will want merge on the state names to obtain only the area for each state, removing the city- and nation-level data. Also, because "district of columbia" appears twice, we’ll use the distinct() function to only include on entry for "district of columbia"
We can then calculate density and remove the
total and land_area columns to only keep state name and density for each state:
popdensity <- left_join(totalPop, landSqMi, by=c("NAME" = "Areaname")) %>%
distinct() %>%
mutate(density = total/land_area) %>%
select(-c(total, land_area))
popdensity
## # A tibble: 51 × 2 ## NAME density ## <chr> <dbl> ## 1 alabama 95.8 ## 2 alaska 1.29 ## 3 arizona 59.9 ## 4 arkansas 57.2 ## 5 california 251. ## 6 colorado 52.5 ## 7 connecticut 742. ## 8 delaware 485. ## 9 district of columbia 11019. ## 10 florida 378. ## # … with 41 more rows
This can now be joined with our growing data frame:
firearms <- left_join(firearms, popdensity, by="NAME")
Subsubsection 3.11.2.7 Firearm Ownership
Last but not least, we calculate firearm ownership as a percent of firearm suicides to all suicides.
ownership_df <- as_tibble(list("NAME" = tolower(suicide_all$State),
"ownership" = suicide_firearm$Deaths/suicide_all$Deaths*100))
ownership_df
## # A tibble: 51 × 2 ## NAME ownership ## <chr> <dbl> ## 1 alabama 70.1 ## 2 alaska 59.9 ## 3 arizona 57.4 ## 4 arkansas 59.5 ## 5 california 37.3 ## 6 colorado 51.0 ## 7 connecticut 27.4 ## 8 delaware 49.8 ## 9 florida 52.0 ## 10 georgia 62.1 ## # … with 41 more rows
This can now be joined onto our tidy data frame:
firearms <- left_join(firearms, ownership_df, by="NAME")
And, with that, we’ve wrangled and tidied all these datasets into a single data frame. This can now be used for visualization and analysis!
Let’s save our new tidy data for case study #2.
save(firearms, file = here::here("data", "tidy_data", "case_study_2_tidy.rda"))
