A question about a new article on police shootings

I. Load Washington Post Data.

First, we have to load the necessary R packages and Washington Post data. I include at the beginning a function for returning state abbreviations from state names regardless of capitalization or punctuation. I also filter to exclude cases after 2016.

getState.abb <- function(x) {
    ## remove punctuation and twim white space, convert to lower case
    x <- tolower(trimws(gsub("[[:punct:]]", "", x)))
    x[which(x == "washington dc")] <- "district of columbia"
    
    states.abb <- c(state.abb, "DC", "DC", "PR")
    states.name <- c(state.name, "District of Columbia", "District of Columbia", 
        "Puerto Rico")
    
    m.abb <- tolower(states.abb)
    m.name <- tolower(gsub("[[:punct:]]", "", states.name))
    
    x[which(x %in% m.abb)] <- states.abb[match(x[which(x %in% m.abb)], m.abb)]
    x[which(x %in% m.name)] <- states.abb[match(x[which(x %in% m.name)], m.name)]
    
    return(x)
}


library(dplyr)
library(readr)

wapo <- read_csv("https://cdn.rawgit.com/washingtonpost/data-police-shootings/master/fatal-police-shootings-data.csv") %>% 
    dplyr::transmute(Race = race, Race = replace(Race, which(Race == "W"), "White"), 
        Race = replace(Race, which(Race == "B"), "Black"), Race = replace(Race, 
            which(Race == "H"), "Hispanic"), Race = replace(Race, which(Race == 
            "A"), "Asian"), Race = replace(Race, which(Race == "N"), "Native American"), 
        Race = replace(Race, which(Race == "O"), "Other"), year = as.numeric(format(date, 
            "%Y")), State = getState.abb(state)) %>% ## Sum across years for state for each race
dplyr::filter(year < 2017) %>% dplyr::group_by(State) %>% dplyr::summarise(killed.White = length(which(Race == 
    "White")), killed.Black = length(which(Race == "Black")), killed.Hispanic = length(which(Race == 
    "Hispanic")), killed.Asian = length(which(Race == "Asian")), killed.NativeAmerican = length(which(Race == 
    "Native American")), killed.Total = n())

knitr::kable(wapo)

II. Load Census Data

Next, we have to load the state-level population data by race. They use 2016 Census data, but we should get similar results from the 2015 Census Data, which would just as legitimate anyway since the Wasghinton Post data covers 2015 and 2016. The census data can be retrieved online:

State Characteristics Datasets: Annual State Resident Population Estimates for 5 Race Groups (5 Race Alone or in Combination Groups) by Age, Sex, and Hispanic Origin: April 1, 2010 to July 1, 2015

Using this data, you can also normalize by gender, age, and race and their cross-tabulations. Below I only examine race. The state names also have to be abbreviated to match the Washington Post data.

library(dplyr)
library(tidyr)
library(readr)

census_state <- read_csv("https://www2.census.gov/programs-surveys/popest/datasets/2010-2015/state/asrh/sc-est2015-alldata5.csv") %>% 
    ## filtering to include both males and females
dplyr::filter(SEX == 0) %>% ## selecting only 2015 data, renaming the columns
dplyr::transmute(State = NAME, Hisp = ORIGIN, Race = RACE, AGE = AGE, population = POPESTIMATE2015) %>% 
    ## Creating Race categories - 'Hispanics' are Hispanic Whites.  But the
## results are roughly the same if you include hispanic blacks, asians, etc.
## since there are so few of them.  I'm also combining Asians and PAcific
## Islanders since they are not distinguished in the WaPo data.
dplyr::mutate(State = getState.abb(State), Race = replace(Race, which(Race == 
    1 & Hisp == 1), "White"), Race = replace(Race, which(Race == 1 & Hisp == 
    2), "Hispanic"), Race = replace(Race, which(Race == 2 & Hisp == 0), "Black"), 
    Race = replace(Race, which(Race == 3 & Hisp == 0), "Native American"), Race = replace(Race, 
        which(Race == 4 & Hisp == 0), "Asian"), Race = replace(Race, which(Race == 
        5 & Hisp == 0), "Asian")) %>% dplyr::select(-c(Hisp)) %>% dplyr::filter(Race %in% 
    c("White", "Hispanic", "Black", "Asian", "Native American")) %>% dplyr::group_by(State, 
    Race) %>% ## Summing across all ages each race within each state
