NJOAG Use of Force and Law Enforcement Officer Diversity data sets

The NJ Office of the Attorney General has started releasing data on police use of force and on law enforcement officer diversity, which I have re-packaged as R data packages NJOAGUOF and NJOAGLEOD.

I’d like to start digging into this data to see what it shows. I’ll start with the municipal agencies – in particular, agency sizes. ¹

install_github("tor-gu/njoaguof")
install_github("tor-gu/njoagleod")

Municipal agencies and census data

library(tidyverse)
library(tidycensus)
library(njoagleod)
# Get a vector of NJ counties
counties <- municipality %>% pull(county) %>% unique() %>% sort()
# Get the population of each municipality using tidycensus::get_estimates
municipality_pop <- counties %>%
  map_df(
    ~ get_estimates(
      geography = "county subdivision",
      state = "NJ",
      county = .,
      year = 2019,
      variables = "POP"
    ) %>%
      separate(
        NAME,
        sep = ", ",
        into = c("municipality", "county", "state")
      )
  ) %>%
  rename(population=value) %>%
  select(-variable, -state)

# Get the officer count for each agency
officer_count_pop <- officer %>%
  count(agency_county, agency_name, name = "officer_count") %>%
  left_join(agency,
            by = c("agency_county", "agency_name")) %>%
  filter(!is.na(municipality)) %>%
  left_join(municipality,
            by = c("agency_county" = "county", "municipality")) %>%
  select(GEOID, agency_name, officer_count) %>%
  left_join(municipality_pop, by = "GEOID")

# We cover nearly 93% of the population of NJ
sum(officer_count_pop$population) / sum(municipality_pop$population)

## [1] 0.93

# About 408 residents per municipal police officer
sum(officer_count_pop$population) / sum(officer_count_pop$officer_count)

## [1] 408

# Get the population for each municipality that is 'missed' and summarize by
# county
missed_county_pop <- municipality_pop %>%
  anti_join(officer_count_pop, by = "GEOID") %>%
  group_by(county) %>%
  summarize(missed_population = sum(population))
# Get the total population by county
county_pop <- municipality_pop %>%
  group_by(county) %>%
  summarize(population = sum(population))
# Combine this to get the county-level coverage -- the proportion of the 
# population of each county with a municipal PD for which we have the data
coverage <- county_pop %>%
  left_join(missed_county_pop, by = "county") %>%
  mutate(missed_population = replace_na(missed_population, 0)) %>%
  mutate(coverage = (population - missed_population) / population) %>%
  select(county, missed_population, coverage)
# Now add the county-level municipal PD officer count
officer_count_pop_by_county <- officer_count_pop %>%
  group_by(county) %>%
  summarize(officer_count = sum(officer_count),
            population = sum(population)) %>%
  mutate(resident_per_officer = population / officer_count) %>%
  left_join(coverage, by = "county")

# Weighted by population, there is no significant correlation between
# residents_per_officer and our coverage.
lm(resident_per_officer ~ coverage,
   officer_count_pop_by_county,
   weights = population) %>%
  summary() %>%
  magrittr::use_series("coefficients")

##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)      430        184   2.334    0.031
## coverage         -11        196  -0.054    0.957

# Plot number of residents per officer by coverage
outliers <- officer_count_pop_by_county %>%
  filter(
    population == min(population) | population == max(population) |
      coverage == min(coverage) | coverage == max(coverage) |
      resident_per_officer == min(resident_per_officer) |
      resident_per_officer == max(resident_per_officer)
  )
officer_count_pop_by_county %>%
  ggplot(aes(x = coverage, y = resident_per_officer)) +
  geom_point(aes(size = population)) +  
  ggrepel::geom_label_repel(data = outliers, aes(label = county)) +
  labs(
    title =
      str_wrap(
        "County population vs proportion of residents with single-municipality PDs"
      ),
    subtitle = "2019 Population Estimates, 2021 Officer Counts",
    caption = "Source: US Census, NJ OAG"
  )

Mapping municipalities with large police departments relative to their population

Now that we have our data set, let’s take a look at the municipalities with the largest police departments relative to their population:

