Linear municipal model

With a linear model, we will find that the best predictors are household_median_income and northern_nj, along with some combination of the racial demographics. We get better results when the the race predictors are transformed – we use \(log(1+t)\) here, though \(1/(1+t)\) also works well.

Below, we construct a linear model with r_white, r_black, r_asian and r_other_races. We use the tidymodels idiom for building up our models in terms of recipes. Our model fit to the training data is fit_lm.

predictors_lm <-
  c(
    "household_median_income",
    "northern_nj",
    "r_white",
    "r_black",
    "r_asian",
    "r_other_races"
  )

mod_lm <- linear_reg()

recipe_lm <- recipe(incident_rate_per_officer ~ ., data = ct_train) %>% 
  step_rm(-incident_rate_per_officer,         # Outcome
          -officer_years,                     # Case weights
          -all_of(predictors_lm)) %>%         # Predictors
  step_log(starts_with("r_"), offset = 1)

wf_lm <- workflow() %>%
  add_model(mod_lm) %>%
  add_recipe(recipe_lm) %>%
  add_case_weights(officer_years)

fit_lm <- wf_lm %>%
  fit(data=ct_train)

All of the predictors are significant in this model, but it is worth noting that the adjusted \(R^2\) only reaches \(.325\). (We will not do much better than this in any of our models.)

summary(fit_lm %>% extract_fit_engine())$adj.r.squared

## [1] 0.3252659

tlm <- fit_lm %>% tidy()

The largest contributions are coming from household_median_income and northern_nj: There are lower rates of use of force in North Jersey and in higher income towns. In this model, the difference between North Jersey and South Jersey is about the same as a difference of $72,000 in household median income. Interpreting the coefficients of the racial demographic predictors is less straightforward. We will return to them in a subsequent post.

A ‘linear-plus’ model incorporating the Jersey Shore and officer count

The shore_town flag does not appear in our linear model, and adding it does not improve the model. However, recalling the interaction we saw in Part I between shore_town and officer_count – as well as our suspicion about the relationship between officer_count and incident_rate_per_officer – we may be tempted to add the interaction shore_town * officer_count to our model. Indeed, we do see an improvement, though we get our best model if we remove officer_count. We call this the ‘linear-plus’ model.

predictors_lm_plus <- c(predictors_lm, "shore_town", "shore_officer_count")
recipe_lm_plus <- recipe(incident_rate_per_officer ~ ., data = ct_train) %>% 
  step_mutate(shore_officer_count = shore_town * officer_count) %>%
  step_rm(-incident_rate_per_officer,         # Outcome
          -officer_years,                     # Case weights
          -all_of(predictors_lm_plus)) %>%    # Predictors
  step_log(starts_with("r_"), offset = 1)

wf_lm_plus <- wf_lm %>%
  update_recipe(recipe_lm_plus)

fit_lm_plus <- wf_lm_plus %>%
  fit(data = ct_train)

summary(fit_lm_plus %>% extract_fit_engine())$adj.r.squared

## [1] 0.3475749

tlm <- fit_lm_plus %>% tidy()

In the linear-plus model, lm_plus, the interaction is captured by shore_officer_count, which is just the officer_count when shore_town is TRUE, and is 0 when shore_town is false. So the lm_plus model has dependence on the officer count for shore towns only.

This model has a marginally better \(R^2 = .35\), but at the cost of adding a predictor, officer_count, that we did not want to add.

Cross validation

Let’s check the robustness of our models by using cross validation.

#### CV of linear models
set.seed(37)
folds <- vfold_cv(ct_train, v = 10, repeats = 5)
metrics_lm <- wf_lm %>%
  fit_resamples(folds) %>% 
  collect_metrics()
metrics_lm_plus <- wf_lm_plus %>%
  fit_resamples(folds) %>% 
  collect_metrics()

mt <- lst(metrics_lm, metrics_lm_plus) %>%
  bind_rows(.id = "model") %>%
  mutate(model = str_sub(model, 9)) %>%
  select(-.config)

We observe that on cross validated data, the linear model slightly beats the linear-plus model on both RMSE (.424 vs .436) and \(R^2\) (.254 vs .233). The advantage the linear-plus model showed before we did cross validation was probably a bit of overfitting. For now, we will take lm as our best model.

