*What is reproducible research?*

*What is reproducible research?*

The so-called “replication crisis” in psychology has left many researchers questioning whether the findings in most published research are true, or if they are simply the by-products of noisy measurements, small sample sizes, and sloppy analysis. Unfortunately, the methods and data of most studies are hidden and obfuscated from view. However, many science journals are now recognizing that transparent results and open data stand to significantly improve the quality of published research. Furthermore, if other researchers can replicate the findings of a study or build upon the author’s original work, the field of science is moved forward.

*Analyzing Police Stop Data*

*Analyzing Police Stop Data*

For this exercise, I’ve used Philadelphia Police Department’s own data on pedestrian and traffic stops. The data represent incident-level information on roughly 1.2 million police-citizen interactions between 2014 and 2016. Ten years ago, analyzing data of this magnitude would be nearly impossible. Today, more and more public agencies are releasing tons of freely open data. The exact data I used are available to download here. For this paper, I am using R for the data analysis, and Rmarkdown for the dissemination.

# Set working directory where data is located setwd("E:/") library(readr) library(lubridate) library(dplyr) # Read data and format dates stops <- data.frame(read_csv("PPD_Vehicle___Pedestrian_Investigations.csv"))

One advantage of reproducible research is that every step of the analysis is transparent. A lot of published research doesn’t discuss how data were recoded or operationalized. However, here you can clearly see how I recoded the data. I used lubridate to fix the dates so R can read them more easily. I used dplyr to select the variables I wanted for the analysis, and then recoded the race categories into more concise groupings. For my own convenience I recoded the outcome variables to “Yes-No” factors. Just a personal preference. What’s important here is that if you disagree, you can do the same analyses I did, but change various specifications. Strong results should be robust to specification changes.

stops$Date <- unlist(lapply(strsplit(stops$Date.Time, " "), "[",1)) stops$Date <- as.Date(stops$Date, format = "%m/%d/%Y") stops$yr_mon <- as.Date(paste0(format(stops$Date, "%Y-%m"),"-01")) stops$yr <- year(stops$Date) # Modeling pedestrian stops ped.stop <- stops %>% filter(Stop.Type == "pedestrian") %>% select(Police.District, Gender:Individual.Contraband, yr_mon, yr) ped.stop$Age[ped.stop$Age == 0] <- NA # Recode categories into factors # Using 'dplyr' makes it easy to see how I recode the data ped.stop$Race <- recode_factor(ped.stop$Race, `White - Non-Latino` = "White", `Black - Non-Latino` = "Black", `White - Latino` = "Latino", `Black - Latino` = "Latino", `Asian` = "Other", `American Indian` = "Other", `Unknown` = "Other") ped.stop$Gender <- recode_factor(ped.stop$Gender, `Male` = "Male", `Female` = "Female") ped.stop$Gender[ped.stop$Gender == "Unknown"] <- NA ped.stop$Individual.Frisked <- recode_factor(ped.stop$Individual.Frisked, `0` = "No", `1` = "Yes") ped.stop$Individual.Searched <- recode_factor(ped.stop$Individual.Searched, `0` = "No", `1` = "Yes") ped.stop$Individual.Arrested <- recode_factor(ped.stop$Individual.Arrested, `0` = "No", `1` = "Yes") ped.stop$Individual.Contraband<- recode_factor(ped.stop$Individual.Contraband, `0` = "No", `1` = "Yes")

ggplot is one of my favorite parts of using R. If you’re not aware, it’s a great graphical package developed by Hadley Wickham which makes visualizing data easy. First, we need to clean the data up a bit by binning the individual incidents into year-race groupings. This is really easy using the dplyr package. The ggthemes package makes visualizing data a lot easier, by providing pretty templates. They even have a fivethirtyeight theme!

# Data Visualization library(ggplot2) library(ggthemes) ######################################## # Percentage of pedestrian interactions by race, by year inter <- ped.stop %>% group_by(Race, yr) %>% summarise(count = n()) %>% group_by(yr) %>% mutate(freq = count / sum(count))

Here we see some information about the distribution of pedestrian stops by year, for each racial category. There is some evidence that blacks and hispanics are involved in pedesterian interactions at disproportinately high rates. For instance, blacks comprise roughly 70% of pedestrian interactions, while only making up about 44% of the population. These trends appear pretty stable from year-to-year, as well.

ggplot(inter) + geom_bar(aes(x = Race, y = freq, fill = Race), stat = "identity") + facet_wrap(~yr)+ labs(title ="Philadelphia Police Department \n Pedestrian Interactions", y = "Percentage of Pedestrian Interactions") + scale_y_continuous(labels = scales::percent) + scale_fill_ptol() + theme_minimal()

Next, we can look at individual outcomes in each of the pedesterian interactions. For the purposes of this study we’ll be focusing on four different areas: whether or not the individual was frisked, searched, arrested, or whether contraband was found. Before we get to the modeling stage, visualing the data is helpful. Let’s look at the distribution of frisks, searches, arrests, and contrand across racial categories by year

