# Forecasting Weekly US Public Mass-Shooting Incidence from Lagged Public Signals

## 1. Abstract

This study tests whether a small collection of publicly observable signals can
forecast weekly United States public mass-shooting incidence. The outcome is the
number of Mother Jones/Codeholics-defined incidents in each Monday-Sunday ISO
week: a public, indiscriminate shooting with at least three fatalities. We build
a 522-week panel from 2014-W01 through 2023-W52 and reserve the final 104 weeks,
2022-W01 through 2023-W52, as a fixed test set. Predictors include trailing
incident counts, seasonality, holiday and school-session indicators, lagged FBI
National Instant Criminal Background Check System (NICS) counts, annual
population, and a completed-month incident-frequency proxy for public attention.
We compare a week-of-year seasonal mean, a lag-only negative-binomial (NB) model,
a regularized full NB model, and regularized logistic models for at least one and
at least two incidents. The full NB obtained test RMSE 0.450, slightly below the
lag-only NB at 0.452 and the seasonal baseline at 0.465, but had worse MAE than
both. For the at-least-one task, logistic average precision was 0.271 versus
0.237 for the NB probability baseline. Only one at-least-two week occurred in
the test set, making discrimination estimates extremely unstable. These are
weak predictive signals, not evidence of causality or a deployable warning
system.

## 2. Research question and contribution

The question is deliberately narrow: at the beginning of week t, how well can
national incidence in that week be forecast using information that precedes its
outcome? The exercise does not attempt to identify individual perpetrators,
locations, or causal mechanisms. It instead treats incidence as a sparse
national time series and asks whether temporal memory, broad firearm-market
activity, calendar structure, and an attention proxy improve on simple
historical baselines.

The contribution is primarily methodological and reproducible. Public
mass-shooting research is unusually sensitive to definition changes, incomplete
records, and retrospective feature construction. This experiment fixes a
definition, commits immutable source snapshots, declares a time split before
modeling, and places every outcome-linked feature behind a temporal lag. It also
reports negative results. With 82 incidents across ten years, model flexibility
is constrained by the number of events rather than the number of weekly rows.
The design therefore emphasizes transparent generalized linear models and
honest held-out evaluation rather than extensive hyperparameter search.

## 3. Outcome definition and the GVA gap

The target follows the Mother Jones database as extended and normalized by the
Codeholics project [@follman2012; @codeholics2026]. An included event has at
least three fatalities in a single public, indiscriminate incident. Shootings
primarily associated with conventional armed crime, gang activity, or exclusively
domestic circumstances are outside this operational target. The threshold
reflects the post-2013 Mother Jones convention and is compatible with the
public-mass-shooting frame used in a substantial part of the academic literature
[@lankfordmadfis2018]. The weekly count is the number of qualifying incidents,
not the number killed or injured. Two binary targets are derived without
relabeling: any incident means count at least one, and high incidence means
count at least two.

Despite the historical filename gva_mass_shootings.csv, these 82 rows are not a
Gun Violence Archive (GVA) extract. GVA commonly reports a broader event class
based on at least four people shot, injured or killed, generally excluding the
shooter, without requiring three deaths or the same public-indiscriminate
screen. That definition answers a different question and produces many more
events [@gva2026]. Calling our outcome “GVA mass shootings” would therefore be
incorrect. This definition gap also means the results cannot be compared
directly with models trained on a four-shot GVA outcome.

Weeks follow ISO Monday-Sunday boundaries. The panel begins at the Monday of
2014-W01, 2013-12-30, and ends at the Monday of 2023-W52, 2023-12-25. An
incident is assigned to the week containing its date. Empty weeks are explicit
zeros. The resulting panel contains 522 weeks and 82 incidents, an average of
0.157 incidents per week. Seventy-seven weeks have at least one incident; five
weeks across the full panel have at least two. The sparsity is intrinsic to the
chosen definition.

## 4. Data sources and the public-attention substitution

Incident dates, case names, locations, fatalities, injuries, and total victims
come from the committed Codeholics enrichment of the Mother Jones series. The
loader retains the study window, validates the minimum-fatality rule, parses
dates, and exposes only analysis columns. Geographic coordinates are not used:
several are absent, and the prediction target is national rather than spatial.

Monthly NICS national sums were prepared from the BuzzFeed News republication
of FBI files [@buzzfeednics2017; @fbinics2026]. We sum handgun, long-gun,
multiple, and other firearm categories. NICS is a background-check activity
signal, not a gun-sales count; state procedures, permits, and multiple-firearm
transactions prevent a one-to-one interpretation [@fbinics2026]. The snapshot
ends in September 2023. For the remaining study weeks, the builder carries the
last observed weekly level forward rather than fabricating October-December
monthly totals. All NICS model columns are subsequently lagged.