Problems with the linear municipal models

One obvious problem with the linear models is that they may predict negative incident rates. In fact, we do have one negative prediction in the training data.

fit_lm %>%
  extract_fit_engine() %>%
  augment(ct_train) %>% 
  ggplot(aes(x = .fitted, y = incident_rate_per_officer)) +
  geom_point(aes(size = officer_count, alpha = partial_year)) +
  geom_abline(color = "red") +
  ylim(0, 2.5) +
  theme(legend.position = "bottom") +
  labs(title = "Linear model fit",
       subtitle = "(Training data only)")

A second problem is that the model may be under-predicting the higher rates, though it is hard to say from this graph because of theheteroskedasticity:

par(mfrow = c(1, 2))
plot(fit_lm %>% extract_fit_engine(), which = c(1, 3))

Log-linear municipal model

We will try to address the linear model problems with a log-linear model, \(\log y = \alpha + \sum\beta_ix_i + \epsilon\).

# Log-linear model
mod_llm <- linear_reg() %>%
  set_engine("glm", family = Gamma("log"))

We can reuse our linear model recipe, except that we will have to filter out the rows with no incidents. We note that we can expect this to cause a bias away from zero for small departments.

recipe_llm <- recipe_lm %>% 
  step_filter(incident_rate_per_officer > 0)
wf_llm <- wf_lm %>%
  update_model(mod_llm) %>%
  update_recipe(recipe_llm)
fit_llm <- wf_llm %>%
  fit(ct_train)

We do not have negative predictions or heteroskedasticity with the log-linear model.

fit_llm %>%
  extract_fit_engine() %>%
  augment(ct_train %>% filter(incident_rate_per_officer > 0),
          type.predict = "response") %>% 
  ggplot(aes(x = .fitted, y = incident_rate_per_officer)) +
  geom_point(aes(size = officer_count, alpha = partial_year)) +
  geom_abline(color = "red") +
  ylim(0, 2.5) + xlim(0, 1.5) +
  theme(legend.position = "bottom") +
  labs(title = "Log-Linear model fit",
       subtitle = "(Training data only)")

par(mfrow = c(1, 2))
plot(fit_llm %>% extract_fit_engine(), which=c(1, 3))

Let’s also try a log-linear-plus variant, similar to linear-plus:

recipe_llm_plus <- recipe_lm_plus %>% 
  step_filter(incident_rate_per_officer > 0)
wf_llm_plus <- wf_llm %>%
  update_recipe(recipe_llm_plus)
fit_llm_plus <- wf_llm_plus %>% 
  fit(ct_train)

Cross validation

Let’s get cross validation metrics for the two log-linear models and compare to the linear models.

metrics_llm <- wf_llm %>%
  fit_resamples(folds) %>% 
  collect_metrics()
metrics_llm_plus <- wf_llm_plus %>%
  fit_resamples(folds) %>% 
  collect_metrics()
mt <- lst(metrics_lm, metrics_llm, metrics_lm_plus, metrics_llm_plus) %>%
  bind_rows(.id = "model") %>%
  mutate(model = str_sub(model, 9)) %>%
  select(-.config)

The differences here are very slight, but the log-linear model llm performs at least as well as the linear model lm, while adding some desirable features (no negative predictions, homoskedasticity). The lm_plus and llm_plus are not adding anything that would justify the compromise we made in constructing them (adding shore_officer_count).

So our best model so far is the log-linear model llm.

ML municipal model

Before turning to the agency predictors, let’s see if machine learning algorithms can turn up any patterns that our linear models are missing.

We are not interested in replacing our interpretable linear model with an opaque ML model that performs only marginally better, so we will not spend a huge amount of time here tuning. We will use the xgboost package to build our model.²

We will begin by constructing a recipe for xgboost. We let xgboost use all the municipal predictors except r_pacific_islander and r_native_american, which are usually very close to zero.

counts <- ct_train %>%
  count(r_native_american > 0.01, r_pacific_islander > 0.005)

We will do the same log transform of the r_ fields as we did in our linear models. Note we have to convert the logical columns (shore_town and northern_nj) to integers.

