Introduction

On July 26, 2020, Propublica released a dataset on police discipline records.

I wanted to get a look at the data and do some exploratory analysis.

library(tidyverse)
## ── Attaching packages ──────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.2     ✓ purrr   0.3.4
## ✓ tibble  3.0.3     ✓ dplyr   1.0.0
## ✓ tidyr   1.1.0     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.5.0
## Warning: package 'ggplot2' was built under R version 3.6.2
## Warning: package 'tibble' was built under R version 3.6.2
## Warning: package 'tidyr' was built under R version 3.6.2
## Warning: package 'purrr' was built under R version 3.6.2
## Warning: package 'dplyr' was built under R version 3.6.2
## ── Conflicts ─────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(gt)
## Warning: package 'gt' was built under R version 3.6.2
library(here)
## here() starts at /Users/samportnow/sam-portnow-website
library(fs)
## Warning: package 'fs' was built under R version 3.6.2
library(skimr)
## Warning: package 'skimr' was built under R version 3.6.2
library(janitor)
## Warning: package 'janitor' was built under R version 3.6.2
## 
## Attaching package: 'janitor'
## The following objects are masked from 'package:stats':
## 
##     chisq.test, fisher.test

Read in Data

data = read_csv(here::here('content', 'post_data', 'allegations_20200726939.csv'))
## Parsed with column specification:
## cols(
##   .default = col_character(),
##   unique_mos_id = col_double(),
##   complaint_id = col_double(),
##   month_received = col_double(),
##   year_received = col_double(),
##   month_closed = col_double(),
##   year_closed = col_double(),
##   mos_age_incident = col_double(),
##   complainant_age_incident = col_double(),
##   precinct = col_double()
## )
## See spec(...) for full column specifications.

Explore Data

skim(data)

Variable type: character

Variable type: numeric

data %>% tabyl(mos_ethnicity) %>% adorn_pct_formatting() %>% gt()

It looks like the officer ethnicity matches the overall NYPD demogrphaic breakdown

Complainant Ethnicity

data %>% tabyl(complainant_ethnicity) %>% adorn_pct_formatting() %>% gt()

As for as the complaintant data is concerned, there isn’t an obvious population to compare to, but I think the closest is the arrest data because arrests can give some approximation of contact with the police.

Here, we see that Black people are overrepresented in their complaints compared to the percentage of arrestees that are Black. Also, Hispanic people are markedly underrepresented. I wonder if this is because of a large undocumented population in New York who may try to avoid contact with the government, or perhaps an issue with the missing ethnicity category.

Cross Tab Complaintant Ethnicity and Officer Ethnicity

data %>% tabyl(complainant_ethnicity, mos_ethnicity) %>%
  adorn_percentages("row") %>%
  adorn_pct_formatting(digits = 2) %>%
  adorn_ns() %>% gt()

The first thing I know is that there doesn’t seem to be much over or underrepresentation in the officer ethnicity relative to the demographic breakdown of the NYPD, with the exception of White, Native American, and Missing Complaitants and White officers. For these groups, which White officers are overrepresnted relative to the demographic breakdown of the NYPD. Actually, it appears that all ethnicity groups overrepresent their officers in complaints when the offending officer is of their same ethnicity, with the exception of Black complaitants. There is a large psychological literature that suggests people on on average display intergroup bias in favor of their own race. Perhaps when someone is mistreated by someone of their own race, it feels like a larger transgression.

Top Level Category of Complaint

data %>% tabyl(fado_type) %>% adorn_pct_formatting() %>% gt()

data %>% tabyl(allegation) %>% adorn_pct_formatting() %>% arrange(desc(n)) %>% gt()

Contact Reason

data %>% tabyl(contact_reason) %>% arrange(desc(n)) %>% adorn_pct_formatting() %>% gt()

data %>% tabyl(outcome_description) %>% arrange(desc(n)) %>% adorn_pct_formatting() %>% gt()

Disposition

data %>% tabyl(board_disposition) %>% adorn_pct_formatting() %>% gt()

So, now I want to know what is the best predictor of having a substantiated claim. I am going to use tidymodels to evaluate a lasso regression to determine the most important predictors. The largely majority of this code is borrowed from Julia Silge

First, I’ll code my outcome

library(stringr)
data = data %>%
  mutate(
    outcome = if_else(str_detect(board_disposition, 'Substantiated'), 1, 0)
  )

Next, I’ll select the columns I want.