top_fifteen <- officer_count_pop %>%
  mutate(residents_per_officer = population / officer_count) %>%
  slice_min(n = 15, order_by = residents_per_officer) %>%
  select(agency_name,
         county,
         residents_per_officer,
         officer_count,
         population)

In case the names of these relatively smaller municipalities don’t mean much to you, let’s take a look at them on a map. We’ll use the tigris package to get the shapefiles for the municipalities from the US census:

library(tigris)
options(tigris_use_cache = TRUE)
nj_municipality_map <- counties %>%
  map_df( ~ county_subdivisions("NJ", county = ., class = "sf"))

Next we join the officer count data to the map and plot.

# Add the residents_per_officer to the map
officer_count_pop <- officer_count_pop %>%
  mutate(residents_per_officer = population / officer_count)
map_with_values <-
  left_join(nj_municipality_map, officer_count_pop, by = "GEOID")

# Select the top values for labeling
labels <- map_with_values %>%
  slice_min(n = 15, order_by = residents_per_officer) %>%
  arrange(desc(INTPTLAT))

# Map the officer density
scale_breaks <- c(30, 40, 60, 120)
ggplot(map_with_values) +
  geom_sf(aes(fill = residents_per_officer, geometry = geometry, size = .2 )) +
  ggrepel::geom_label_repel(
    data = labels,
    mapping = aes(label = NAME, geometry = geometry),
    stat = "sf_coordinates",
    nudge_x = -.7,
    min.segment.length = 0
  ) +
  scale_size_identity() +
  scale_fill_gradient(
    na.value = "lightgrey",
    low = "white",
    high = "red",
    trans = "reciprocal",
    breaks = scale_breaks,
    name = "Residents Per Officer"
  ) +
  labs(
    title = str_wrap(
      "NJ municipalities with large police departments relative to population"
    ),
    subtitle = "2019 Population Estimates, 2021 Officer Counts",
    caption = "Source: US Census, NJ OAG"
  ) +
  theme(
    axis.ticks = element_blank(),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    panel.background = element_rect(fill = "lightblue"),
    panel.grid.major = element_line(color = "lightblue"),
    legend.position = "bottom"
  )

So, a map of the largest percapita municipal PDs in NJ is a map of the Jersey Shore, plus Hi-Nella.

Part-time residents and part-time officers

Isolating the shore towns

The NJ map includes several landless non-municipal county subdivisions off the coast, which we can identify by filtering on ALAND == 0.

# Select the landless regions and plot the GEOIDs
water <- nj_municipality_map %>% filter(ALAND == 0)
ggplot(nj_municipality_map) +
  geom_sf(aes(geometry = geometry)) +
  geom_sf_label(data = water,
                mapping = aes(label = GEOID, geometry = geometry)) +
  labs(title = str_wrap("Landless county subdivision GEOIDs in NJ"),
       caption = "Source: US Census") +
  theme(
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank()
  )

We can take a stab at defining the “Jersey Shore” as those municipalities that border any of these regions, excluding the Delaware Bay region 3401100000.

# Build the jersey shore map
jersey_shore_water <- water %>% filter(GEOID != "3401100000")
jersey_shore_indices <-
  sf::st_intersects(jersey_shore_water, nj_municipality_map) %>%
  unlist() %>% unique()
jersey_shore <-
  nj_municipality_map[jersey_shore_indices, ] %>% setdiff(water)

# Plot the jersey_shore towns
ggplot(jersey_shore) + geom_sf() +
  ggrepel::geom_label_repel(
    mapping = aes(label = NAME, geometry = geometry),
    stat = "sf_coordinates",
    nudge_x = -.7,
    min.segment.length = 0,
    size = 3
  ) +
  labs(title = str_wrap("52 \"Jersey Shore\" municipalities"),
       caption = "Source: US Census") +
  theme(
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank()
  )

Now we can add a flag, shore_town, to the officer_count_pop table:

officer_count_pop <- officer_count_pop %>% 
  mutate(shore_town=GEOID %in% jersey_shore$GEOID)

