ABSTRACT
Background Since 2012 fentanyl-involved fatal overdoses have risen from 4% of all fatal overdoses in Connecticut to 82% in 2019. We investigated the geographic and temporal trends in fentanyl-involved deaths in Connecticut during 2009-2019 and examined the relationship of fentanyl-involved deaths and overall trends of opioid overdose deaths.
Methods Data on the dates and locations of accidental/undetermined opioid-involved fatalities that occurred during 2009-2019 were obtained from Connecticut Office of the Chief Medical Examiner. Using a Bayesian space-time binomial model adjusted for demographic covariates, we estimated spatial and temporal trends in the proportion of fentanyl-involved deaths, as well as overall opioid-involved overdose deaths.
Results During 2009-2019, a total of 6,632 opioid-involved overdose deaths were identified. Among these, 3,234 (49%) were fentanyl-involved. The modeled spatial patterns suggested that opioid-involved deaths in northeastern Connecticut had higher probability of being fentanyl-involved, while New Haven and its neighboring towns and the southwestern region of Connecticut, primarily Greenwich, had a lower risk. Model estimates also suggested the geographic patterns of fentanyl-involved deaths gradually overtook the preceding non-fentanyl opioid overdose deaths. The estimated temporal trend showed the probability of fentanyl involvement increased substantially since 2014.
Conclusion Our findings suggest that geographic variation exists in the probability of fentanyl-involved deaths, and areas at heightened risk are identified. Further studies are warranted to explore the factors contributing to the geographic heterogeneity and continuing dispersion of fentanyl-involved fatal overdoses in Connecticut. Such efforts could inform overdose prevention and intervention strategies.
The three waves of the American opioid crisis (i.e., prescription opioids, heroin, and synthetic opioids other than methadone) continue to be a significant public health emergency in the United States (US).1 Beginning in late 2013, synthetic opioid-involved overdoses (chiefly illicitly produced fentanyl and fentanyl analogs) have dominated the ongoing third wave of the opioid crisis.2 Illicit drug supply is the key driver of the current opioid crisis.3 Illicitly manufactured fentanyl and fentanyl analogs, which are less costly to produce and distribute (but are 50 to 100 times more potent than morphine),4 have increased in the US drug market.5 Reports from the US Drug Enforcement Administration highlighted the dramatic increase in drug seizures that tested positive for fentanyl, from 4,642 in 2014 to 98,954 in 2019 in the US.6,7 The national rate of drug overdose deaths involving fentanyl increased dramatically from 0.5 per 100,000 persons in 2011 to 5.9 per 100,000 in 2016.8 From 2014 to 2015, overdose deaths attributed to synthetic opioids (primarily illicitly manufactured fentanyl) increased by 72% to nearly 10,000.9 Most recently, the rate of drug overdose deaths involving synthetic opioids rose by 27%, from 9.0 per 100,000 in 2017 to 11.4 in 2019.10,11 In 2018, fentanyl and its analogs were involved in approximately two thirds of opioid deaths.12
Abrupt changes in illicit drug supply may have contributed to the heightened risk of opioid overdoses since 2013 and the influx of fentanyl.3 Examining the spatial and temporal trends of fentanyl-involved deaths may help public health researchers and policymakers better understand geographic and temporal variation in fentanyl availability. While several studies have been conducted to investigate the small-area geographic distribution of fentanyl-involved overdose deaths13,14, to our knowledge, no research has explored whether the fentanyl-involved overdose deaths simply take over the rising overall opioid-involved overdoses in the same geographic region or whether fentanyl has established new regions of risk for overdose where, previously, overdoses were less common.
To further inform overdose prevention and intervention efforts, we investigated the geographic and temporal trends of fentanyl-involved overdose deaths at the city/town level in Connecticut, a state that is highly affected by the opioid crisis. We examined whether spatial and temporal distributions of fentanyl-involved overdose deaths covary with overall trends of opioid-involved overdose deaths, and whether the rising tide of fentanyl-involved overdose deaths follows a geographic pattern different from that of the preceding wave of opioid-involved overdose deaths.
METHODS
Study Sample and Data
We obtained data on all opioid-involved overdose deaths from the records of the Connecticut Office of the Chief Medical Examiner (OCME), including the cause and manner of death, toxicological test results in available specimens, demographic information (i.e., age, sex, race/ethnicity), and the addresses of residence, injury, and/or death for each case.15,16
All opioid-involved overdose deaths from 2009 to 2019 in Connecticut that the OCME determined to be of accidental or undetermined manner were included in this analysis. Non-residents in Connecticut were excluded. We dichotomized the opioids involved in the overdose deaths (i.e., fentanyl-involved versus non-fentanyl-involved). Other substances of interest included heroin/morphine, pharmaceutical opioids (i.e., di-hydrocodeine, hydromorphone, hydrocodone, hydromorphone, oxymorphone, and tramadol), methadone/buprenorphine, benzodiazepines, cocaine, ethanol, methamphetamine, xylazine, mitragynine (kratom), and gabapentin. This research project has been reviewed by the Connecticut OCME and deemed not human subjects research by the Yale University Human Investigations Committee.
Geospatial Data
All residential, injury, and death addresses were geocoded within ArcGIS (ESRI, Redlands, CA). Unmatched addresses were then reviewed manually and geocoded using ArcGIS or the R package tidygeocoder17. In cases where decedents were listed as homeless or where no address was recorded, these cases were not geocoded. Decedents with residential addresses outside Connecticut were excluded. Injury addresses were used as primary locales for the space-time model and mapping the spatial trends. If injury addresses were unknown, residential addresses were used instead. If both addresses were unknown, death addresses were used. The geocoded overdose deaths were then assigned to one of the 169 corresponding cities or towns in Connecticut.
Covariates
We used demographic information available for residential population at the city/town level in the Connecticut from the American Community Survey (ACS).18 We linked these with the city/town-level ACS 5-year estimates from 2015-2019 data and 2010-2014 data to the opioid-involved overdose death records in 2015-2019 and in 2009-2014, respectively. The covariates were the total population size, proportion of the population by demographic covariates (i.e., age groups, sex, race, ethnicity and education level), proportion foreign born, the proportion with home ownership, median household income, and poverty rate.
Statistical Analyses
We first described demographic information of the decedents with opioid-involved fatal overdoses in Connecticut, stratified by fentanyl-involved overdose deaths and non-fentanyl overdose deaths. Then we examined spatial and temporal trends for all opioid-involved overdose deaths in Connecticut. We modelled the number of opioid-involved overdose deaths nit in town i during year t as independently and identically Poisson distributed variables with mean λit,
The logarithm of the mean number of opioid-involved overdose deaths (λit) was then modeled as
where xit is the vector of covariates for town i at year t (including the time-varying and space-varying variables from the ACS), and θ is a vector of fixed effect coefficients for xit. In addition, αi is the town-level spatial main effect, γt is the yearly temporal main effect, and δit is the interaction term between space (town level) and time (year). The population size in each city/town was included as an offset log(pop_n) for the Poisson model. The spatial term αi is a random effect that follows the conditional autoregressive model proposed by Besag, York and Mollie.19 The random effect can be further decomposed into two components, an intrinsic conditional autoregressive term that smooths each city/town-level estimate by forming a weighted average with all adjacent jurisdictions, plus a spatially unstructured component that models independent location-specific error and is assumed to be independently, identically, and normally distributed across cities or towns. The temporal trend γt is modeled by the sum of two components, a first-order random walk-correlated time component, and a temporally unstructured component that models independent year-specific error and is independently, identically, and normally distributed across calendar years. The space-time interaction term δit, is modelled as an independent noise term for each city/town and time period and allows for temporal trends in a given city/town to deviate from the overall spatial and temporal trends given by αi and γt, so that local patterns can emerge across time and space. Details of the model specification are described in Appendix S1.
Further, we modelled the probability of an opioid-involved overdose death being fentanyl-involved at the city/town level, using a Bayesian space-time binomial model.20 We considered a binomial model for the number of fentanyl-involved overdose deaths conditional on the total number of opioid-involved overdose deaths in the town areas. Let pit be the probability of an opioid-involved overdose death being fentanyl-involved in town i during year t. We assumed that the number of fentanyl-involved overdose deaths yit in town i during year t is distributed as
and the corresponding likelihood is
A logistic regression was used to model the probability pit, as
where, similar to model for opioid-involved overdose deaths, xit is the vector of covariates for city/town i at year t (including the time-varying and space-varying variables from the ACS), and θ is a vector of fixed effect coefficients for xit. In addition, μi in the model is the city/town-level spatial main effect, φt is the yearly temporal main effect, and σit is the interaction term between space (city/town level) and time (year). In addition, a Poisson model for counts of fentanyl-involved overdose deaths among opioid overdose deaths was examined in a sensitivity analysis.
Finally, we used likelihood ratio tests to examine whether the space-time interaction terms were significantly different from zero in all models. When this term is not significantly different from zero, it suggests that the geographic pattern of fentanyl-involved overdose deaths does not vary significantly over the study period.
To estimate Bayesian model parameters, we employed integrated nested Laplace approximations (INLAs) which approximate the full posterior distribution and are a computationally efficient alternative to Markov Chain Monte Carlo.21 We used the R-INLA package for model fitting.22 Model comparison was performed, and details can be found in Appendix S2. All analyses were performed using R Statistical Software (version 4.0.2; R Foundation for Statistical Computing, Vienna, Austria).
RESULTS
A total of 6,632 opioid-involved overdose deaths were identified by the Connecticut OCME between 2009 and 2019, after excluding 27 deaths with residential addresses outside Connecticut. Among these, 3,234 (49%) were fentanyl-involved, and 3,398 (51%) were non-fentanyl-involved fatalities. The characteristics of these overdose deaths are described in Table 1. People who died of fentanyl-involved overdose deaths were more likely to be male, Black, Hispanic, or involve at least one of the following substances: cocaine, xylazine, gabapentin or mitragynine, compared with those who died of non-fentanyl-involved overdose.
Characteristics of fentanyl-involved overdose deaths among Connecticut residents, 2009-2019
Figure 1 shows the observed temporal trend of fentanyl-involved overdose deaths and non-fentanyl overdose deaths in Connecticut during the study period 2009-2019. Fentanyl-involved overdose deaths increased significantly since 2014. In 2019, 977 (86%) of the 1,138 opioid-involved overdose deaths were fentanyl-involved, compared to 15 (6%) of 266 opioid-involved overdose deaths in 2009 and 36 (8%) of 429 opioid-involved overdose deaths in 2013. The yearly geographic patterns of observed proportion being fentanyl-involved among opioid-involved overdose deaths at the town level are shown in Appendix S3.
Monthly temporal trend fentanyl-involved overdose deaths among Connecticut residents, 2009-2019. This stacked figure shows, as the number of overdose deaths increased, the share of overdose deaths involving fentanyl increased as well.
The marginal posterior distributions of spatial (αi) and temporal (γt) random effects of Bayesian space-time Poisson model for overall opioid-involved overdose deaths are shown in Figure 2. Some towns – Torrington and Southington – were at heightened risk of opioid-involved overdose deaths. The temporal trends showed the risk of opioid-involved overdose deaths among residents in Connecticut increased substantially from 2009 to 2019, although there was a slight decrease in 2018.
Marginal posterior means of the autoregressive spatial (αi) and temporal (γt) random effects (relative risk) from the Bayesian space-time Poisson model for overall opioid-involved overdose deaths. The left panel depicts the spatial patterns of overall opioid-involved overdose deaths when holding temporal terms constant. The right panel depicts the temporal patterns of overall opioid-involved overdose deaths when holding spatial terms constant. The random effects were exponentiated. Larger values of relative risk indicate the higher risk of opioid-involved overdose deaths at town level in Connecticut.
Based on the fitted Bayesian space-time binomial model for the proportion of fentanyl-involved overdose deaths, posterior distributions of the fixed-effect coefficients (θ) of covariates are summarized in Figure 3. Increased proportion of residents aged 55-64 in a city/town was associated with decreased probability being fentanyl-involved, given an opioid-involved overdose death in that town, with an odds ratio (OR) of 0.92 (95% compatibility interval [CI], 0.86, 0.99). In addition, the proportion in poverty was negatively associated with the probability of a fatal overdose being fentanyl-involved (OR=0.95 [95% CI, 0.91, 0.99]).
Summary of the posterior means and 95% compatibility intervals (CIs) of the fixed effect coefficients of covariates. The regression coefficients are exponentiated to represent exponentiated coefficients from binomial logistic regression.
The marginal posterior distributions of spatial (μi) and temporal (φt) random effects are summarized in Figure 4, along with geographic patterns across towns. Some areas had increased probability of fentanyl involvement, and the spatial modeling suggests that even after adjusting for the city/town level demographic covariates, the northeastern region of Connecticut had a higher probability of fentanyl-involved deaths, with Hartford, East Hartford and Manchester at the highest risk. In contrast, New Haven and its surrounding towns and the southwestern Connecticut—primarily Greenwich—had a lower probability of fentanyl-involved deaths. The temporal trends in Figure 4 describe main temporal trends similar to those observed in Figure 1, that is, that the probability of an opioid-involved overdose death being fentanyl-involved increased monotonically since 2014.
Marginal posterior means of the autoregressive spatial (μi) and temporal (φt) random effects from the Bayesian space-time binomial model for fentanyl-involved overdose deaths. The left panel depicts the spatial patterns of fentanyl-involved overdose deaths when holding temporal terms constant. The right panel depicts the temporal patterns of fentanyl-involved overdose deaths when holding spatial terms constant. The random effects were exponentiated. Larger values of exponentiated coefficients indicate the higher probability being fentanyl involved given an opioid overdose death at town level in Connecticut.
The city/town-level predicted probabilities of being a fentanyl-involved overdose death from the Bayesian space-time binomial model are shown in Figure 5. A higher probability being fentanyl-involved, given an opioid-involved overdose death, was first found in several towns in the northeastern region of Connecticut during 2009-2013. Beginning in 2014, the fentanyl “hotspots” started to spread across Connecticut.
The predicted probability of an overdose death being fentanyl-involved from the Bayesian space-time binomial model, stratified by town and calendar year, Connecticut. Color scales differ by year to show the distribution within each year.
A likelihood ratio test shows that the space-time interaction term of the Bayesian space-time binomial model for the proportion of fentanyl-involved opioid overdose deaths was not significantly different from zero, suggesting the geographic patterns of fentanyl-involved overdose deaths gradually overtook the preceding non-fentanyl opioid overdose deaths.
The results of the sensitivity analysis using Poisson model for fentanyl-involved overdose deaths in Connecticut were similar to those using a binomial logistic model and are described in Appendix S4.
DISCUSSION
In the present study we examined the town-level geographic patterns and yearly temporal trends of fentanyl-involved overdose deaths among Connecticut residents during 2009-2019 and evaluated the relationship of fentanyl-involved overdose deaths with overall trends of opioid-involved overdose deaths. To our knowledge, this is the first study to map the spatial and temporal distributions of fentanyl-involved overdose deaths in Connecticut. We identified regions of Connecticut, particularly in the northeastern part of the state, as having relatively high probability of an opioid overdose death being fentanyl-involved. The temporal trends show the fentanyl-involved overdose deaths in Connecticut increased significantly since 2014. We also found that the geographic pattern of fentanyl-involved overdose deaths was relatively constant over time and gradually replaced the preceding wave of opioid-involved overdose deaths.
Our findings show that, compared with non-fentanyl overdose deaths, people who died of fentanyl-involved overdose are more likely to be Black or Hispanic in Connecticut, suggesting that these groups are disproportionately affected by the most recent waves of opioid crisis. These results are consistent with the findings from previous studies that examined opioid-involved overdose deaths in other regions and across the US.23–25 More research is warranted to address the social determinants of fentanyl-involved overdose deaths among Black and Hispanic communities, and locally informed harm reduction efforts and services should be targeted to these populations.
Our results indicated that the northeastern region of Connecticut had higher probability of an opioid overdose death being fentanyl-involved whereas New Haven, its surrounding towns, and the southwestern Connecticut had a lower probability. This suggests that the Connecticut supply of fentanyl and its analogues may have originally entered Connecticut from the north, gradually diffused from the northeastern to southwestern Connecticut, and most recently has dispersed across the state. Our findings also suggest that the fentanyl-involved overdose deaths might simply replace the preceding waves of opioid overdose deaths (e.g., heroin), rather than create new overdose risk following distinct geographic patterns. These findings were in line with the facts that fentanyl and its analogues entered into illicit drug supply as an “adulterant” in powder heroin and counterfeit pills in Northeast.26 The growing presence and strong potency of fentanyl and its analogues then led to the significant increase in opioid overdose deaths across Connecticut in recent years.
Our study is subject to several limitations. First, although we adjusted for city/town-level demographic covariates in our Bayesian space-time models, it is possible that this ecological analysis obscures individual-level relationships between time, place, and overdose risk. We described the geographic patterns and temporal trends of fentanyl-involved overdose deaths at the city/town level only and could not examine in detail the potential social determinants for an opioid overdose death being fentanyl-involved. Second, we primarily used injury addresses for the geographic locales and used residential or death addresses only when injury address was missing, which may lead to geographical misclassification. However, most deaths occurred within person’s own residence so that the injury and residential addresses were the same. Moreover, death addresses were often in the hospital near the injury place. Therefore, we expect that the influence of geographic misclassification on our study findings was minimal. Lastly, while we investigated the spatial and temporal trends of fentanyl-involved overdose deaths, we did not explore the trends of overdose deaths that involved both fentanyl and other substances such as methamphetamine. Further studies are needed to investigate the spatial and temporal trends of polysubstance use with fentanyl in Connecticut.
In conclusion, given the high probability of an overdose death being fentanyl-involved, more research should be devoted to exploring the social determinants and supply-side drivers of fentanyl-involved overdose deaths in order to better inform future harm reduction services and policies to halt the current opioid crisis.
Data Availability
All data produced in the present study are available upon reasonable request to the authors
Acknowledgements
We are also grateful to Dr. James R. Gill and his staff at the Connecticut Office of the Chief Medical Examiner for providing access to data. In addition, we thank Joshua L. Warren, Thomas A. Thornhill, Julia Dennett, and A Ram for helpful comments.
Appendix
Appendix S1 Priors on random effects in Bayesian space-time Poisson model for overall opioid-involved overdose deaths
where 
represents the spatial structured random effect and is modeled under the class of intrinsic Gaussian Markov random fields models. Q denotes the precision matrix (neighboring matrix), and Q− is the generalized inverse of the matrix Q. The marginal variances are 
, which are dependent on the matrix 
is the spatial unstructured random effect and 
is the marginal variance. Gamma priors with small rate parameters are commonly assigned to τu and τv. Here, Gamma(1, 0.0005) is considered.
is first order random walk temporal random effect defined as a random step at each point in time (Δπi). All random steps are independent and identically distributed. 
is the temporal unstructured random effect and 
is the marginal variance. Gamma priors with small rate parameters are commonly assigned to τπ and τρ. Here, Gamma(1, 0.00005) is considered.
Appendix S2: Model Comparison
Appendix S2.1: Model Comparison for overall opioid-involved overdose deaths
For overall opioid-involved overdose deaths, Poisson models were used to estimate the mean number of opioid-involved overdose deaths in each town at each year from 2009-2019. Concerning some towns with 0 opioid-involved overdose deaths. Zero-inflated Poisson models were also examined. Deviance information criterion (DIC) and Watanabe-Akaike information criterion were used to compare Poisson model and zero-inflated Poisson model.
Model Comparison for overall opioid-involved overdose deaths
Based on DIC and WAIC, Poisson model is selected.
Appendix S2.2: Model Comparison for fentanyl-involved overdose deaths
For fentanyl-involved overdose death among all opioid-involved overdose deaths, binomial model and logistic regression were used to estimate the probability being fentanyl-involved given an opioid overdose death in each town at each year from 2009-2019. Several parsimonious models, in addition to the logistic model mentioned in the main text were examined. DIC and WAIC were used for model comparison.
Model Comparison for fentanyl-involved overdose deaths
Based on DIC and WAIC, the binomial model including conditional autoregressive term at town level, first-order random walk structured temporal term for calendar year, space-time interaction term, and demographic covariates, without unstructured temporal term, was selected.
Appendix S3
Spatial patterns of observed proportion being fentanyl-involved among opioid overdose deaths at the city/town level, stratified by year, 2009-2019. Towns with grey color represent no opioid-related deaths.
Appendix S4
Poisson model was used as an alternative to binomial model to estimate the risk of fentanyl-involved overdose deaths among opioid-involved overdose deaths at town level in Connecticut, 2009-2019. Specifically, we modelled the number of fentanyl-involved overdose deaths yit in town i during year t as independently and identically Poisson distributed variables with mean μit,
Then the logarithm of the mean number of fentanyl-involved overdose deaths (μit) is modeled as
where xit is the vector of covariates for town i at year t (including the time-varying and space-varying variables from ACS as aforementioned), and β is a vector of fixed effect coefficients for xit. In addition, αi in the model is town-level spatial main effect, φt is the yearly temporal main effect, and δit is the interaction term between space (town level) and time (year). The number of overall opioid-involved overdose deaths in each town in each year was used as an offset log(od_n).
Specifically, the spatial term αi is a random effect that follows the conditional autoregressive model proposed by Besag, York and Mollie. The random effect can be further decomposed into two components, an intrinsic conditional autoregressive term that smooths each town-level estimate by forming a weighted average with all adjacent census tracts, plus spatially unstructured component that models independent location-specific error and is assumed to be independently identically normally distributed across towns. The temporal trend φt, is modeled by the sum of two components, a first-order random walk-correlated time component, and a temporally unstructured component that models independent year-specific error and is independently identically and normally distributed across calendar years. The space-time interaction term δit, is modelled as an independent noise term for each town and time period, and allows for temporal trends in a given towns to deviate from the overall spatial and temporal trends given by αi and φt, such that local patterns can emerge across time and space. The results were shown in Appendix S4 Figure 1.
Marginal posterior means of the autoregressive spatial (αi) and temporal (φt) random effects from the Bayesian space-time Poisson model for fentanyl-involved overdose deaths. The random effects were exponentiated. Larger values of exponentiated coefficients indicate the higher risk of being fentanyl involved among all opioid overdose death at town level in Connecticut.
Footnotes
Funding sources: National Institutes of Health/National Institute on Drug Abuse DP2-DA049282 and R37DA15612); National Institutes of Health/Eunice Kennedy Shriver National Institute of Child Health and Human Development DP2-HD091799
Conflict of interest: None declared.
Data and code: Data are available on request.
REFERENCES