library(tidymodels)
## ── Attaching packages ─────────────────────────────────────────────────────────────────────── tidymodels 0.1.1 ──
## ✓ broom     0.7.0      ✓ recipes   0.1.13
## ✓ dials     0.0.8      ✓ rsample   0.0.7 
## ✓ infer     0.5.3      ✓ tune      0.1.1 
## ✓ modeldata 0.0.2      ✓ workflows 0.1.2 
## ✓ parsnip   0.1.2      ✓ yardstick 0.0.7
## Warning: package 'broom' was built under R version 3.6.2
## Warning: package 'dials' was built under R version 3.6.2
## Warning: package 'scales' was built under R version 3.6.2
## Warning: package 'infer' was built under R version 3.6.2
## Warning: package 'modeldata' was built under R version 3.6.2
## Warning: package 'recipes' was built under R version 3.6.2
## Warning: package 'rsample' was built under R version 3.6.2
## Warning: package 'tune' was built under R version 3.6.2
## Warning: package 'yardstick' was built under R version 3.6.2
## ── Conflicts ────────────────────────────────────────────────────────────────────────── tidymodels_conflicts() ──
## x scales::discard() masks purrr::discard()
## x dplyr::filter()   masks stats::filter()
## x recipes::fixed()  masks stringr::fixed()
## x dplyr::lag()      masks stats::lag()
## x yardstick::spec() masks readr::spec()
## x recipes::step()   masks stats::step()
data = data %>% select(command_now, command_at_incident, year_received, year_closed, rank_abbrev_incident, rank_abbrev_now,
                       mos_ethnicity:contact_reason, outcome)

I want to see what predicts a founded complaint, and the allegation of the complaint seems too proximally related to the complaints, so I am going to drop those variables.

data = data %>%
  select(- c(fado_type, allegation))

Create a length of investigation variable

data$length_of_investigation = data$year_closed - data$year_received

Change years to factors

data$year_closed = factor(data$year_closed)
data$year_received = factor(data$year_received)

Recode character variables to factor

data = data %>% mutate_if(is.character, as.factor)

Recoding missing in factor variables as missing

data = data %>% mutate_if(is.factor, ~ fct_explicit_na(., 'Missing'))

Change outcome to a factor

data$outcome = as.factor(data$outcome)

Look at data again

skim(data)

Variable type: factor

Variable type: numeric

data = data %>% select(-complainant_age_incident)
data$precinct = factor(data$precinct)
data$precinct = fct_explicit_na(data$precinct, 'Missing')

skim(data)

Variable type: factor

Variable type: numeric

Finally, if any factor has an n less than 10, I went to recode that to ther

recode_to_other = function(.x){
  
  keeps = which(table(.x) > 10) %>% names()
  fct_other(.x, keep = keeps, other_level = 'other')
}

data = data %>% mutate_if(is.factor, recode_to_other)
## Warning: Unknown levels in `f`: other

## Warning: Unknown levels in `f`: other

## Warning: Unknown levels in `f`: other

## Warning: Unknown levels in `f`: other

## Warning: Unknown levels in `f`: other

Before fitting the model, I want to know the distribution of outcomes.

prop.table(table(data$outcome))
## 
##         0         1 
## 0.7511541 0.2488459

75% of complains are unfounded. So, if we fit a model without any predictors, we’ll already be at 75% accuracy. I confirm this below with a logistic regression with no predictors. The predict probability of an sustantiated claim is ~25%.

glm.mod = glm(outcome ~ 1, data = data, family = binomial(link = 'logit'))
predictions = predict(glm.mod, type = 'response')
table(predictions)
## predictions
## 0.248845854068537 
##             33358
data_split <- initial_split(data)
data_train <- training(data_split)
data_test <- testing(data_split)
data_rec <- recipe(outcome~ ., data = data_train) %>%
  step_zv(all_numeric(), -all_outcomes()) %>%
  step_dummy(all_predictors(), -all_numeric()) %>%
  step_normalize(all_numeric(), -all_outcomes())


data_prep <- data_rec %>%
  prep(strings_as_factors = FALSE)

summary(data_prep)
## # A tibble: 805 x 4
##    variable                type    role      source  
##    <chr>                   <chr>   <chr>     <chr>   
##  1 mos_age_incident        numeric predictor original
##  2 length_of_investigation numeric predictor original
##  3 outcome                 nominal outcome   original
##  4 command_now_X001.PCT    numeric predictor derived 
##  5 command_now_X005.DET    numeric predictor derived 
##  6 command_now_X005.PCT    numeric predictor derived 
##  7 command_now_X006.PCT    numeric predictor derived 
##  8 command_now_X007.PCT    numeric predictor derived 
##  9 command_now_X009.DET    numeric predictor derived 
## 10 command_now_X009.PCT    numeric predictor derived 
## # … with 795 more rows
lasso_spec <- logistic_reg(mode = 'classification', penalty = 0.1, mixture = 1) %>%
  set_engine("glmnet")

wf <- workflow() %>%
  add_recipe(data_rec)

lasso_fit <- wf %>%
  add_model(lasso_spec) %>%
  fit(data = data_train)

lasso_fit %>%
  pull_workflow_fit() %>%
  tidy()