municipal_predictors <- ct %>% select(population:northern_nj) %>% names()
recipe_xgb <- recipe(incident_rate_per_officer ~ ., data = ct_train) %>%
  step_rm(-incident_rate_per_officer, -officer_years,
          -all_of(municipal_predictors)) %>%
  step_rm(r_pacific_islander, r_native_american) %>%
  step_log(starts_with("r_"), offset = 1) %>%
  step_integer(northern_nj, shore_town)

Here is the best tune we found for this recipe:

# XGBoost model with tuned parameters
mod_xgb <- boost_tree(
  mtry = .48,
  trees = 600,            # nrounds
  min_n = 24,             # min_child_weight
  tree_depth = 2,         # max_depth
  learn_rate = .0185,     # eta
  loss_reduction = .001,  # gamma
  sample_size = .715,
  stop_iter = 50          # early_stopping_rounds
) %>% 
  set_engine("xgboost",
             alpha = .1,
             lambda = 17500,
             nthread = 12,
             counts = FALSE) %>%
  set_mode("regression")  
# XGBoost workflow
wf_xgb <- workflow() %>%
  add_model(mod_xgb) %>%
  add_recipe(recipe_xgb) %>%
  add_case_weights(officer_years)

Comparison with log-linear model

If we look at the cross validation metrics and compare the XGBoost model to our existing models, we see that we have a small but noticeable improvement in both \(R^2\) and RMSE.

metrics_xgb <- wf_xgb %>%
  fit_resamples(folds) %>% 
  collect_metrics()
mt <- lst(metrics_lm, 
    metrics_llm, 
    metrics_lm_plus, 
    metrics_llm_plus,
    metrics_xgb) %>%
  bind_rows(.id = "model") %>%
  mutate(model = str_sub(model, 9)) %>%
  select(-.config)

Next, let’s compare the shape of the predictions generated by the log-linear and ML models:

fit_xgb <- wf_xgb %>%
  fit(ct_train)
ct_train_aug <- bind_cols(
  ct_train,
  fit_llm %>% predict(ct_train) %>% rename(llm.pred = .pred),
  fit_xgb %>% predict(ct_train) %>% rename(xgb.pred = .pred)
)
ct_train_aug %>%
  pivot_longer(llm.pred:xgb.pred, names_to = "model", values_to = "pred") %>%
  ggplot(aes(x = incident_rate_per_officer, y = pred)) +
  geom_point(aes(size = officer_count, alpha = partial_year)) +
  geom_abline() +
  xlim(0,2.1) +
  facet_wrap("model") +
  labs(title = "Log-linear and XGBoost predictions by incident rate",
        subtitle = "(Training data only)")

ct_train_aug %>%
  ggplot(aes(y = llm.pred, x = xgb.pred)) +
  geom_point(aes(size = officer_count, alpha = partial_year)) +
  geom_abline() +
  xlim(.3, 1) +
  labs(title = "Log-linear predictions by XGBoost predictions",
        subtitle = "(Training data only)")

Here we see that that the ML model is far more conservative than the log linear model, never making predictions outside of the band \([.35,1]\), though the two models are pretty tightly correlated.

Next, let’s find out what the XGBoost model is seeing that the linear models are missing, by first taking a look at variable importance.

library(vip)
library(gridExtra)
grid.arrange(
  vip(fit_llm %>% extract_fit_engine()) + 
    labs(title = "Log-linear model vip"),
  vip(fit_xgb %>% extract_fit_engine()) + 
    labs(title = "XGBoost model vip"),
  ncol = 2
)

The ML model didn’t find much use for the variables omitted from our linear models, but it does rate r_asian more highly. Let’s see if we can find out where fit_xgb is making better predictions than fit_llm by breaking the training set down by the top three variables, r_asian, household_median_income, and northern_nj:

tbl <- ct_train_aug %>%
  mutate(region = if_else(northern_nj, "NORTH", "SOUTH"),
         income = if_else(household_median_income < 89000, "LOW", "HIGH"),
         r_asian = case_when(
           r_asian < .023 ~ "LOW",
           r_asian < .064 ~ "MEDIUM",
           TRUE ~ "HIGH"
         )) %>%
  group_by(r_asian, region, income) %>%
  summarize(count = n(),
            rate.mean = weighted.mean(incident_rate_per_officer,
                                      as.double(officer_years)),
            llm.mean = weighted.mean(abs(llm.pred), 
                                     as.double(officer_years)),
            xgb.mean = weighted.mean(abs(xgb.pred), 
                                     as.double(officer_years)),
            llm.mae = weighted.mean(abs(llm.pred-incident_rate_per_officer), 
                                    as.double(officer_years)),
            xgb.mae = weighted.mean(abs(xgb.pred-incident_rate_per_officer), 
                                    as.double(officer_years)),
            .groups="drop"
  ) %>%
  arrange(xgb.mae - llm.mae)

labeller = as_labeller(
  c("LOW" = "Low Asian density", "MEDIUM" = "Medium Asian density", 
    "HIGH" = "High Asian density",
    "NORTH" = "North Jersey", "SOUTH" = "South Jersey")
)
ct_train_aug %>%
  mutate(region = if_else(northern_nj, "NORTH", "SOUTH"),
         r_asian = factor(case_when(
           r_asian < .023 ~ "LOW",
           r_asian < .064 ~ "MEDIUM",
           TRUE ~ "HIGH"
         )),
         r_asian = fct_relevel(r_asian, "LOW", "MEDIUM", "HIGH"),
         xgb_v_llm =
           abs(xgb.pred - incident_rate_per_officer) -
           abs(llm.pred - incident_rate_per_officer)) %>%
  ggplot(aes(x = household_median_income, y = incident_rate_per_officer,
             color = xgb_v_llm)) +
  geom_point(aes(size = officer_count, alpha = partial_year)) +
  geom_smooth(aes(weight = as.double(officer_years)), method = "lm", se = FALSE,
              color = "grey") +
  annotate("rect", xmin = 0, xmax = 250000, ymin = .37, ymax = 1, 
           alpha = .1, fill = "green") +
  scale_color_gradient2(low = "red", mid = "white", high = "blue",
                        limits = c(-.3, .3), oob = squish) +
  scale_x_continuous(labels = label_number(suffix = "K", scale = .001)) +
  ylim(0, 2.2) +
  facet_grid(rows = vars(region), cols = vars(r_asian), labeller = labeller) +
  theme(legend.position = "none") +
  labs(title = "Model performance by Asian population density, region and income",
       subtitle = "Red = ML model better, Blue = Log-linear model better")

(The green bar represents the prediction range for the ML model)

General conclusions that we can draw from this table and graph:

• Household income becomes less important as the Asian population density increases.
• The log-linear model performs best relative to the ML model when the Asian population density is low.
• The log-linear model does better with high incident rates in South Jersey and low rates in North Jersey.
• Conversely, the ML model does better with low incident rates rates in South Jersey and high incident rates in North Jersey.

The South Jersey, low-income, low-Asian-density sector is particularly interesting:

• It has the highest mean rate of use of force per officer (1.12).
• It has the highest mean prediction for both ML (1.008) and log-linear models (.912).
• It is where the log-linear model (.345 MAE) has its best performance relative to ML model (.387 MAE).
• It has the highest MAE for the ML model in an absolute sense.
• It contains 6 of the 10 cases where the ML model outperforms the log-linear model.
• It contains 6 of the 10 cases where the log-linear model outperforms the ML model.

We could plausibly gain a marginal improvement in the linear models by adding some interaction terms between these three variables, but we will not pursue this here.

Linear agency models

Let’s start very simply with a linear model. First we set up the table with the residuals added.

ct_train_agcy <-
  bind_cols(
    ct_train,
    wf_llm %>% fit(ct_train) %>% predict(ct_train) %>% rename(llm.pred = .pred),
  ) %>%
  mutate(llm.resid = incident_rate_per_officer - llm.pred)

Now we set up a linear model of the residual in terms of the agency predictors.

mod_agcy_lm <- linear_reg()
agency_predictors <- ct_train_agcy %>%
  select(officer_mean_age:officer_r_race_na) %>% names()
