This document runs every regression reported in the paper: the five-specification Table 1 progression, the FE-2SLS model, the failed firefighter instrument, the Hausman and Durbin–Wu–Hausman tests, and the 2005–2024 robustness check. It is meant to be knit directly, code and output together, as the reproducible record behind the paper’s results.
PANEL_PATH <- "data_out/city_panel_v5.csv"
panel <- read_csv(PANEL_PATH, show_col_types = FALSE)
# Sanity check: this should be the cleaned, rescaled 2015-2024 panel.
stopifnot(
nrow(panel) == 2756,
n_distinct(panel$ori) == 297,
max(panel$poverty_rate, na.rm = TRUE) > 1 # confirms percentage-point scale
)
mv <- c("lcrmrte_total", "lpolpc", "lmedinc", "poverty_rate",
"unemp_rate", "pct_youngmale", "lpolpc_lag")
d <- panel %>% filter(if_all(all_of(mv), is.finite))
cat("Full panel: ", nrow(panel), " rows / ", n_distinct(panel$ori), " cities\n")## Full panel: 2756 rows / 297 citiescat("Modeling sample: ", nrow(d), " rows / ", n_distinct(d$ori), " cities\n")## Modeling sample: 1863 rows / 253 citiespd <- pdata.frame(d, index = c("ori", "year"))
ctrl <- "lmedinc + poverty_rate + unemp_rate + pct_youngmale"
clse <- function(m) coeftest(m, vcov = vcovHC(m, type = "sss", cluster = "group"))
f <- function(rhs, ...) plm(as.formula(paste("lcrmrte_total ~", rhs)), data = pd, ...)Five nested specifications of the same equation, adding controls and then fixed effects in sequence.
m1 <- f("lpolpc", model = "pooling")
m2 <- f(paste("lpolpc +", ctrl), model = "pooling")
m3 <- f(paste("lpolpc +", ctrl), model = "within", effect = "individual") # entity FE
m4 <- f(paste("lpolpc +", ctrl), model = "within", effect = "time") # time FE
m5 <- f(paste("lpolpc +", ctrl), model = "within", effect = "twoways") # baselinecat("===== Evolution of the police coefficient =====\n")## ===== Evolution of the police coefficient =====labs <- c("(1) police only", "(2) + controls", "(3) + entity FE",
"(4) + time FE", "(5) two-way FE [baseline]")
mods <- list(m1, m2, m3, m4, m5)
for (i in seq_along(mods)) {
ct <- clse(mods[[i]])
cat(sprintf("%-28s %+.3f (se %.3f) p = %.3f\n",
labs[i], ct["lpolpc", "Estimate"], ct["lpolpc", "Std. Error"],
ct["lpolpc", "Pr(>|t|)"]))
}## (1) police only +0.511 (se 0.083) p = 0.000
## (2) + controls +0.184 (se 0.099) p = 0.064
## (3) + entity FE -0.250 (se 0.244) p = 0.306
## (4) + time FE +0.250 (se 0.096) p = 0.009
## (5) two-way FE [baseline] -0.358 (se 0.234) p = 0.127se_list <- lapply(list(m1, m2, m3, m4, m5), function(m)
sqrt(diag(vcovHC(m, type = "sss", cluster = "group"))))
stargazer(m1, m2, m3, m4, m5,
type = "text",
se = se_list,
title = "Table 1: The Effect of Police on Crime, U.S. Cities 2015-2024",
dep.var.labels = "ln(Total crime rate per 100k)",
covariate.labels = c("ln(Police per capita)", "ln(Median income)",
"Poverty rate (pct. pts.)", "Unemployment rate",
"Share male 15-24 (pct. pts.)"),
column.labels = c("Police only", "+ Controls", "+ Entity FE",
"+ Time FE", "Two-way FE"),
add.lines = list(c("City FE", "No", "No", "Yes", "No", "Yes"),
c("Year FE", "No", "No", "No", "Yes", "Yes")),
omit.stat = c("f", "ser"),
notes = "Standard errors clustered by city. Column (5) is the baseline.",
digits = 3) Dependent variable:
------------------------------------------------------------------
ln(Total crime rate per 100k)
Police only + Controls + Entity FE + Time FE Two-way FE
(1) (2) (3) (4) (5) | ln(Police per capita) 0.511*** 0.184* -0.250 0.250*** -0.358 (0.083) (0.099) (0.244) (0.096) (0.234) |
| ln(Median income) -0.418*** -0.563*** -0.381** 0.027 (0.136) (0.114) (0.188) (0.321) |
| Poverty rate (pct. pts.) 0.019*** 0.025** 0.014* 0.018* (0.007) (0.010) (0.008) (0.010) |
| Unemployment rate -0.105*** -0.220*** -0.050* -0.090 (0.033) (0.054) (0.030) (0.074) |
| Share male 15-24 (pct. pts.) -0.031** 0.062*** -0.022* 0.029 (0.014) (0.024) (0.013) (0.023) |
| Constant 5.399*** 12.081*** (0.418) (1.790) |
City FE No No Yes No Yes
Year FE No No No Yes Yes
Observations 1,863 1,863 1,863 1,863 1,863
R2 0.056 0.107 0.074 0.093 0.009
Adjusted R2 0.056 0.104 -0.074 0.087 -0.155
===============================================================================================
Note: p<0.1; p<0.05; p<0.01 Standard
errors clustered by city. Column (5) is the baseline.
Every control loses significance once both fixed effects are included in column (5); the within-\(R^2\) falls sharply for the same reason, since the fixed effects absorb most of the cross-sectional variation by construction.
iv <- plm(as.formula(paste("lcrmrte_total ~ lpolpc +", ctrl,
"| lpolpc_lag +", ctrl)),
data = pd, model = "within", effect = "twoways")
cat("===== FE-2SLS (instrument = lagged police) =====\n")## ===== FE-2SLS (instrument = lagged police) =====print(clse(iv)["lpolpc", , drop = FALSE])## Estimate Std. Error t value Pr(>|t|)
## lpolpc -0.6237713 0.4161592 -1.498877 0.1341032fs <- plm(as.formula(paste("lpolpc ~ lpolpc_lag +", ctrl)),
data = pd, model = "within", effect = "twoways")
cat("First-stage strength (want > 10):\n")## First-stage strength (want > 10):print(linearHypothesis(fs, "lpolpc_lag = 0",
vcov = vcovHC(fs, type = "sss", cluster = "group")))##
## Linear hypothesis test:
## lpolpc_lag = 0
##
## Model 1: restricted model
## Model 2: lpolpc ~ lpolpc_lag + lmedinc + poverty_rate + unemp_rate + pct_youngmale
##
## Note: Coefficient covariance matrix supplied.
##
## Res.Df Df Chisq Pr(>Chisq)
## 1 1599
## 2 1598 1 135.51 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1d_aug <- d
d_aug$vhat <- residuals(fs)
pd_aug <- pdata.frame(d_aug, index = c("ori", "year"))
aug <- plm(as.formula(paste("lcrmrte_total ~ lpolpc +", ctrl, "+ vhat")),
data = pd_aug, model = "within", effect = "twoways")
ct_aug <- clse(aug)
cat("Coefficient on first-stage residual:\n")## Coefficient on first-stage residual:cat(sprintf(" %+.3f (se %.3f) p = %.4f\n",
ct_aug["vhat", "Estimate"], ct_aug["vhat", "Std. Error"],
ct_aug["vhat", "Pr(>|t|)"]))## +0.434 (se 0.519) p = 0.4031cat(if (ct_aug["vhat", "Pr(>|t|)"] < 0.10)
"-> Reject exogeneity: prefer the IV estimate.\n"
else
"-> Cannot reject exogeneity: OLS and IV are not statistically distinguishable.\n")## -> Cannot reject exogeneity: OLS and IV are not statistically distinguishable.Following Levitt (2002), we attempted a second instrument: metropolitan-area firefighter employment. Included here for completeness and transparency, not because we rely on it.
if ("lff_employment" %in% names(d) && sum(is.finite(d$lff_employment)) > 50) {
sub <- d %>% filter(is.finite(lff_employment))
pdf2 <- pdata.frame(sub, index = c("ori", "year"))
fs3 <- plm(as.formula(paste("lpolpc ~ lff_employment +", ctrl)),
data = pdf2, model = "within", effect = "twoways")
cat("First stage on firefighters (expected to fail):\n")
print(linearHypothesis(fs3, "lff_employment = 0",
vcov = vcovHC(fs3, type = "sss", cluster = "group")))
cat("\nGiven the weak first stage, we do not report a second-stage estimate.\n")
} else {
cat("Firefighter data not present in this panel version; skipping.\n")
}## First stage on firefighters (expected to fail):
##
## Linear hypothesis test:
## lff_employment = 0
##
## Model 1: restricted model
## Model 2: lpolpc ~ lff_employment + lmedinc + poverty_rate + unemp_rate +
## pct_youngmale
##
## Note: Coefficient covariance matrix supplied.
##
## Res.Df Df Chisq Pr(>Chisq)
## 1 1376
## 2 1375 1 0.5717 0.4496
##
## Given the weak first stage, we do not report a second-stage estimate.fe <- f(paste("lpolpc +", ctrl), model = "within", effect = "twoways")
re <- f(paste("lpolpc +", ctrl), model = "random", effect = "twoways")
cat("===== Hausman FE vs RE =====\n")## ===== Hausman FE vs RE =====print(phtest(fe, re))##
## Hausman Test
##
## data: as.formula(paste("lcrmrte_total ~", rhs))
## chisq = 15.87, df = 5, p-value = 0.007224
## alternative hypothesis: one model is inconsistentTests whether a longer panel changes the result, and whether the police-crime relationship differs before and after 2015.
PANEL_2005_PATH <- "data_out/city_panel_2005_v5.csv"
if (file.exists(PANEL_2005_PATH)) {
panel_2005 <- read_csv(PANEL_2005_PATH, show_col_types = FALSE)
d2 <- panel_2005 %>% filter(if_all(all_of(mv), is.finite))
cat("2005-2024 modeling sample:", nrow(d2), "rows /",
n_distinct(d2$ori), "cities\n")
} else {
cat("city_panel_2005_v5.csv not found; skipping the extended-panel check.\n")
d2 <- NULL
}## 2005-2024 modeling sample: 4304 rows / 257 citiespd2 <- pdata.frame(d2, index = c("ori", "year"))
m5_ext <- plm(as.formula(paste("lcrmrte_total ~ lpolpc +", ctrl)),
data = pd2, model = "within", effect = "twoways")
iv_ext <- plm(as.formula(paste("lcrmrte_total ~ lpolpc +", ctrl,
"| lpolpc_lag +", ctrl)),
data = pd2, model = "within", effect = "twoways")
cat("===== 2005-2024 baseline vs. 2015-2024 baseline =====\n")## ===== 2005-2024 baseline vs. 2015-2024 baseline =====ct5 <- clse(m5)
ct5e <- clse(m5_ext)
cat(sprintf("2015-2024 (n=%d): %+.3f (se %.3f) p=%.3f\n",
nrow(d), ct5["lpolpc","Estimate"], ct5["lpolpc","Std. Error"],
ct5["lpolpc","Pr(>|t|)"]))## 2015-2024 (n=1863): -0.358 (se 0.234) p=0.127cat(sprintf("2005-2024 (n=%d): %+.3f (se %.3f) p=%.3f\n",
nrow(d2), ct5e["lpolpc","Estimate"], ct5e["lpolpc","Std. Error"],
ct5e["lpolpc","Pr(>|t|)"]))## 2005-2024 (n=4304): -0.014 (se 0.107) p=0.899Rather than eyeballing split-sample estimates, this tests directly whether the police coefficient differs before and after 2015 using an interaction term in a single pooled regression on the full 2005–2024 sample.
d2$post2015 <- as.integer(d2$year >= 2015)
pd2b <- pdata.frame(d2, index = c("ori", "year"))
m_int <- plm(as.formula(paste(
"lcrmrte_total ~ lpolpc + lpolpc:post2015 +", ctrl)),
data = pd2b, model = "within", effect = "twoways")
cat("===== Police x post-2015 interaction =====\n")## ===== Police x post-2015 interaction =====print(clse(m_int)[c("lpolpc", "lpolpc:post2015"), ])## Estimate Std. Error t value Pr(>|t|)
## lpolpc 0.002623913 0.10582728 0.0247943 0.9802203
## lpolpc:post2015 -0.037697857 0.03123217 -1.2070202 0.2274953cat("\nIf 'lpolpc:post2015' is significant (p<0.10), the relationship changed\n")##
## If 'lpolpc:post2015' is significant (p<0.10), the relationship changedcat("between eras. If not, we cannot reject a stable coefficient across the\n")## between eras. If not, we cannot reject a stable coefficient across thecat("full 2005-2024 period.\n")## full 2005-2024 period.d_early <- d2 %>% filter(year <= 2014)
d_late <- d2 %>% filter(year >= 2015)
fit_era <- function(df, label) {
pdf <- pdata.frame(df, index = c("ori", "year"))
m <- plm(as.formula(paste("lcrmrte_total ~ lpolpc +", ctrl)),
data = pdf, model = "within", effect = "twoways")
ct <- clse(m)
cat(sprintf("%-12s %+.3f (se %.3f) p=%.3f [n=%d, cities=%d]\n",
label, ct["lpolpc","Estimate"], ct["lpolpc","Std. Error"],
ct["lpolpc","Pr(>|t|)"], nrow(df), n_distinct(df$ori)))
}
cat("===== Split-sample regressions (sanity check on the interaction test) =====\n")## ===== Split-sample regressions (sanity check on the interaction test) =====fit_era(d_early, "2005-2014:")## 2005-2014: +0.054 (se 0.055) p=0.331 [n=2201, cities=257]fit_era(d_late, "2015-2024:")## 2015-2024: -0.338 (se 0.197) p=0.086 [n=2103, cities=256]| Specification | Coefficient | SE |
|---|---|---|
| (1) Police only | 0.511 | 0.083 |
| (2) + controls | 0.184 | 0.099 |
| (3) + entity FE | -0.250 | 0.244 |
| (4) + time FE | 0.250 | 0.096 |
| (5) Two-way FE [baseline] | -0.358 | 0.234 |
| FE-2SLS (lagged police) | -0.624 | 0.416 |