## # A tibble: 30,881 x 5
##    term         step estimate lambda dev.ratio
##    <chr>       <dbl>    <dbl>  <dbl>     <dbl>
##  1 (Intercept)     1    -1.11 0.0295  9.99e-14
##  2 (Intercept)     2    -1.11 0.0269  1.68e- 3
##  3 (Intercept)     3    -1.11 0.0245  4.34e- 3
##  4 (Intercept)     4    -1.11 0.0224  6.95e- 3
##  5 (Intercept)     5    -1.11 0.0204  9.13e- 3
##  6 (Intercept)     6    -1.11 0.0186  1.15e- 2
##  7 (Intercept)     7    -1.11 0.0169  1.38e- 2
##  8 (Intercept)     8    -1.12 0.0154  1.60e- 2
##  9 (Intercept)     9    -1.12 0.0140  1.81e- 2
## 10 (Intercept)    10    -1.12 0.0128  2.01e- 2
## # … with 30,871 more rows
set.seed(1234)
data_boot <- bootstraps(data_train)

tune_spec <- logistic_reg(mode = 'classification', penalty = tune(), mixture = 1) %>%
  set_engine("glmnet")

lambda_grid <- grid_regular(penalty(), levels = 50)
doParallel::registerDoParallel()

set.seed(2020)
lasso_grid <- tune_grid(
  wf %>% add_model(tune_spec),
  resamples = data_boot,
  grid = lambda_grid
)

What results did we get?

lasso_grid %>%
  collect_metrics()
## # A tibble: 100 x 7
##     penalty .metric  .estimator  mean     n std_err .config
##       <dbl> <chr>    <chr>      <dbl> <int>   <dbl> <fct>  
##  1 1.00e-10 accuracy binary     0.739    14 0.00104 Model01
##  2 1.00e-10 roc_auc  binary     0.622    14 0.00102 Model01
##  3 1.60e-10 accuracy binary     0.739    14 0.00104 Model02
##  4 1.60e-10 roc_auc  binary     0.622    14 0.00102 Model02
##  5 2.56e-10 accuracy binary     0.739    14 0.00104 Model03
##  6 2.56e-10 roc_auc  binary     0.622    14 0.00102 Model03
##  7 4.09e-10 accuracy binary     0.739    14 0.00104 Model04
##  8 4.09e-10 roc_auc  binary     0.622    14 0.00102 Model04
##  9 6.55e-10 accuracy binary     0.739    14 0.00104 Model05
## 10 6.55e-10 roc_auc  binary     0.622    14 0.00102 Model05
## # … with 90 more rows
lasso_grid %>%
  collect_metrics() %>%
  ggplot(aes(penalty, mean, color = .metric)) +
  geom_errorbar(aes(
    ymin = mean - std_err,
    ymax = mean + std_err
  ),
  alpha = 0.5
  ) +
  geom_line(size = 1.5) +
  facet_wrap(~.metric, scales = "free", nrow = 2) +
  scale_x_log10() +
  theme(legend.position = "none")

highest_auc <- lasso_grid %>%
  select_best("roc_auc", maximize = FALSE)
## Warning: The `maximize` argument is no longer needed. This value was ignored.
final_lasso <- finalize_workflow(
  wf %>% add_model(tune_spec),
  highest_auc
)
library(vip)
## Warning: package 'vip' was built under R version 3.6.2
## 
## Attaching package: 'vip'
## The following object is masked from 'package:utils':
## 
##     vi
lasso_table = final_lasso %>%
  fit(data_train) %>%
  pull_workflow_fit() %>%
  vi(lambda = highest_auc$penalty) %>%
  mutate(
    Importance = abs(Importance),
    Variable = fct_reorder(Variable, Importance)
  )
lasso_table %>% mutate(Importance = round(Importance, 2)) %>% filter(Importance > .1) %>% arrange(desc(Importance)) %>% gt()

Probably the most interesting finding here is that the excecution of a search warrant predicts a lower likelihood fo having a substantiated claim, which makes sense.

Next, let’s look at how the model actually performed.

last_fit(
  final_lasso,
  data_split
) %>%
  collect_metrics()
## # A tibble: 2 x 3
##   .metric  .estimator .estimate
##   <chr>    <chr>          <dbl>
## 1 accuracy binary         0.748
## 2 roc_auc  binary         0.643

Unfortunately, the model did not perform well at all. We would just as accurate if we had went without any predictors at all.

Conclusion

This looks like a situation where a predictive model doesn’t tell us much of anything, which makes sense. I don’t really know much of anything with regards to substantiated claims of police misconduct, and I am probably missing lots of important variables that can describe the process. A predictive model that is applied to data without any idea on the data generation process is bound to be useless.

However, the descriptive statistics are very interesting here. In particular, I found it interesting that the following categories make up over 50% of all claims of police misconduct:

• PD suspected C/V of violation/crime – street
• Other
  
• PD suspected C/V of violation/crime – auto
  
• PD suspected C/V of violation/crime – bldg
  
• Moving violation