The impact of fentanyl on state- and county-level psychostimulant and cocaine overdose death rates by race in Ohio from 2010 to 2020: a time series and spatiotemporal analysis

Background Over the past decade in the USA, increases in overdose rates of cocaine and psychostimulants with opioids were highest among Black, compared to White, populations. Whether fentanyl has contributed to the rise in cocaine and psychostimulant overdoses in Ohio is unknown. We sought to measure the impact of fentanyl on cocaine and psychostimulant overdose death rates by race in Ohio. Methods We conducted time series and spatiotemporal analyses using data from the Ohio Public Health Information Warehouse. Primary outcomes were state- and county-level overdose death rates from 2010 to 2020 for Black and White populations. Measures of interest were overdoses consisting of four drug involvement classes: (1) all cocaine overdoses, (2) cocaine overdoses not involving fentanyl, (3) all psychostimulant overdoses, and (4) psychostimulant overdoses not involving fentanyl. We fit a time series model of log standardized mortality ratios (SMRs) using a Bayesian generalized linear mixed model to estimate posterior median rate ratios (RR). We conducted a spatiotemporal analysis by modeling the SMR for each drug class at the county level to characterize county-level variation over time. Results In 2020, the greatest overdose rates involved cocaine among Black (24.8 deaths/100,000 people) and psychostimulants among White (10.1 deaths/100,000 people) populations. Annual mortality rate ratios were highest for psychostimulant-involved overdoses among Black (aRR = 1.71; 95% CI (1.43, 2.02)) and White (aRR = 1.60, 95% CI (1.39, 1.80)) populations. For cocaine not involving fentanyl, annual mortality rate ratios were similar among Black (aRR = 1.04; 95% CI (0.96,1.16)) and White (aRR = 1.02; 95% CI (0.87, 1.20)) populations. Within each drug category, change over time was similar for both racial groups. The spatial models highlighted county-level variation for all drug categories. Conclusions Without the involvement of fentanyl, cocaine overdoses remained constant while psychostimulant overdoses increased. Tailored harm reduction approaches, such as distribution of fentanyl test strips and the removal of punitive laws that influence decisions to contact emergency services, are the first steps to reduce cocaine overdose rates involving fentanyl among urban populations in Ohio. In parallel, harm reduction policies to address the increase in psychostimulant overdoses are warranted. Supplementary Information The online version contains supplementary material available at 10.1186/s12954-024-00936-9.


Statistical Models
For reference and completeness, the Bayesian hierarchical models are detailed here.These models were fit independently for each of the four drug categories present in the analysis using state-level and county-level overdose death count data obtained from the Ohio Public Health Information Warehouse.In the following exposition let s indicate county (s = 1, . . ., 88), t year (t = 2010, ..., 2020) and r race (r = 0 is White and r = 1 is Black).
Therefore the total expected number of overdose deaths E tr for year t and race r, based upon the age-specific overdose death rate in the 2010 Ohio population, is It follows that λ tr represents the state-level relative risk of overdose death for race r in year t compared to the corresponding risk for race r in year t expected from the age-specific overdose death rate in the 2010 Ohio population, commonly referred to as the standardized mortality ratio (SMR).Using the log link, we model the SMR as where β r is a vector of race-specific fixed effects corresponding to the vector of covariates X t containing an intercept, years since 2010 (i.e.t − 2010), an indicator of race (with White as the reference level) and the interaction of years since 2010 and race.We also include a term for the remaining race-specific heterogeneity ϵ tr to capture variation not accounted for by the fixed effects.We assume an auto-regression of order 1 structure to account for temporal autocorrelation.Thus, we assume where ϕ r is a race-specific auto-regressive parameter and σ 2 r is a variance.The fixed effects are assigned independent flat priors.Auto-regressive parameters ϕ r are assumed to be uniformly distributed over (0,1).All variance parameters are assigned independent inverse gamma distributions with shape and scale parameters both equal to 0.5.The models were fit using Markov chain Monte Carlo implemented in NIMBLE [1] in R. The algorithm was run for 1,000,000 iterations, discarding the first 500,000 as burn-in, and then thinned by keeping every 100 th iteration.We assessed convergence visually using trace plots.