The file named gtrends_mass_shooting.csv is not Google Trends data. Authenticated,
stable Google Trends retrieval was unavailable for the static offline build.
The committed substitute is the monthly number of qualifying incidents from the
same incident source. It is interpreted as a direct proxy for the volume of
events capable of generating public awareness and news coverage. This is a
material limitation: it overlaps conceptually with lagged outcomes and does not
measure searches, story counts, sentiment, or audience exposure. To avoid using
a completed month's future total within that same month, each monthly proxy
value becomes available only after month end; weekly lags are applied after
that availability rule. The variable should therefore be read as delayed
incident-driven attention, not as Google search interest.

Annual population comes from World Bank indicator SP.POP.TOTL [@worldbank2026].
It enters as a slow-moving scale control using the ISO week-year. All four raw
CSV snapshots remain unchanged by the pipeline.

## 5. Temporal design and leakage controls

Development data span 2014-W01 through 2021-W52, 418 weeks. The final two
calendar years, 2022-W01 through 2023-W52, contain exactly 104 weeks and are
held out from parameter fitting. Within development data, 2021 is the declared
last chronological validation year. The backtest additionally creates expanding
origins with full-year validation blocks in 2019, 2020, and 2021. Those folds
diagnose temporal stability; they are not mixed into the final test tables.

Outcome rolling features shift the count by one week before summation. Thus a
four-week trailing count for week t contains weeks t-1 through t-4 and never
week t. The same-week-of-year feature is a 52-week lag. The seasonal baseline
uses up to five earlier matching ISO weeks and falls back to the development
grand mean when fewer than two exist. During a rolling forecast, previously
observed outcomes may enter later lags, which matches a real weekly updating
setting, but test outcomes are never used to estimate coefficients.

Each observed NICS month is allocated to days with a piecewise-linear triangular
window whose explicit anchors are the month start, midpoint, and end. Weights
sum to the monthly total and never interpolate to an adjacent month's value.
Daily allocations are then summed within ISO weeks. This satisfies the requested
within-month interpolation, although it should not be interpreted as observed
daily transactions. Lags of 1, 4, 13, 26, and 52 weeks are formed. The proxy
availability delay described in Section 4 is followed by 1, 2, 4, 8, 13, 26,
and 52 week lags. Initial unavailable lag history is set to a neutral zero, and
NICS boundary values are carried, so no row after 2014-W03 has missing values.

## 6. Predictive features

Incident-memory predictors are trailing totals over 1, 2, 4, 8, 13, 26, and 52
weeks, plus the count 52 weeks earlier. The specified 1, 4, 13, 26, and 52 week
windows are available to every full model; the added 2 and 8 week windows make
the lag-only NB specification complete. These variables represent persistence,
clustering, or regression to a recent rate without assuming any mechanism.

Calendar predictors encode month, day of year, and ISO week of year with sine
and cosine pairs. Holiday flags identify the ISO week containing January 1,
July 4, Thanksgiving, and Christmas/New Year. The school-session variable is
the fraction of Monday-Friday dates in a week that fall between August 15 and
May 31. It is therefore one for most school-year weeks, zero for most summer
weeks, and fractional at boundaries. Population and all lagged NICS and proxy
variables complete the full design.

Before GLM estimation, numeric predictors are converted to finite values using
training medians, constant columns are removed, and exact duplicate columns are
dropped. Means and standard deviations are estimated on the training partition
only, then applied to the forecast partition. Standardized values are clipped
at ten standard deviations to limit numerical leverage. Removing exact calendar
aliases is a numerical de-duplication step, not outcome-based feature selection.

## 7. Models

The seasonal baseline predicts the average count for matching week of year over
the prior five years, requiring at least two historical matches. Its predictive
quantiles use a Poisson distribution around that mean. The lag-only model is a
statsmodels generalized linear model with a log link and negative-binomial NB2
family. Its seven predictors are the trailing incident windows. Dispersion is a
bounded method-of-moments estimate from the fitting outcome. NB predictive
quantiles follow variance mu plus alpha times mu squared.

The full NB uses every available feature and light ridge regularization to
stabilize estimation in a sparse panel. Two separate binomial GLMs estimate the
probability of count at least one and count at least two. Their slopes receive
ridge regularization while the intercept remains unpenalized, preserving the
rare-event base rate. For binary baselines, the lag-only NB predictive
distribution is converted to probabilities of at least one or at least two.

A random forest was conditionally specified but not included. The build rule
allowed it only if importing scikit-learn completed in under one second; the
measured cold import was 1.54 seconds. The omission was made before results were
examined and is recorded in the README and model metadata.