*Frisks*

frisk <- ped.stop %>% group_by(Race, yr, Individual.Frisked) %>% summarise(count = n()) %>% group_by(yr, Race) %>% mutate(freq = count / sum(count)) # Line ggplot(frisk[frisk$Individual.Frisked == "Yes",]) + geom_line(aes(x = Race, y = freq, group = yr), linetype = 2) + geom_point(aes(x = Race, y = freq)) + facet_wrap(~yr) + labs(title ="Pedestrian Interactions Resulting in a Frisk", y = "Percentage of Pedestrian Interactions") + scale_y_continuous(labels = scales::percent) + theme_minimal()

*Searches*

#Percentage of pedestrian interactions resulting in a search search <- ped.stop %>% group_by(Race, yr, Individual.Searched) %>% summarise(count = n()) %>% group_by(yr, Race) %>% mutate(freq = count / sum(count)) # Line ggplot(search[search$Individual.Searched == "Yes",]) + geom_line(aes(x = Race, y = freq, group = yr), linetype = 2) + geom_point(aes(x = Race, y = freq)) + facet_wrap(~yr) + labs(title ="Pedestrian Interactions Resulting in a Search", y = "Percentage of Pedestrian Interactions") + scale_y_continuous(labels = scales::percent) + theme_minimal()

*Arrests*

# Percentage of pedestrian interactions resulting in an arrest arrest <- ped.stop %>% group_by(Race, yr, Individual.Arrested) %>% summarise(count = n()) %>% group_by(yr, Race) %>% mutate(freq = count / sum(count)) # Line ggplot(arrest[arrest$Individual.Arrested == "Yes",]) + geom_line(aes(x = Race, y = freq, group = yr), linetype = "dashed") + geom_point(aes(x = Race, y = freq)) + facet_wrap(~yr) + labs(title ="Pedestrian Interactions Resulting in an Arrest", y = "Percentage of Pedestrian Interactions") + scale_y_continuous(labels = scales::percent) + scale_fill_ptol() + theme_minimal()

*Contraband*

#Percentage of pedestrian interactions where contraband is found cont <- ped.stop %>% group_by(Race, yr, Individual.Contraband) %>% summarise(count = n()) %>% group_by(yr, Race) %>% mutate(freq = count / sum(count)) # Line ggplot(cont[cont$Individual.Contraband == "Yes",]) + geom_line(aes(x = Race, y = freq, group = yr), linetype = "dashed") + geom_point(aes(x = Race, y = freq)) + facet_wrap(~yr) + labs(title ="Pedestrian Interactions Resulting in Contraband", y = "Percentage of Pedestrian Interactions") + scale_y_continuous(labels = scales::percent) + scale_fill_ptol() + theme_minimal()

What do these visuals suggest? As we see in the first plot, blacks and Hispanics are disproportionately more likely to be frisked during citizen-police interactions. When stopped, about 15 to 18% of interactions result in the citizen being frisked, while whites are frisked between 7 to 9%. While blacks are slightly more often searched or arrested, Hispanics are dramatically more often searched, arrested, and have contraband found. Between 12 to 14% of Hispanic pedestrian interactions result in a search or arrest.

*Modeling Probability of Finding Contraband*

What else is interesting is that while blacks are frisked and searched more often than whites, the percentage of cases where contraband is found is comparable or less than whites. This might suggest officers are frisking based on race, and not necessarily a sound belief that the individual has contraband. However, Hispanics are found with contraband almost twice as often as whites or blacks. This might be an interesting area to look into. However, its unclear whether Hispanics are more likely to be carrying contraband, or maybe officers are simply stopping younger people more often, who happen to be Hispanic. Let’s explore the probability of finding contraband using a logistic regression, which will allow us to control for each individual variable