Spatio-temporal model
For a detailed examination of the county-level overdose death count data, we modify the multivariate Poisson linear mixed effects model introduced by Kline, Pan, and Hepler (2021) [2] to fit our current problem setting.This model complements the state-level analysis by providing insights into county-level overdose death patterns.It assumes that where Y str is number of overdose deaths for county s, year t and race r and E str is the corresponding expected number of overdose deaths for county s, year t and race r based upon the age-specific overdose death rate in the 2010 Ohio population.As in the time-series analysis, this formulation assumes the expected overdose rates between racial groups are homogeneous.E str is defined in an analogous manner to E tr in the time-series model, but for the county-level.Follow the procedure detailed in the time series model and obtain E stri by summing over the age groups: It follows that λ str represents the relative risk of overdose death for race r in county s and year t compared to the corresponding risk for race r in county s and year t expected from the age-specific death rate in the 2010 Ohio population.
Using the log link, we model the SMR as where γ r is a vector of race-specific fixed effects corresponding to the vector of covariates W st containing an intercept and years since 2010.We also include a spatio-temporal shared component ν st with race-specific loading δ r and a term for the remaining race-specific heterogeneity ϵ str .For identification, we assume the loading for White Ohioans is equal to 1 (i.e.δ 0 = 1).Since racial groups within a county share a common environment, we include a shared component ν st in the model to reflect this.We also use this term to account for spatial structure with an intrinsic conditional auto-regressive model and shared temporal auto-correlation using an auto-regression of order 1.That is, we assume the following conditional distributions: where η is a temporal auto-regressive parameter, w sℓ is an indicator of whether counties s and ℓ are neighbors as defined by adjacency, w s+ is the number of neighbors of county s, and τ 2 is a variance parameter.We constrain s ν st = 0 for every year t which allows us to use equation ( 9) as a valid process model [3].
Acknowledging that racial groups may experience the same environment differently, we include a race-specific loading δ r to allow for differing effects of the shared component by race.We also model race-specific variation ϵ str to capture variation not accounted for by the shared component or fixed effects.We assume an auto-regression of order 1 structure to account for temporal autocorrelation within a county.Thus, we assume where θ r is a race-specific auto-regressive parameter and ω 2 r is a variance.The fixed effects and δ r are assigned independent flat priors.Auto-regressive parameters η and θ r are assumed to be uniformly distributed over (0,1).All variance parameters are assigned independent inverse gamma distributions with shape and scale parameters both equal to 0.5.The models were fit using Markov chain Monte Carlo implemented in NIMBLE in R. For the spatio-temporal model the algorithm was run for 1,000,000 iterations, discarding the first 500,000 as burn-in, and then thinned by keeping every 100 th iteration.We assessed convergence visually using trace plots.

1 2 ws+
ws+ ℓ w sℓ ν ℓt , τ 2 ws+ for t = 2010 ν st |ν −st , ν s(t−1) ∼ N ην s(t−1) + 1 ws+ ℓ w sℓ (ν ℓt − ην ℓ(t−1) ), τ 1.1 Time-series modelA multivariate Poisson linear mixed effects model was used to perform the timeseries analysis of the state-level overdose death count data.Assume that tr is number of overdose deaths in Ohio for year t and race r and E tr is the corresponding expected number of overdose deaths in Ohio for year t and race r based upon the age-specific overdose death rate in the 2010 Ohio population.In other words, E tr is the number of age-adjusted overdose deaths in Ohio for year t and race r obtained via indirect-standardization with the 2010 Ohio population serving as the reference population.This formulation assumes the expected overdose rates between racial groups are homogeneous.E tr is defined as follows.First, denote the 2010 state-level overdose death rate for age group i as ρ i where i ∈ {≤ 1, 1-4, 5-14, 15-24,..., 75-84, 85+}.Let P stri be the resident population indexed by county, year, race and age group.The corresponding expected number of overdose deaths for county s, year t, race r and age group i is