Municipal agency diversity data

We start by getting a summary of the officer diversity data in njoagleod.

library(tidyverse)
library(njoagleod)
officer_info <- officer %>% group_by(agency_name, agency_county) %>%
  summarize(
    officer_count = n(),
    officer_r_male = sum(officer_gender == "Male", na.rm = TRUE) / n(),
    officer_mean_age = mean(officer_age, na.rm = TRUE),
    officer_r_race_na = sum(is.na(officer_race) / n()),
    officer_r_white = sum(officer_race == "White", na.rm = TRUE) / n(),
    officer_r_black = sum(officer_race == "Black", na.rm = TRUE) / n(),
    officer_r_asian = sum(officer_race == "Asian", na.rm = TRUE) / n(),
    officer_r_hispanic_or_latino =
      sum(officer_race == "Hispanic", na.rm = TRUE) / n(),
    officer_r_race_other = 1 - sum(c_across(
      officer_r_race_na:officer_r_hispanic_or_latino
    )),
    .groups = "drop"
  ) %>%
  left_join(agency, by = c("agency_name", "agency_county")) %>%
  filter(!is.na(municipality)) %>%
  left_join(municipality, by = c("agency_county"="county", "municipality")) %>%
  select(
    GEOID,
    officer_count,
    officer_mean_age,
    officer_r_male,
    officer_r_white:officer_r_hispanic_or_latino,
    officer_r_race_other,
    officer_r_race_na
  )

Note that the municipalites are identified here by US Census GEOID. We will repeat this pattern for each of the tables, joining them together at the end using the GEOID as the key.

Jersey Shore flag

We create a table with a shore_town flag, as in Part I. The flag is TRUE if the municipality has a border with the Atlantic Ocean.

library(tigris)
library(sf)
options(tigris_use_cache = TRUE)

counties <- municipality %>% pull(county) %>% unique() %>% sort()

# Get the map from tigris
nj_municipality_map <- counties %>%
  map_df(~ county_subdivisions("NJ", county = ., class = "sf"))

# Get the towns that border the Atlantic Ocean
jersey_shore_water <- nj_municipality_map %>%
  filter(ALAND == 0) %>%
  filter(GEOID != "3401100000")
jersey_shore_indices <- st_intersects(jersey_shore_water,
                                      nj_municipality_map) %>%
  unlist() %>% unique()

# Now build a table of shore towns
municipality_shore <- nj_municipality_map %>%
  as_tibble() %>%
  mutate(shore_town = FALSE)
municipality_shore$shore_town[jersey_shore_indices] = TRUE
municipality_shore <- municipality_shore %>%
  filter(ALAND != 0) %>%
  select(GEOID, shore_town)

Population estimates and race and ethnicity data

Next we gather the race and ethnicity data from the US Census. We are using the tidycensus package to get the data from the 2010 decennial census.

library(tidycensus)
# Get race/enthnicity ratios (2010 / get_decennial)
municipality_race_ethnicity_ratios <- counties %>%
  map_df(
    ~ get_decennial(
      geography = "county subdivision",
      state = "NJ",
      county = .,
      year = 2010,
      output = "wide",
      variables = c(
        population = "P001001",
        white = "P003002",
        black = "P003003",
        native_american = "P003004",
        asian = "P003005",
        pacific_islander = "P003006",
        other_races = "P003007",
        two_or_more_races = "P003008",
        non_hispanic_or_latino = "P004002",
        hispanic_or_latino = "P004003"
      )
    )
  ) %>%
  filter(population != 0) %>%
  # Convert totals to ratios
  mutate(across(white:hispanic_or_latino, ~ .x / population)) %>%
  rename_with( ~ paste0("r_", .), white:hispanic_or_latino) %>%
  select(GEOID, r_white:r_hispanic_or_latino)

Next we gather the the population and population density from the 2019 estimates, again using the tidycensus package.

# Get population and density  (2019 / get_estimates)
municipality_population <- counties %>%
  map_df(
    ~ get_estimates(
      geography = "county subdivision",
      state = "NJ",
      county = .,
      year = 2019,
      output = "wide",
      product = "population"
    )
  ) %>%
  select(GEOID, population = POP, density = DENSITY)

Household median income

Next, we get figures for the household median income from the 2015-2019 ACS.