recipe_agcy_lm <- recipe(llm.resid ~ ., data = ct_train_agcy) %>% 
  step_rm(-llm.resid,                         # Outcome
          -officer_years,                     # Case weights
          -all_of(agency_predictors))         # Predictors

wf_agcy_lm <- workflow() %>%
  add_model(mod_agcy_lm) %>%
  add_recipe(recipe_agcy_lm) %>%
  add_case_weights(officer_years)

fit_agcy_lm <- wf_agcy_lm %>%
  fit(ct_train_agcy)

tbl <- fit_agcy_lm %>% tidy()

The model is not promising. It has \(p=.423\) and adjusted \(R^2=.0003\). Log-transforming the officer_r_ fields, as we did for the r_ fields in the municipal models, does not turn out to be helpful. If we iteratively remove the worst predictor and refit, we arrive at our best linear model, which has a single predictor, officer_r_asian.

recipe_agcy_lm_asian <- recipe(llm.resid ~ officer_r_asian + officer_years,
                               data = ct_train_agcy)
wf_agcy_lm_asian <- wf_agcy_lm %>%
  update_recipe(recipe_agcy_lm_asian)
fit_agcy_lm_asian <- wf_agcy_lm_asian %>%
  fit(ct_train_agcy)
tbl <- fit_agcy_lm_asian %>% tidy()

This is not a very strong model, but on cross validation it looks worse, as it is beaten by the empty model (RMSE .418 vs .417 for the empty model). Visually, it is hard to see how the model could be much improved by non-linear terms without overfitting.

ct_train_agcy %>%
  ggplot(aes(x = officer_r_asian, y = llm.resid)) +
  geom_point(aes(size = officer_count, alpha = partial_year)) +
  geom_smooth(aes(weight = as.double(officer_years)), color = "lightyellow") +
  scale_x_log10() +
  theme(legend.position = "bottom") +
  ylim(-1, 1.8) + 
  labs(title = "Log Linear residuals by officer_r_asian",
       subtitle = "(Training data only)")

Let’s now see if we can find an ML algorithm that will turn up anything that we missed.

ML agency models

We will use XGBoost again. Here is our model and best tune:

recipe_agcy_xgb <-
  recipe(llm.resid ~ ., data = ct_train_agcy) %>%
  step_rm(GEOID:partial_year, incident_rate_est, population:llm.pred)

mod_agcy_xgb <- boost_tree(
  trees = 3000,           # nrounds
  min_n = 14,             # min_child_weight
  tree_depth = 11,        # max_depth
  learn_rate = .008,      # eta
  loss_reduction = 44,    # gamma
  sample_size = 1,
  stop_iter = 50          # early_stopping_rounds
) %>% 
  set_engine("xgboost",
             alpha = .003,
             lambda = .003,
             colsample_bytree = 1,
             nthread = 12,
             counts = FALSE) %>%
  set_mode("regression")

wf_agcy_xgb <- workflow() %>%
  add_model(mod_agcy_xgb) %>%
  add_recipe(recipe_agcy_xgb) %>%
  add_case_weights(officer_years)

Unfortunately this model barely outperforms the empty model (RMSE .416 vs .417 for the empty model, \(R^2\) = .03). The top variables in this model are officer_r_asian, which we examined when evaluating linear models, and officer_mean_age.

fit_agcy_xgb <- wf_agcy_xgb %>%
  fit(ct_train_agcy)
fit_agcy_xgb %>%
  extract_fit_engine() %>% vip()

Here is now fit_agcy_xgb uses officer_mean_age.

augment(fit_agcy_xgb, ct_train_agcy) %>%
  pivot_longer(c(".pred", "llm.resid"), 
               names_to = "residual_type", 
               values_to = "residual") %>%
  ggplot(aes(x = officer_mean_age, y = residual, color = residual_type)) +
  geom_point() +
  scale_color_hue(labels = c("Predicted", "Actual")) +
  ylim(-1.1, 2) +
  guides(color = guide_legend(title = "Residual Type")) +
  theme(legend.position = "bottom") +
  labs(title = "ML model predictions vs officer_mean_age")

Possibly there is something going on with younger agencies, but the prediction spike in the narrow band around age 36 appears to be an overfit driven by a few extreme values.

In any case, the ML model is barely distinguishable from the empty model.