## 8. Evaluation measures

Count forecasts are evaluated using root mean squared error (RMSE), mean absolute
error (MAE), stabilized mean absolute percentage error (MAPE), and pinball loss
at quantiles 0.5, 0.8, and 0.9. Because most true weekly counts are zero, ordinary
MAPE is undefined. The reported stabilized version uses max of the absolute
outcome and one as denominator and is labeled explicitly in code.

Binary forecasts are assessed with area under the precision-recall curve
(average precision), area under the ROC curve, and Brier score. Precision-recall
performance is especially relevant under imbalance, while Brier score evaluates
probability accuracy. With only one positive high-incidence test week, neither
ROC nor average precision supports a stable ranking conclusion. No threshold
accuracy is reported because it would reward an always-negative classifier.

## 9. Held-out results

| Count model | RMSE | MAE | MAPE (%) | Pinball .50 | Pinball .80 | Pinball .90 |
|---|---:|---:|---:|---:|---:|---:|
| Seasonal | 0.465 | 0.310 | 30.00 | 0.115 | 0.187 | 0.109 |
| Lag-only NB | 0.452 | 0.319 | 31.02 | 0.115 | 0.194 | 0.083 |
| Full NB | 0.450 | 0.344 | 33.50 | 0.115 | 0.167 | 0.082 |

The full NB's RMSE improvement over the seasonal model is 0.015 in absolute
terms, while its MAE is worse by 0.034. It improves the upper-quantile losses but
not the median loss. This mixed pattern does not support a broad claim that the
full feature set dominates simple history.

| Any-incident model | AUC-PR | AUC-ROC | Brier |
|---|---:|---:|---:|
| Lag-only NB probability | 0.237 | 0.557 | 0.1793 |
| Full logistic | 0.271 | 0.553 | 0.1769 |

The logistic model modestly improves average precision and Brier score, while
ROC area is essentially unchanged and only slightly above chance. There were 23
positive any-incident weeks among 104 test weeks.

| At-least-two model | AUC-PR | AUC-ROC | Brier |
|---|---:|---:|---:|
| Lag-only NB probability | 0.0213 | 0.553 | 0.00969 |
| Full logistic | 0.0130 | 0.262 | 0.00963 |

There was one high-incidence test week. The logistic Brier score is marginally
lower because its probabilities stay near the rare base rate, but its ranking
places the single positive poorly. These numbers are descriptive of one event,
not reliable comparative estimates.

## 10. Negative findings, limitations, and ethics

Several requested signals did not produce decisive gains. The full NB's extra
features reduced RMSE only slightly and worsened MAE and stabilized MAPE. The
logistic any model improved probability error but not ROC ranking. The high
classifier is underidentified in practical terms: four positives occur in
development and one in test. A first regularized logistic implementation also
penalized the intercept, which pulled the high-event forecast toward an
implausible 12% weekly mean. The implementation was corrected during model
development to leave the intercept unpenalized; the resulting mean test
probability is about 1.2%.

The news proxy is the most important construct-validity weakness. It is a
delayed transformation of the outcome source and can capture recent event tempo,
but not attention independently of events. NICS data are monthly, and triangular
daily allocation creates timing detail not present in the source. The late-2023
carry-forward avoids invented totals but removes contemporary variation.
Population changes little over ten years. Holiday and school calendars are
national simplifications that cannot express local schedules or exposure.

Database definitions inevitably encode editorial decisions, and missing or
reclassified incidents could change a five-positive high-incidence analysis.
Forecast association must not be interpreted as a causal effect of firearm
checks, media attention, holidays, or schools. National weekly probabilities
should not be used to profile people or communities, allocate punitive policing,
or claim that a specific attack is predictable. Responsible use is limited to
auditable methodological research and aggregate preparedness discussion.

## 11. Reproducibility and conclusion

Running make all executes six tests and rebuilds every artifact from the four
committed raw snapshots. results/features.csv contains the 522-week panel;
results/predictions_test.csv contains only the 104 fixed-test forecasts and
outcomes; three metrics files separate count, any, and high tasks.
results/cv_metrics.csv preserves rolling-origin development results rather than
mixing them with test scores. Coefficient and metadata files record fitted
parameters and the conditional forest decision. Figures show observed versus
predicted counts, full-NB coefficient magnitude, residuals, any-event ROC, and
calibration. The source loaders and metrics have direct unit coverage.

The central finding is modest: recent incidence and broad public signals allow
small improvements on some held-out scoring rules, but no model is consistently
superior and the at-least-two target is too rare for stable evaluation. The
most defensible result is therefore the reproducible benchmark itself, including
its negative findings and clear separation between a narrowly defined public
mass-shooting outcome and broader GVA-style event counts.