If gun ownership causally increases homicide, it should predict firearm homicides more strongly than non-firearm homicides. Non-firearm homicides share the same confounders (poverty, inequality, demographics) but should not be mechanistically affected by gun availability. This is a standard falsification/placebo test design.
We test this at three levels:
Data prepared by scripts/04_falsification_prep.R.
Homicide counts from CDC WONDER (1999–2020, pooled). All homicides:
ICD-10 X85–Y09, Y87.1, U01–U03. Firearms homicides: X93–X95, U01.4.
Non-firearm = subtraction. Suppressed counties imputed using state-level
totals.
state_d <- readRDS(file.path(data_dir, "falsification_state.rds"))
mmsa_d <- readRDS(file.path(data_dir, "falsification_mmsa.rds"))
county_d <- readRDS(file.path(data_dir, "falsification_county.rds"))
cat("State N:", nrow(state_d |> filter(!is.na(rand_hfr))), "\n")## State N: 50cat("MMSA N:", nrow(mmsa_d), "\n")## MMSA N: 119cat("County N:", nrow(county_d), "\n")## County N: 3074For each level, we regress total, firearm, and non-firearm homicide rates on race shares + gun ownership (all standardized). If gun ownership causally increases homicide, it should have a larger coefficient for firearm than non-firearm homicide.
ms_table <- function(models, title, filename) {
modelsummary(models,
estimate = "{estimate}{stars}",
statistic = "({std.error})",
stars = c("*" = .05, "**" = .01, "***" = .001),
gof_map = c("adj.r.squared", "nobs"),
coef_rename = c("black_share" = "Black %", "hispanic_share" = "Hispanic %",
"native_share" = "Native %", "asian_share" = "Asian %",
"firearm_ownership" = "HH firearm own."),
title = title,
output = file.path(fig_dir, filename))
}
fit_models <- function(df, gun_var = "firearm_ownership") {
f <- as.formula(paste("y ~ black_share + hispanic_share + native_share + asian_share +", gun_var))
df_z <- df |> mutate(across(c(all_of(gun_var), black_share, hispanic_share, native_share, asian_share), scale))
list(
"Total" = lm(update(f, rate_all ~ .), data = df_z),
"Firearm" = lm(update(f, rate_fire ~ .), data = df_z),
"Non-firearm" = lm(update(f, rate_nonfire ~ .), data = df_z),
"Fire share" = lm(update(f, fire_share ~ .), data = df_z)
)
}mods_state <- fit_models(state_d |> filter(!is.na(rand_hfr)) |> rename(firearm_ownership = rand_hfr))
ms_table(mods_state, "State-level (RAND HFR): Homicide by weapon type", "table_falsification_state_rand.png")mods_mmsa <- fit_models(mmsa_d |> rename(firearm_ownership = brfss_mmsa))
ms_table(mods_mmsa, "MMSA-level (BRFSS): Homicide by weapon type", "table_falsification_mmsa.png")mods_county <- fit_models(county_d)
ms_table(mods_county, "County-level (SVM): Homicide by weapon type", "table_falsification_county.png")The proper test stacks firearm and non-firearm rates as a long dataset, interacting weapon type with all predictors. A significant gun ownership × weapon type interaction means gun ownership differentially predicts firearm vs non-firearm homicide. Clustered standard errors account for within-unit dependence.
run_interaction_test <- function(df, gun_var, cluster_var) {
df_long <- bind_rows(
df |> mutate(rate = rate_fire, is_firearm = 1L),
df |> mutate(rate = rate_nonfire, is_firearm = 0L)
) |> mutate(across(c(all_of(gun_var), black_share, hispanic_share, native_share, asian_share), scale))
f <- as.formula(paste("rate ~ (black_share + hispanic_share + native_share + asian_share +",
gun_var, ") * is_firearm"))
m <- lm(f, data = df_long)
vcov_cl <- vcovCL(m, cluster = df_long[[cluster_var]])
coeftest(m, vcov = vcov_cl)
}cat("=== State: RAND HFR ===\n")## === State: RAND HFR ===run_interaction_test(state_d |> filter(!is.na(rand_hfr)),
"rand_hfr", "state_fips") |> round(4)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.8545 0.0500 37.11 <2e-16 ***
## black_share 0.3735 0.0480 7.78 <2e-16 ***
## hispanic_share 0.1371 0.0678 2.02 0.0463 *
## native_share 0.4019 0.1011 3.97 0.0001 ***
## asian_share 0.0938 0.0327 2.87 0.0051 **
## rand_hfr 0.1692 0.0728 2.32 0.0225 *
## is_firearm 1.8740 0.1334 14.05 <2e-16 ***
## black_share:is_firearm 1.4976 0.1664 9.00 <2e-16 ***
## hispanic_share:is_firearm 0.4094 0.1638 2.50 0.0143 *
## native_share:is_firearm -0.2912 0.1970 -1.48 0.1428
## asian_share:is_firearm 0.0011 0.0745 0.01 0.9882
## rand_hfr:is_firearm 0.5830 0.1778 3.28 0.0015 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cat("=== State: BRFSS direct ===\n")## === State: BRFSS direct ===run_interaction_test(state_d |> filter(!is.na(brfss_direct)),
"brfss_direct", "state_fips") |> round(4)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.8529 0.0497 37.28 <2e-16 ***
## black_share 0.3767 0.0479 7.86 <2e-16 ***
## hispanic_share 0.1743 0.0683 2.55 0.0125 *
## native_share 0.3847 0.0964 3.99 0.0001 ***
## asian_share 0.0213 0.0443 0.48 0.6325
## brfss_direct 0.2131 0.0731 2.91 0.0046 **
## is_firearm 1.9003 0.1430 13.29 <2e-16 ***
## black_share:is_firearm 1.5207 0.1726 8.81 <2e-16 ***
## hispanic_share:is_firearm 0.2780 0.1451 1.92 0.0588 .
## native_share:is_firearm -0.1787 0.2212 -0.81 0.4213
## asian_share:is_firearm -0.1676 0.1597 -1.05 0.2970
## brfss_direct:is_firearm 0.3607 0.1687 2.14 0.0354 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cat("=== MMSA: BRFSS survey ===\n")## === MMSA: BRFSS survey ===run_interaction_test(mmsa_d, "brfss_mmsa", "mmsa_code") |> round(4)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.8338 0.0370 49.56 <2e-16 ***
## black_share 0.4922 0.0374 13.14 <2e-16 ***
## hispanic_share 0.2049 0.0411 4.99 <2e-16 ***
## native_share 0.4237 0.0279 15.16 <2e-16 ***
## asian_share 0.0938 0.0254 3.70 0.0003 ***
## brfss_mmsa 0.1791 0.0441 4.06 0.0001 ***
## is_firearm 2.2144 0.1399 15.83 <2e-16 ***
## black_share:is_firearm 2.0335 0.2093 9.72 <2e-16 ***
## hispanic_share:is_firearm 0.1910 0.1198 1.59 0.1122
## native_share:is_firearm -0.2204 0.0974 -2.26 0.0246 *
## asian_share:is_firearm -0.0277 0.1063 -0.26 0.7949
## brfss_mmsa:is_firearm 0.1034 0.1436 0.72 0.4721
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The interaction test above uses raw homicide rates as the outcome. However, firearm homicide rates are roughly twice as large as non-firearm rates on average. This means that even if gun ownership had an equal proportional effect on both outcomes, the absolute coefficient for firearm homicide would be ~2x larger, producing a spurious interaction.
state_sub <- state_d |> filter(!is.na(rand_hfr))
cat("Mean firearm rate:", mean(state_sub$rate_fire), "\n")## Mean firearm rate: 3.73cat("Mean non-firearm rate:", mean(state_sub$rate_nonfire), "\n")## Mean non-firearm rate: 1.85cat("Ratio:", mean(state_sub$rate_fire) / mean(state_sub$rate_nonfire), "\n")## Ratio: 2.01To address this, we re-run the interaction test with outcomes expressed as rate ratios (each state’s rate divided by the national population-weighted mean), putting both outcomes on a comparable proportional scale.
run_interaction_rr <- function(df, gun_var, cluster_var, pop_var) {
nat_fire <- weighted.mean(df$rate_fire, df[[pop_var]])
nat_nonfire <- weighted.mean(df$rate_nonfire, df[[pop_var]])
df_long <- bind_rows(
df |> mutate(rate = rate_fire / nat_fire, is_firearm = 1L),
df |> mutate(rate = rate_nonfire / nat_nonfire, is_firearm = 0L)
) |> mutate(across(c(all_of(gun_var), black_share, hispanic_share, native_share, asian_share), scale))
f <- as.formula(paste("rate ~ (black_share + hispanic_share + native_share + asian_share +",
gun_var, ") * is_firearm"))
m <- lm(f, data = df_long)
vcov_cl <- vcovCL(m, cluster = df_long[[cluster_var]])
coeftest(m, vcov = vcov_cl)
}
cat("=== State: RAND HFR (rate ratios) ===\n")## === State: RAND HFR (rate ratios) ===run_interaction_rr(state_sub, "rand_hfr", "state_fips", "hom_pop") |> round(4)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.0047 0.0271 37.11 <2e-16 ***
## black_share 0.2023 0.0260 7.78 <2e-16 ***
## hispanic_share 0.0743 0.0368 2.02 0.0463 *
## native_share 0.2178 0.0548 3.97 0.0001 ***
## asian_share 0.0508 0.0177 2.87 0.0051 **
## rand_hfr 0.0917 0.0395 2.32 0.0225 *
## is_firearm -0.1124 0.0335 -3.35 0.0012 **
## black_share:is_firearm 0.2454 0.0352 6.97 <2e-16 ***
## hispanic_share:is_firearm 0.0565 0.0483 1.17 0.2448
## native_share:is_firearm -0.1913 0.0713 -2.68 0.0087 **
## asian_share:is_firearm -0.0281 0.0214 -1.32 0.1917
## rand_hfr:is_firearm 0.0884 0.0519 1.70 0.0923 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cat("=== State: BRFSS direct (rate ratios) ===\n")## === State: BRFSS direct (rate ratios) ===run_interaction_rr(state_d |> filter(!is.na(brfss_direct)),
"brfss_direct", "state_fips", "hom_pop") |> round(4)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.9885 0.0265 37.28 <2e-16 ***
## black_share 0.2010 0.0256 7.86 <2e-16 ***
## hispanic_share 0.0930 0.0364 2.55 0.0125 *
## native_share 0.2053 0.0514 3.99 0.0001 ***
## asian_share 0.0114 0.0236 0.48 0.6325
## brfss_direct 0.1137 0.0390 2.91 0.0046 **
## is_firearm -0.1036 0.0351 -2.95 0.0041 **
## black_share:is_firearm 0.2464 0.0357 6.91 <2e-16 ***
## hispanic_share:is_firearm 0.0137 0.0404 0.34 0.7360
## native_share:is_firearm -0.1567 0.0737 -2.13 0.0364 *
## asian_share:is_firearm -0.0459 0.0415 -1.11 0.2723
## brfss_direct:is_firearm 0.0216 0.0455 0.47 0.6364
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1We also confirm with fully standardized betas (both predictors and outcomes z-scored), which directly compare effect sizes across outcomes:
state_zz <- state_sub |>
mutate(across(c(rand_hfr, black_share, hispanic_share, native_share, asian_share,
rate_fire, rate_nonfire), scale))
m_fire_z <- lm(rate_fire ~ black_share + hispanic_share + native_share + asian_share +
rand_hfr, data = state_zz)
m_nonfire_z <- lm(rate_nonfire ~ black_share + hispanic_share + native_share + asian_share +
rand_hfr, data = state_zz)
modelsummary(list("Firearm (std)" = m_fire_z, "Non-firearm (std)" = m_nonfire_z),
estimate = "{estimate}{stars}",
statistic = "({std.error})",
stars = c("*" = .05, "**" = .01, "***" = .001),
gof_map = c("adj.r.squared", "nobs"),
coef_rename = c("black_share" = "Black %", "hispanic_share" = "Hispanic %",
"native_share" = "Native %", "asian_share" = "Asian %",
"rand_hfr" = "HH firearm own. (RAND)"),
title = "Fully standardized coefficients (state-level, RAND HFR)",
output = file.path(fig_dir, "table_falsification_fully_standardized.png"))results <- tribble(
~Level, ~`Gun measure`, ~N, ~`Raw interaction p`, ~`RR interaction p`,
"State", "RAND HFR", 50, 0.0015, 0.092,
"State", "BRFSS direct", 48, 0.035, NA,
"MMSA", "BRFSS survey", 119, 0.47, NA,
"County", "SVM imputed", 3074, NA, NA
)
knitr::kable(results, caption = "Gun ownership x weapon type interaction (clustered SEs)")| Level | Gun measure | N | Raw interaction p | RR interaction p |
|---|---|---|---|---|
| State | RAND HFR | 50 | 0.002 | 0.092 |
| State | BRFSS direct | 48 | 0.035 | NA |
| MMSA | BRFSS survey | 119 | 0.470 | NA |
| County | SVM imputed | 3074 | NA | NA |
The falsification test yields inconclusive results:
A key complication is that racial composition dominates the weapon-type split: the Black share x weapon type interaction is massive at every level (beta ~ 1.5, p < .0001), reflecting that areas with larger Black populations have a disproportionately higher share of firearm (vs. non-firearm) homicides. This confounding structure leaves little residual variance for gun ownership to explain, and limits the power of this falsification design.
sessionInfo()## R version 4.5.3 (2026-03-11)
## Platform: x86_64-pc-linux-gnu
## Running under: Linux Mint 21.1
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0 LAPACK version 3.10.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_DK.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_DK.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_DK.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Europe/Brussels
## tzcode source: system (glibc)
##
## attached base packages:
## [1] grid stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] haven_2.5.5 survey_4.4-8 survival_3.8-6
## [4] foreign_0.8-91 sandwich_3.1-1 lmtest_0.9-40
## [7] zoo_1.8-14 readxl_1.4.5 tidycensus_1.7.5
## [10] kirkegaard_2025-11-19 robustbase_0.99-6 psych_2.5.6
## [13] assertthat_0.2.1 weights_1.1.2 magrittr_2.0.4
## [16] lubridate_1.9.4 forcats_1.0.1 stringr_1.6.0
## [19] purrr_1.2.0 tidyr_1.3.1 tibble_3.3.0
## [22] tidyverse_2.0.0 gt_1.3.0 effectsize_1.0.1
## [25] spdep_1.4-1 spatialreg_1.4-2 spData_2.3.4
## [28] car_3.1-3 carData_3.0-5 fixest_0.13.2
## [31] e1071_1.7-16 xgboost_1.7.11.1 glmnet_4.1-10
## [34] Matrix_1.7-4 readr_2.1.5 modelsummary_2.5.0
## [37] ggplot2_4.0.1.9000 tigris_2.2.1 sf_1.0-21
## [40] dplyr_1.1.4
##
## loaded via a namespace (and not attached):
## [1] splines_4.5.3 later_1.4.4 tinytable_0.14.0
## [4] cellranger_1.1.0 datawizard_1.3.0 rpart_4.1.24
## [7] lifecycle_1.0.4 httr2_1.2.1 Rdpack_2.6.4
## [10] ellmer_0.4.0 globals_0.18.0 processx_3.8.6
## [13] lattice_0.22-7 vroom_1.6.6 MASS_7.3-65
## [16] insight_1.4.2 backports_1.5.0 btw_1.1.0
## [19] Hmisc_5.2-4 sass_0.4.10 rmarkdown_2.30
## [22] jquerylib_0.1.4 yaml_2.3.10 sp_2.2-0
## [25] minqa_1.2.8 chromote_0.5.1 DBI_1.2.3
## [28] RColorBrewer_1.1-3 multcomp_1.4-28 abind_1.4-8
## [31] rvest_1.0.5 nnet_7.3-20 coro_1.1.0
## [34] TH.data_1.1-4 rappdirs_0.3.3 webshot2_0.1.2
## [37] dreamerr_1.5.0 listenv_0.9.1 gdata_3.0.1
## [40] units_1.0-0 performance_0.15.2 parallelly_1.45.1
## [43] codetools_0.2-20 xml2_1.4.0 tidyselect_1.2.1
## [46] shape_1.4.6.1 farver_2.1.2 lme4_1.1-37
## [49] base64enc_0.1-3 jsonlite_2.0.0 mitml_0.4-5
## [52] Formula_1.2-5 iterators_1.0.14 emmeans_1.11.2-8
## [55] systemfonts_1.3.1 foreach_1.5.2 tools_4.5.3
## [58] ragg_1.5.0 Rcpp_1.1.0 glue_1.8.0
## [61] mnormt_2.1.1 pan_1.9 gridExtra_2.3
## [64] xfun_0.57 mcptools_0.2.0 mgcv_1.9-3
## [67] websocket_1.4.4 withr_3.0.2 numDeriv_2016.8-1.1
## [70] fastmap_1.2.0 mitools_2.4 boot_1.3-32
## [73] fansi_1.0.6 digest_0.6.39 timechange_0.3.0
## [76] R6_2.6.1 estimability_1.5.1 mice_3.18.0
## [79] colorspace_2.1-2 textshaping_1.0.4 wk_0.9.4
## [82] LearnBayes_2.15.1 gtools_3.9.5 utf8_1.2.6
## [85] generics_0.1.4 data.table_1.17.8 class_7.3-23
## [88] htmlwidgets_1.6.4 httr_1.4.7 parameters_0.28.2
## [91] pkgconfig_2.0.3 gtable_0.3.6 rsconnect_1.5.1
## [94] S7_0.2.1 brio_1.1.5 htmltools_0.5.8.1
## [97] scales_1.4.0 nanonext_1.8.0 reformulas_0.4.1
## [100] knitr_1.50
## [ reached 'max' / getOption("max.print") -- omitted 48 entries ]