First, we’ll fit individual models for each year. Our dependent variable is whether or not contraband was found during the interaction (Yes/No), so we’ll use a logistic regression. Our independent variables are the subjects gender, their race, and natural log of their age. There are some missing data in here, but for now we’ll assume they are missing completely at random (a pretty strong assumption, but that's a discussion for another time!), and carry on with the analysis.

######################################## #Predicting Contraband Recovered library(broom) # Models cont.14<- glm(Individual.Contraband ~ Gender + Race + log(Age), data = ped.stop[ped.stop$yr == 2014,], family = binomial) cont.15<- glm(Individual.Contraband ~ Gender + Race + log(Age), data = ped.stop[ped.stop$yr == 2015,], family = binomial) cont.16<- glm(Individual.Contraband ~ Gender + Race + log(Age), data = ped.stop[ped.stop$yr == 2016,], family = binomial)

2014 Incidents

summary(cont.14)

## ## Call: ## glm(formula = Individual.Contraband ~ Gender + Race + log(Age), ## family = binomial, data = ped.stop[ped.stop$yr == 2014, ]) ## ## Deviance Residuals: ## Min 1Q Median 3QMax ## -0.5141-0.2656-0.2433-0.2192 3.0569 ## ## Coefficients: ##Estimate Std. Error z value Pr(>|z|) ## (Intercept)-1.957090.11376 -17.203< 2e-16 *** ## GenderFemale -0.357690.03946-9.064< 2e-16 *** ## RaceBlack-0.033260.03384-0.9830.32564 ## RaceLatino0.807620.0423419.073< 2e-16 *** ## RaceOther-0.557790.14578-3.8260.00013 *** ## log(Age) -0.439060.03223 -13.623< 2e-16 *** ## --- ## Signif. codes:0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## (Dispersion parameter for binomial family taken to be 1) ## ## Null deviance: 57770on 204736degrees of freedom ## Residual deviance: 56900on 204731degrees of freedom ## (831 observations deleted due to missingness) ## AIC: 56912 ## ## Number of Fisher Scoring iterations: 6

2015 Incidents

summary(cont.15)

## ## Call: ## glm(formula = Individual.Contraband ~ Gender + Race + log(Age), ## family = binomial, data = ped.stop[ped.stop$yr == 2015, ]) ## ## Deviance Residuals: ## Min 1Q Median 3QMax ## -0.6099-0.2539-0.2289-0.2020 3.0138 ## ## Coefficients: ##Estimate Std. Error z value Pr(>|z|) ## (Intercept)-1.587440.11364 -13.969< 2e-16 *** ## GenderFemale -0.338830.03701-9.154< 2e-16 *** ## RaceBlack-0.247150.03160-7.822 5.19e-15 *** ## RaceLatino0.684680.0405316.895< 2e-16 *** ## RaceOther-0.514890.12541-4.106 4.03e-05 *** ## log(Age) -0.528840.03220 -16.426< 2e-16 *** ## --- ## Signif. codes:0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## (Dispersion parameter for binomial family taken to be 1) ## ## Null deviance: 60224on 231751degrees of freedom ## Residual deviance: 59189on 231746degrees of freedom ## (575 observations deleted due to missingness) ## AIC: 59201 ## ## Number of Fisher Scoring iterations: 6

2016 Incidents

summary(cont.16)

## ## Call: ## glm(formula = Individual.Contraband ~ Gender + Race + log(Age), ## family = binomial, data = ped.stop[ped.stop$yr == 2016, ]) ## ## Deviance Residuals: ## Min 1Q Median 3QMax ## -0.7139-0.2737-0.2370-0.2055 3.1143 ## ## Coefficients: ##Estimate Std. Error z value Pr(>|z|) ## (Intercept)-1.237130.17209-7.189 6.52e-13 *** ## GenderFemale -0.401430.05879-6.828 8.60e-12 *** ## RaceBlack 0.153910.05067 3.0370.00239 ** ## RaceLatino0.670220.0661110.137< 2e-16 *** ## RaceOther-0.352600.19602-1.7990.07205 . ## log(Age) -0.696200.04885 -14.252< 2e-16 *** ## --- ## Signif. codes:0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## (Dispersion parameter for binomial family taken to be 1) ## ## Null deviance: 26184on 97034degrees of freedom ## Residual deviance: 25786on 97029degrees of freedom ## (295 observations deleted due to missingness) ## AIC: 25798 ## ## Number of Fisher Scoring iterations: 6

*Closing Remarks*

Based on the models above, we can observe a few interesting patterns. Controlling for age and sex, Latinos are consistently more often found with contraband relative to whites - roughly twice as likely. Strikingly, blacks are just about equally likely to be found with contraband as whites. The above plot shows yearly estimates which hover around 0. What do these suggest about race and police interactions?

It appears that officers tend to stop, frisk, and search black individuals more often than whites. However, the probability of finding contraband is roughly equal or less than that of whites. The “hit rate”, so to speak, is quite low among blacks. This might indicate that officers are stopping many black individuals who are not carrying contraband. There is already evidence that officers in urban departments have engaged in racially biased stop, frisk, and searches. Additionally, the results here may suggest that Philadelphia police officers are more targeted against Latinos, and are more often searching individuals who are actually involved in criminal activity.

Do these data definitely suggest that Philadelphia Police are racially biased in pedestrian stops? It is difficult to say. While minorities are quite obviously more often stopped and searched by police relative to their proportion of the population, these data alone do not provide stone-cold evidence. In a way, this is why I am a proponent of open data and reproducible research. Rather than make decisions based a single study, science should move forward based on the current data and prior knowledge (something the Bayesians have known for a long time).I would encourage any skeptic, critic, or interested researcher to take these data and examine the results for themselves. Healthy disagreement is the heart of science.

As a side note, I think R and R markdown have a great opportunity for being used as a teaching tool. I like the idea of being able to teach complex concepts (ie: regression, data analytic techniques, spatial statistics) in a more visual manner, and being able to walk the class, point-by-point, through the steps. I have to thank Ashton Shortridge from MSU’s geography department for giving me the idea.