municipality_income <- counties %>%
  map_df(
    ~ get_acs(
      geography = "county subdivision",
      state = "NJ",
      county = .,
      year = 2019,
      output = "wide",
      variables = c(household_median_income = "B19013_001")
    )
  ) %>%
  select(GEOID, household_median_income = household_median_incomeE) %>%
  filter(!is.na(household_median_income))

Use of force incident counts

Finally, we gather 2021 counts of use of force incidents from the njoaguof package. This task is somewhat complicated by the fact that we (probably) do not have the full year of data from every municipality. The earliest incidents in the dataset are from October 2020, but many individual agencies do not have incidents until 2021. For example, here are all the incidents reported by the Watchung PD, from October 2020 through May 2022.

library(njoaguof)
incident %>%
  filter(
      agency_county == "Somerset County" &
      agency_name == "Watchung PD"
  ) %>% 
  ggplot() + geom_histogram(aes(x=incident_date_1), binwidth=30) +
  scale_x_date(limits=lubridate::ymd(c("2020-10-01","2022-05-31"))) +
  xlab("Date") + ylab("Number of incidents") +
  labs(title=str_wrap("Watchung PD use of force incidents"),
       subtitle = str_wrap("Oct 2020 - May 2022"),
       caption = "Source: NJ OAG")

One suspects that there may have been use of force incidents between October 2020 and July 2021 that went unreported by the Watchung PD, simply because they had not yet started reporting such incidents at the time. It is not clear the earliest date at which the Watchung PD would have reported a use of force incident, but it may have been well into 2021.

For this reason, we will attempt to estimate when each department started reporting use of force incidents. We will make the simplifying assumption that for each department, the incidents represent a homogeneous Poisson process with an unknown starting point and continuing through the end of 2021. ¹ For agencies which reported their first incident in 2021, we estimate the rate \(\lambda\) of incidents by looking at the period after the first incident to the end of the year: \[\lambda_{est}=n/d\] where \(n\) is the number of incidents after the first one, and \(d\) is the number of days from the first incident to the end of 2021. For a Poisson process with rate \(\lambda\), the expected value of the time until the first incident is \(e^{-\lambda}/(1-e^{-\lambda})\). Thus we estimate the start of the reporting period \(t_0\) as \[{t_0}_{est} = t_1 - e^{-\lambda_{est}}/(1-e^{-\lambda_{est}})\] where \(t_1\) is the date of the first incident. ² Note that in order to use this estimate, there must be at least two incidents in 2021.

Let us now calculate \(\lambda_{est}\) and \({t_0}_{est}\) and the estimated proportion of the year 2021 this represents for every municipal agency who reported their first incident in 2021:

initial_incident_2021_estimates <- incident %>%
  # Find the first incident and the number of incidents after
  group_by(agency_county, agency_name) %>%
  summarize(
    first_incident = min(incident_date_1),
    incidents_after_first = sum(lubridate::year(incident_date_1) == 2021) - 1,
    .groups = "drop"
  ) %>%
  filter(lubridate::year(first_incident) == 2021) %>%
  mutate(
    # Estimate lambda
    lambda_est_period = lubridate::ymd("2022-01-01") - first_incident,
    lambda_est = incidents_after_first / as.double(lambda_est_period),
    # Estimate t_0
    t_0_est = case_when(
      incidents_after_first == 0 ~ as.Date(NA),
      TRUE ~ first_incident - exp(-lambda_est) / (1 - exp(-lambda_est)),
    ),
    # Estimate the partial year
    reporting_days_est = lubridate::ymd("2022-01-01") -
      pmax(t_0_est, lubridate::ymd("2021-01-01")),
    partial_year = as.double(reporting_days_est) /
      as.double(lubridate::ymd("2022-01-01") - lubridate::ymd("2021-01-01"))
  ) %>%
  # Filter out non-municipal agencies
  left_join(agency, by = c("agency_county", "agency_name")) %>%
  filter(!is.na(municipality)) %>%
  select(agency_county, municipality, lambda_est, t_0_est, partial_year)

There are 117 municipal agencies that reported their first incident in 2021. Of these, only 13 reported a single incident in 2021; we cannot estimate the the start of the reporting period for these. Of the remainder, we reckon that we have a full year of reporting for 49 and a partial year for 55:

initial_incident_2021_estimates %>% count(partial_year< 1)

## # A tibble: 3 × 2
##   `partial_year < 1`     n
##   <lgl>              <int>
## 1 FALSE                 49
## 2 TRUE                  55
## 3 NA                    13

Now we can build a table of all municipal agencies with an incident in either 2020 or 2021 and estimate the yearly rate of incidents for 2021.

• For municipalities with an incident in 2020, this rate is exactly the number of incidents in 2021.
• For municipalities that reported their first incident in 2021, the estimated rate is the number of incidents reported 2021 divided by partial_year.

incident_rate_2021 <- incident %>%
  # Get the 2021 incident count
  filter(lubridate::year(incident_date_1) <= 2021) %>%
  group_by(agency_county, agency_name) %>%
  summarize(incident_count = sum(lubridate::year(incident_date_1) == 2021),
            .groups = "drop") %>%
  # Filter out the non-municipal agencies
  left_join(agency, by = c("agency_county", "agency_name")) %>%
  filter(!is.na(municipality)) %>%
  select(agency_county, municipality, incident_count) %>%
  # Add a partial_year field, initialed to 1.0, and then update
  # with the initial_incident_2021_estimates
  mutate(partial_year = 1.0) %>%
  rows_update(
    select(
      initial_incident_2021_estimates,
      agency_county,
      municipality,
      partial_year
    ),
    by = c("agency_county", "municipality")
  ) %>%
  # Estimate the 2021 yearly rate
  mutate(incident_rate_est = incident_count / partial_year)

Finally, let’s switch to using the GEOID as the primary key.

municipality_incident_rate <- incident_rate_2021 %>%
  left_join(municipality, 
            by = c("agency_county" = "county", "municipality")) %>%
  select(GEOID, incident_count, partial_year, incident_rate_est)

Princeton Borough and Princeton Township

The former Princeton Township is present in the municipality_race_ethnicity_ratios table, since it is derived from the 2010 census, but absent from the other tables. We will handle this by updating municipality_race_ethnicity_ratios, combining the two Princetons, weighted by their 2010 populations.

# Combine the two Princetons.
#   Township:  Population 16265,  GEOID 3402160915
#   Borough:   Population 12307   GEOID 3402160900
combined_princeton <- municipality_race_ethnicity_ratios %>%
  filter(GEOID %in% c("3402160900", "3402160915")) %>%
  arrange(GEOID) %>%
  mutate(weight = c(12307, 16265)) %>%
  summarize(across(r_white:r_hispanic_or_latino,
                   ~ sum(.x * weight) / sum(weight))) %>%
  mutate(GEOID = "3402160900") %>%
  relocate(GEOID)

# Update the table with new ratios for Princeton Borough
municipality_race_ethnicity_ratios <- municipality_race_ethnicity_ratios %>% 
  rows_update(combined_princeton, by="GEOID")

US Census and Officer Diversity Race And Ethnicity Differences

Before combining the data, we should note here that the US Census and NJ OAG treat Hispanic and Latino origin differently. In the NJ OAG Law Enforcement Officer Diversity dataset, Hispanic and Latino origin is considered a race, so for officers we will have \[ r_{male} + r_{white} + \cdots + r_{non\_hispanic\_or\_latino} + \cdots + r_{two\_or\_more\_races} + r_{race\_na} = 1 \]

On the other hand, the US Census treats Hispanic and Latino origin as an ethnicity independent of race. Thus, for the general population, we will have \[ r_{white} + r_{black} + \cdots + r_{two\_or\_more\_races} = 1 \] and

\[ r_{non\_hispanic\_or\_latino} + r_{hispanic\_or\_latino} = 1 \]

The combined table

Now we can build the combined table.

combined_table <- officer_info %>% 
  inner_join(municipality_population, by="GEOID") %>%
  inner_join(municipality_race_ethnicity_ratios, by="GEOID") %>%
  inner_join(municipality_income, by="GEOID") %>%
  inner_join(municipality_shore, by="GEOID") %>%
  inner_join(municipality_incident_rate, by="GEOID")

This table has 407 rows, of which 395 have an estimate of the incident rate. These 395 rows cover 90% of the population of New Jersey.

sum(combined_table[!is.na(combined_table$incident_rate_est),]$population) /
  sum(municipality_population$population)

## [1] 0.9