Forecasting state overdose deaths through the surge

Can state economic and demographic data from 2014–2018 forecast overdose deaths in 2019–2022? Poisson GLMs and spline GAMs against a hierarchical Bayesian negative-binomial model, all fit in base R.

Robert M. ShollKansas State University

Summary

How well do state-level data from 2014–2018 forecast overdose deaths in 2019–2022, the years of the fentanyl surge? I compare five likelihood models (Poisson GLMs and spline GAMs with state fixed effects) with a hierarchical Bayesian negative-binomial model.

Every model is trained on 2014–2018 and scored on 2019–2022. The best forecaster is GLM 2, with state effects and the LASSO-selected covariates (RMSE 508 deaths, median absolute error 13% of a state's count). The Bayesian model with one shared national trend is close behind (RMSE 521) and has the best-calibrated intervals (77% of outcomes inside its 90% band). Giving every state its own pre-2019 trend makes forecasts worse (RMSE 959): five years are too few to extrapolate a state-specific slope four years ahead. No model reaches 90% coverage. Every one underpredicts 2020–2022, because nothing in these covariates anticipates fentanyl saturating the drug supply. State identity carries most of the signal; the economic covariates add little.

Data

overdose_states.csv: 50 states (no DC) × 2014–2022. The fields are overdose deaths, the published age-adjusted rate per 100,000, population (thousands), GDP, per-capita income, homeownership, percent with a bachelor's degree, and state tax collections. The rates match CDC all-drug overdose death rates (for example West Virginia 2021 at 90.9), so the response is all overdose deaths, not fentanyl-only deaths.

The source file swaps the 2014 death counts for Arizona and Arkansas; the values below are corrected.
StateYearPublished rateDeaths in fileImplied by rateDeaths used
AR201418.201,211540356
AZ201412.603568481,211

The crude rate over the published rate ranges 0.66–1.43 across all 450 rows (median 0.97). That is the expected gap between crude and age-adjusted rates.

2014 2015 2016 2017 2018 2019 2020 2021 2022 0 20 40 60 80 test period Nebraska California median state Delaware Tennessee West Virginia Deaths per 100,000, age-adjusted Overdose death rates by state, 2014–2022
Figure 1. Grey lines are all 50 states; the black line is the median state. The 2020 jump falls entirely inside the test period.

Variable selection

Candidates are per-capita GDP, per-capita income, per-capita state tax revenue, homeownership, bachelor's degree share, and year. The response is log deaths per resident in the training years. Folds are whole states, so no state appears in both the fitting and scoring folds. The one-standard-error rule picks the sparsest model within one SE of the best.

−9 −8 −7 −6 −5 −4 −3 log λ (penalty falls to the right) −0.02 0.00 0.02 0.04 homeownership GDP tax revenue income bachelor’s % year Standardized coefficient Coefficient paths −9 −8 −7 −6 −5 −4 −3 log λ (penalty falls to the right) 12.5 15.0 17.5 20.0 22.5 25.0 27.5 one-SE choice Mean squared error × 1,000, ±1 SE Cross-validated error
Figure 2. Left: coefficient paths as the penalty relaxes. Right: grouped five-fold cross-validated error with ±1 SE. The dotted line marks the minimum; the dashed line marks the one-SE choice, which keeps year and bachelor’s share.
Standardized LASSO coefficients for log death rate.
VariableCoefficient at λ(1se)Coefficient at λ(min)Selected
GDP per capita00
Income per capita00.0198
Tax revenue per capita00
Homeownership %0-0.0235
Bachelor's degree %0.02060.0319yes
Year0.04140.0457yes

The LASSO keeps bachelor’s % and year. A second set read off a causal diagram (GDP, bachelor’s %, year) is kept as a theory-driven alternative.

Models

All models use log(population) as an offset, so coefficients describe death rates.
ModelStructureFamily
GLM 1 · intercept onlylog μ = log pop + αquasi-Poisson
GLM 2 · state + LASSO setlog μ = log pop + α_state + LASSO covariatesquasi-Poisson
GLM 3 · state + DAG setlog μ = log pop + α_state + GDP + bachelor's + yearquasi-Poisson
GAM 1 · state + splines (LASSO)as GLM 2 with natural cubic splines (3 df)quasi-Poisson
GAM 2 · state + splines (DAG)as GLM 3 with splines on GDP and bachelor'squasi-Poisson
Bayes · NB, state trendslog μ = log pop + a_s + b_s·year + GDP + bachelor's; a_s and b_s partially poolednegative binomial
Bayes · NB, common trendlog μ = log pop + a_s + β·year + GDP + bachelor's; a_s partially poolednegative binomial

The Bayesian models give every state its own level. The first also gives every state its own pre-2019 trend. Both are drawn toward the national average in proportion to how noisy the state's own data are. Covariates are centered within state, so their coefficients describe within-state change and do not compete with the state intercepts. Both are fit with a Metropolis-within-Gibbs sampler written for this project: two chains, 16,000 iterations, 4,000 burn-in, thinned by 6.

State-trend model: posterior summaries and convergence (R-hat near 1 and large ESS mean the chains agree). Common-trend model R-hat for its fixed effects: 1.003, 1.005, 1.009.
Parametermeanq05q95RhatESS
mu_b (mean yearly log-rate change)0.07590.04450.1071.02151
sigma_b (between-state spread of trend)0.06490.05260.07931.002,122
sigma_a (between-state spread of level)0.3990.3370.4711.004,000
beta: GDP per capita (per SD)-0.104-0.2150.001681.02228
beta: bachelor's degree % (per SD)0.0922-0.08280.2621.01174
phi (NB size)1481171831.001,990
a[WV] (state level)-7.73-7.80-7.671.002,771
0 250 500 750 1000 1250 1500 1750 2000 Saved draw 0.00 0.05 0.10 0.15 Log-rate change per year Mean yearly trend, μb 0 250 500 750 1000 1250 1500 1750 2000 Saved draw 100 125 150 175 200 225 Dispersion parameter Negative-binomial size, φ
Figure 3. Two chains (navy and ochre) for two parameters after burn-in. They overlap without drifting, consistent with the R-hat values in the table above.
WY OK MT MS UT NM AK OR SD KS CO TX ID AL IA NE WA GA RI NV CA AR KY MN NH AZ ND WI HI TN WV LA IN MO MI SC VA OH NC VT ME MA IL FL PA CT NY DE MD NJ −10% 0% 10% 20% 30% national mean 8% a year Yearly change in death rate, 2014–2018: median and 90% interval State trends before 2019
Figure 4. State-specific pre-2019 trends from the state-trend model, sorted. New Jersey, Maryland, Delaware and much of the Northeast were already climbing fastest before the test period.

Out-of-sample forecasts, 2019–2022

Forecast accuracy on the 200 held-out state-years (deaths). The best RMSE is highlighted.
ModelRMSEMAEMedian APE %90% PI coverage %
GLM 1 · intercept only1,12874038.565.5
GLM 2 · state + LASSO set50828113.476
GLM 3 · state + DAG set56029613.274.5
GAM 1 · state + splines (LASSO)51329613.075
GAM 2 · state + splines (DAG)62435316.574
Bayes · NB, state trends95949818.364
Bayes · NB, common trend52128913.377
2019 2020 2021 2022 0 200 400 600 800 1,000 1,200 Bayes, common trend GAM 1, lasso set GLM 2, lasso set GLM 3, DAG set GAM 2, DAG set GLM 1, intercept only Bayes, state trends Root-mean-square error, deaths per state Forecast error grows through 2021 for every model
Figure 5. Error grows every year from 2019 to 2021 for every model as the surge pulls deaths away from the 2014–2018 pattern.
2014 2018 2022 0 500 1,000 1,500 West Virginia 2014 2018 2022 0 2,000 4,000 6,000 8,000 10,000 Ohio 2014 2018 2022 0 2,000 4,000 6,000 8,000 10,000 California 2014 2018 2022 0 1,000 2,000 3,000 4,000 5,000 Texas 2014 2018 2022 0 200 400 600 800 1,000 Delaware 2014 2018 2022 0 200 400 600 Kansas Deaths
Figure 6. Observed deaths (dots), the state-trend Bayesian model’s posterior predictive median and 90% band (navy), and the GLM 2 forecast (ochre dashed). Training years are left of the vertical line. State trends fail both ways: states already climbing before 2019 are overshot (Ohio, Delaware), and western states that were flat before 2019 are badly undershot when fentanyl arrived (California, Texas, Kansas).