dplyr::summarise(population = sum(population, na.rm = T)) %>% tidyr::spread(Race, 
    population) %>% ## Creating a 'Total' population category
dplyr::mutate(Total = White + Hispanic + Black + Asian + `Native American`)

knitr::kable(census_state)

III. Combine Data and Calculate Rates per 100,000 Population

wapoRates <- dplyr::left_join(wapo, census_state) %>% dplyr::group_by(State) %>% 
    dplyr::transmute(Total = 1e+05 * killed.Total/Total, White = 1e+05 * killed.White/White, 
        Black = 1e+05 * killed.Black/Black, Hispanic = 1e+05 * killed.Hispanic/Hispanic, 
        `Native American` = 1e+05 * killed.NativeAmerican/`Native American`, 
        Asian = 1e+05 * killed.Asian/Asian)

knitr::kable(wapoRates)

IV. Calculate Means of State Rates for each Race

wapoRates_sum <- wapoRates %>% dplyr::ungroup() %>% dplyr::select(-State) %>% 
    dplyr::summarise_all(mean, na.rm = T) %>% tidyr::gather(Race, `Mean Killed Rate`, 
    1:6) %>% dplyr::arrange(desc(`Mean Killed Rate`))

knitr::kable(wapoRates_sum)

As shown above, the Black Rate is higher. Next, I’ll plot the results. Note that Washington D.c. Results are not included. I’ve calculated this elsewhere and it doesn’t much change the results.

V. Graph the results.

library(ggplot2)
library(ggsci)
library(ggthemes)
library(ggrepel)

## First change to long form
wapoRates_L <- wapoRates %>% dplyr::group_by(State) %>% tidyr::gather(Race, 
    KilledR, 2:7) %>% dplyr::ungroup() %>% dplyr::filter(!is.na(Race)) %>% dplyr::mutate(Race = factor(Race, 
    levels = c("Total", "White", "Black", "Hispanic", "Asian", "Native American"))) %>% 
    dplyr::group_by(Race) %>% dplyr::mutate(KilledR_mean = mean(KilledR, na.rm = T))

ggplot(data = wapoRates_L, aes(x = Race, y = KilledR, color = Race, fill = Race, 
    shape = Race)) + 
geom_boxplot(alpha = 0.5, position = "dodge", color = "black") + geom_dotplot(data = subset(wapoRates_L, 
    KilledR > 0), alpha = 0.5, position = "dodge", color = "black", binaxis = "y", 
    stackdir = "centerwhole", dotsize = 0.75, binwidth = 0.1) + xlab(label = "") + 
    ylab(label = "") + geom_label(aes(x = Race, y = KilledR_mean, group = Race, 
    label = paste0(round(KilledR_mean, 2))), size = 3, fill = "lightgray", alpha = 0.1, 
    color = "black", position = position_dodge(width = 0.9)) + ggsci::scale_fill_d3(palette = "category10") + 
    ggsci::scale_color_d3(palette = "category10") + theme_gdocs() + theme(text = element_text(size = 12, 
    face = "bold"), legend.position = c(0.12, 0.9), legend.text = element_text(size = 12, 
    face = "bold"), legend.background = element_rect(fill = alpha("gray", 0.2)), 
    axis.title.y = element_text(size = 12), axis.text = element_text(size = 10)) + 
    guides(size = FALSE, stat = FALSE, fill = FALSE, shape = FALSE, color = FALSE) + 
    ggtitle(label = "Race of Civilians Killed by Police per 100,000 Residents", 
        subtitle = "State-Level Data from  Washington Post (2015-2016) & US Census (2015).\nSource: johnbradford.github.io") + 
    coord_cartesian(ylim = c(0, 5)) + ggrepel::geom_text_repel(data = subset(wapoRates_L, 
    KilledR > 2.5), aes(x = Race, y = KilledR, label = State), color = "black", 
    size = 3)

These results contravene those of Shane, Lawton, and Swenson (2017) and are well known. I thought that maybe they had calculated their mean fatality rate using total population rather than race-specific population as the denominator. This seems doubtful. Re-calculating using total population as the denominator, the White rate is higher than the black rate, but the results are still far off. My results are 0.37 white killings per 100,000 people (of all races) and 0.14 killings of blacks per 100,000 people (of all races), rather than 0.74 and 0.63, respectively. Moreover, I could not replicate the ‘mean fatality base rate’ of 0.897 for Hispanics.

Was this a mistake on their part, or did I miss something?