Comparing the Jersey Shore to the rest of the state

As expected, the small Jersey Shore towns (like Mantoloking and Deal) have noticeably larger police forces relative to their population than the rest of the state. But this distinction starts to break down with larger shore towns:

# Plot the officer count vs the population for small and medium towns,
# separating the shore towns from the rest of the state
officer_count_pop %>%
  mutate(town_size = factor(
    case_when(
      population < 5000 ~ "Small",
      population < 50000 ~ "Medium",
      TRUE ~ "Large"
    ),
    levels = c("Small", "Medium", "Large")
  )) %>%
  filter(town_size != "Large") %>%
  ggplot(aes(x = population, y = officer_count, color = shore_town)) +
  geom_point() +
  geom_smooth(se = FALSE, method="loess") +
  facet_wrap("town_size", scales = "free") +
  labs(title = "NJ Residents Per Municipal Law Enforcement Officer",
       subtitle = "Small and Medium Sized Towns, Jersey Shore vs Rest of State",
       caption = "Source: US Census, NJ OAG") +
  theme(legend.position = "bottom")

If we look at just the 26 small shore towns – population at most 5000 – the number of residents per officer is 111, compared to 415 for the rest of the state:

shore_vs_non_shore <- officer_count_pop %>%
  mutate(small_shore_town = (population <= 5000 & shore_town)) %>%
  group_by(small_shore_town) %>%
  summarise(officer_count = sum(officer_count),
            population = sum(population)) %>%
  mutate(res_per_officer = population / officer_count)

With this in mind, let’s redo our look at the municipalities with the largest police departments relative to population after excluding the small shore towns. Our new top fifteen is:

officer_count_pop_core <- officer_count_pop %>%
  filter(population > 5000 | shore_town == FALSE)

top_fifteen_core <- officer_count_pop_core %>% 
  slice_min(n = 15, order_by = residents_per_officer) %>%
  select(agency_name,
         county,
         residents_per_officer,
         officer_count,
         population)

Hi-Nella is the new leader, followed by Far Hills and Lower Alloways Creek. It is not clear if either Far Hills or Lower Alloways Creek is including part-time officers, like Hi-Nella does. The size of the Lower Alloways Creek PD may be related to the presence of the nuclear power plant. West Wildwood and Ocean Gate should probably be considered part of the Jersey Shore, though they don’t actually border the ocean.

Geographically, no clear pattern appears, except perhaps that South Jersey is somewhat over-represented:

# Add the residents_per_officer to the map, excluding shore towns
map_with_values <-
  left_join(nj_municipality_map, officer_count_pop_core, by = "GEOID")

# Map the officer density
labels <- map_with_values %>%
  slice_min(n = 15, order_by = residents_per_officer) %>%
  arrange(desc(INTPTLAT))

# Plot the map
scale_breaks <- c(90, 120, 180, 360)
ggplot(map_with_values) +
  geom_sf(aes(fill = residents_per_officer, geometry = geometry, size = .2 )) +
  ggrepel::geom_label_repel(
    data = labels,
    mapping = aes(label = NAME, geometry = geometry),
    stat = "sf_coordinates",
    nudge_x = -.7,
    min.segment.length = 0
  ) +
  scale_size_identity() +
  scale_fill_gradient(
    na.value = "lightgrey",
    low = "white",
    high = "red",
    trans = "reciprocal",
    breaks = scale_breaks,
    name = "Residents Per Officer"
  ) +
  labs(
    title = str_wrap(
      "NJ municipalities with large police departments relative to population, small shore towns excluded"
    ),
    subtitle = "2019 Population Estimates, 2021 Officer Counts",
    caption = "Source: US Census, NJ OAG"
  ) +
  theme(
    axis.ticks = element_blank(),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
    panel.background = element_rect(fill = "lightblue"),
    panel.grid.major = element_line(color = "lightblue"),
    legend.position = "bottom"
  )

Next

In subsequent posts, we dig in further, looking at:

• Municipality demographics (US Census)
• Department diversity (njoagleod package)
• Use of force (njoaguof package)