This paper is in the following e-collection/theme issue:

Statistical Methods for Surveillance and Population Health (156) Surveillance Systems (462) Clinical Information and Decision Making (1407) Clinical Informatics (1022)

Published on in Vol 9 (2023)

Preprints (earlier versions) of this paper are available at https://preprints.jmir.org/preprint/41450, first published .
Small Area Forecasting of Opioid-Related Mortality: Bayesian Spatiotemporal Dynamic Modeling Approach

Small Area Forecasting of Opioid-Related Mortality: Bayesian Spatiotemporal Dynamic Modeling Approach

Small Area Forecasting of Opioid-Related Mortality: Bayesian Spatiotemporal Dynamic Modeling Approach

Authors of this article:

Cici Bauer1 Author Orcid Image ;   Kehe Zhang1 Author Orcid Image ;   Wenjun Li2 Author Orcid Image ;   Dana Bernson3 Author Orcid Image ;   Olaf Dammann4, 5 Author Orcid Image ;   Marc R LaRochelle6, 7 Author Orcid Image ;   Thomas J Stopka4, 8, 9 Author Orcid Image
Article Authors Cited by (3) Tweetations (2) Metrics

    Original Paper

    1Department of Biostatistics and Data Science, School of Public Health, The University of Texas Health Science Center at Houston, Houston, TX, United States

    2Department of Public Health, University of Massachusetts Lowell, Lowell, MA, United States

    3Office of Population Health, Department of Public Health, The Commonwealth of Massachusetts, Boston, MA, United States

    4Department of Public Health and Community Medicine, Tufts University School of Medicine, Boston, MA, United States

    5Department of Gynecology and Obstetrics, Hannover Medical School, Hannover, Germany

    6Clinical Addiction Research and Education Unit, Section of General Internal Medicine, Department of Medicine, Boston University School of Medicine, Boston, MA, United States

    7Grayken Center for Addiction, Boston Medical Center, Boston, MA, United States

    8Department of Urban and Environmental Policy and Planning, Tufts University, Medford, MA, United States

    9Department of Community Health, Tufts University, Medford, MA, United States

    Corresponding Author:

    Cici Bauer, PhD

    Department of Biostatistics and Data Science

    School of Public Health

    The University of Texas Health Science Center at Houston

    1200 Pressler Street, Room E819

    Houston, TX, 77030

    United States

    Phone: 1 713 500 9581

    Email: cici.x.bauer@uth.tmc.edu


    Background: Opioid-related overdose mortality has remained at crisis levels across the United States, increasing 5-fold and worsened during the COVID-19 pandemic. The ability to provide forecasts of opioid-related mortality at granular geographical and temporal scales may help guide preemptive public health responses. Current forecasting models focus on prediction on a large geographical scale, such as states or counties, lacking the spatial granularity that local public health officials desire to guide policy decisions and resource allocation.

    Objective: The overarching objective of our study was to develop Bayesian spatiotemporal dynamic models to predict opioid-related mortality counts and rates at temporally and geographically granular scales (ie, ZIP Code Tabulation Areas [ZCTAs]) for Massachusetts.

    Methods: We obtained decedent data from the Massachusetts Registry of Vital Records and Statistics for 2005 through 2019. We developed Bayesian spatiotemporal dynamic models to predict opioid-related mortality across Massachusetts’ 537 ZCTAs. We evaluated the prediction performance of our models using the one-year ahead approach. We investigated the potential improvement of prediction accuracy by incorporating ZCTA-level demographic and socioeconomic determinants. We identified ZCTAs with the highest predicted opioid-related mortality in terms of rates and counts and stratified them by rural and urban areas.

    Results: Bayesian dynamic models with the full spatial and temporal dependency performed best. Inclusion of the ZCTA-level demographic and socioeconomic variables as predictors improved the prediction accuracy, but only in the model that did not account for the neighborhood-level spatial dependency of the ZCTAs. Predictions were better for urban areas than for rural areas, which were more sparsely populated. Using the best performing model and the Massachusetts opioid-related mortality data from 2005 through 2019, our models suggested a stabilizing pattern in opioid-related overdose mortality in 2020 and 2021 if there were no disruptive changes to the trends observed for 2005-2019.

    Conclusions: Our Bayesian spatiotemporal models focused on opioid-related overdose mortality data facilitated prediction approaches that can inform preemptive public health decision-making and resource allocation. While sparse data from rural and less populated locales typically pose special challenges in small area predictions, our dynamic Bayesian models, which maximized information borrowing across geographic areas and time points, were used to provide more accurate predictions for small areas. Such approaches can be replicated in other jurisdictions and at varying temporal and geographical levels. We encourage the formation of a modeling consortium for fatal opioid-related overdose predictions, where different modeling techniques could be ensembled to inform public health policy.

    JMIR Public Health Surveill 2023;9:e41450

    doi:10.2196/41450

    Keywords

    opioid-related mortality; small area estimation; spatiotemporal models; Bayesian; forecasting 


    Opioid-related overdoses continue to be at crisis levels in communities across the United States, with more than 75,673 fatal overdoses in the 12-month period ending in April 2021 [1-7], worsening during the COVID-19 pandemic [8]. Opioid-related deaths increased more than 5-fold in Massachusetts between 2000 and 2016, with more than 2000 per year from 2016 to 2021 [9,10]. Fatal opioid-related overdose rates above the national level have been ubiquitous across communities in Massachusetts [11]. Despite this crisis, public health responses to the opioid overdose epidemic have been limited by an inability to rapidly identify current fatal overdose patterns, predict future local clusters, and evaluate the effectiveness of interventions.

    Identification and prediction of local fatal opioid overdoses require a comprehensive and high-quality surveillance system that provides data sources to capture granular geographic information of the fatal opioid overdose cases across space and time. Ideally, additional information such as individual demographics, past medical history (particularly mental health history), local drug supply, and other risk factors should be included as well to enhance our understanding of the opioid crisis. Many states have established surveillance systems to monitor opioid-related morbidity and mortality to inform planning and evaluate control efforts. Some surveillance systems include unlinked individual data sources for different opioid-related reporting (eg, vital records and prescription drug monitoring programs), and some aim to provide an individually linked database across various data sources [12]. While the latter provides much enhanced data capacity for a wide range of opioid-related research, creating such linked databases takes substantial time and financial resources.

    Accurate identification and prediction of fatal opioid-related overdose trends also requires sophisticated spatial and predictive analytical approaches, as noted in a recent review of methodological approaches for the prediction of opioid use–related epidemics in the United States [13]. Most published literature has focused on the identification of community- or neighborhood-level risk factors for opioid-related overdose. For example, Bozorgi et al [14] explored different machine learning and spatial analytical approaches to identify leading contextual risk factors for drug overdose at the block group level in South Carolina. A similar study by Schell et al [15] also focused on identifying new neighborhood-level predictors of opioid-related overdose deaths in Rhode Island at the census block group level, using least absolute shrinkage and selection operator and random forest algorithms. Abell-Hart et al [16] identified counties with a high number of underserved opioid overdose patients in New York state. Basak et al [17] detected spatiotemporal hot spots (ie, counties) at high risk of prescription opioid misuse and overdose using regression models that relied on Medicare claims data in Virginia, North Carolina, and West Virginia.

    Statistical models or machine learning techniques proposed specifically for prediction purposes are limited, with most previous studies performed at the US county or state level [18,19], lacking the spatial granularity needed to guide local public health departments for preemptive actions. Bayesian spatiotemporal models have received substantial attention in the past several years in opioid-related research, given their ability to include temporal and spatial correlations and improved precision in small area estimation [20-22]. Sumetsky et al [20] for instance, developed a Bayesian logistic growth model for opioid overdose mortality predictions for 146 counties in North and South Carolina.

    In this study, we developed and validated several Bayesian spatiotemporal dynamic predictive models designed for small area forecasting of opioid-related overdose mortality at the ZIP Code Tabulation Area (ZCTA) level. We investigated the benefits of including various area-level demographic and socioeconomic factors in improving the prediction. The predictive performance was assessed using opioid-related overdose mortality data from Massachusetts between 2005 and 2019. We identified the top ZCTAs with the highest predicted fatal opioid overdose rates or counts and by urban and rural areas. Our prediction results can help inform local public health departments’ planning and targeting of resource allocations.


    Opioid-Related Mortality Data

    Opioid-related mortality data were obtained from the Massachusetts Department of Public Health’s (MDPH’s) Registry of Vital Records and Statistics (RVRS) for 2005 through 2019. International Classification of Disease, Tenth Revision, codes for mortality were used to select from the underlying cause of death field by RVRS staff and identify poisonings or overdoses: X40-X44, X60-X64, X85, and Y10-Y14. All multiple underlying cause of death fields were then used to define opioid-related death with T40.0, T40.1, T40.2, T40.3, T40.4, and T40.6. We excluded opioid-related deaths for individuals younger than 20 years, since fatal opioid-related death was rare in this age group (n=191). This age grouping also allows for comparison of opioid-related overdose rates with other drug overdose studies [23]. We used the ZCTA of the reported residential address for decedents at the time of the fatal overdose and excluded those either with missing addresses or those with addresses outside of Massachusetts’ boundaries (n=642). The final analytic sample included 16,377 fatal opioid-related overdoses from 2005 to 2019. The flowchart for deriving the analytic samples used for this analysis, summarized by each year, is presented in Figure S1 in Multimedia Appendix 1. We also developed a web-based dashboard to present the map and summaries of these data [24]. Requests for the opioid-related decedent data should be made directly to the MDPH’s RVRS.

    ZCTA Demographic and Socioeconomic Factors

    ZCTA-level demographic and socioeconomic variables were obtained from 2005 to 2019 from the American Community Survey 5-year estimates (2005-2009, 2010-2014, and 2015-2019, respectively) [25]. We used the total population counts and characteristics for people aged 20 years or older at the ZCTA level to account for differences in population size. The ZCTA-level demographic variables included race or ethnicity (proportions of White and Hispanic individuals), education (proportion of individuals with a bachelor’s degree), employment (proportion of unemployed individuals), poverty level (proportion of individuals living under the federal poverty level), income (per capita income in US $), living conditions (proportion of renters), transportation (proportion of individuals without vehicles), and English speaking (proportions of individuals with limited English). A detailed summary of these covariates, for all ZCTAs and those stratified by urban or rural areas, is presented in Table S1 in Multimedia Appendix 1. These covariates were considered from a larger set of community-level covariates and were selected on the basis of the literature and our prior research. Following calculation of descriptive statistics and bivariate analyses, we further narrowed the list to a subset of population-level measures at the ZCTA level. The final set of covariates included in this analysis was chosen after fitting a simple Poisson regression model and selecting those with a significant association with opioid-related overdose mortality (see Figure S5 in Multimedia Appendix 1).

    Urban and Rural ZCTAs

    The Massachusetts State Office of Rural Health uses a composite scoring system to categorize rurality for each ZCTA. Among the 537 ZCTAs included in this analysis, 390 (73%) were classified as urban areas, and the rest as rural areas. A map of the ZCTAs, color-coded by their urban or rural status, is provided in Figure S2 in Multimedia Appendix 1.

    Statistical Models

    Our choice of dynamic spatiotemporal model was motivated by a previous study of drug overdose rates across US counties, which showed that the strongest predictor of overdose rates was overdose rates in nearby counties from the previous year. We considered the following Bayesian dynamic spatiotemporal models [26,27]. Yit denotes the number of opioid-related fatal overdose counts in ZCTA i (i=1, …, I=537) and year t (t=2005, …, t=2019). We assumed a Poisson distribution with mortality rate μit; that is, Yit|μit ~ Poisson(Nitμit), with Nit representing the population size as an offset. The mortality rate μit was then decomposed using the following models:

    Model 1: log (μit) = αi + S(t) + ηt, ηt ~ AR(1),

    Model 2: log (μit) = αi + S(t) + ωit, ωit ~ GMRF(τΣAR(1) ⊗︀ΣI), and

    Model 3: log (μit) = α + S(t) + ωit, ωit ~ GMRF(τΣAR(1) ⊗︀ΣCAR)

    In all 3 models, S(t) captured the overall (or marginal) temporal trend and could be modeled via linear or quadratic terms or spline functions. The major differences among the 3 models were how random effects (η or ω) were modeled in the spatiotemporal interactions. In model 1, αi was the random effect accounting for the difference of the opioid-related mortality at the ZCTA level, and ηt was a first-order autogressive latent effect AR(1) for the temporal trend. Essentially, model 1 regressed the opioid-related mortality rate at time t, on the log scale, over the previous time point t–1; however, instead of running one model for each ZCTA, the Bayesian framework allowed pooling of the ZCTAs to improve prediction. In models 2 and 3, α was the overall intercept, and the ZCTA-level differences of the opioid-related mortality rates were absorbed in the spatiotemporal interaction term ωit, which now accounted for the spatial dependency between ZCTAs. Models 2 and 3 were inseparable spatiotemporal models, where the interaction term ωit had a Gaussian Markov Random Field (GMRF) with mean 0 and covariance matrix τΣT⊗︀ΣS [28]. The interaction term was the Kronecker product of the temporal structure ΣT and the spatial structure ΣS, and τ was the precision parameter (ie, reciprocal of the variance parameter). We again used the AR(1) temporal structure in models 2 and 3 but considered different structures for the spatial dependency. In model 2, the spatial structure was an I×I identity matrix. Although the identity matrix may appear to assume independence of the areas, the hierarchical structure of the model imposed borrowing information from the neighboring areas. In model 3, we assumed a conditional autoregressive model (CAR) [29] as the spatial structure, which assumed that ZCTAs that were geographically adjacent were more similar than those that were far away. We note that other spatial or temporal structures (eg, random walk models) can also be considered in the proposed model framework.

    Models 1, 2, and 3 were referred to as the respective base models. We then added ZCTA-level demographic and socioeconomic variables xit and the urban or rural indicator variable described above, to each of the base models to assess potential improved prediction performance. All models were fit in a Bayesian framework, with the priors chosen as the noninformative priors commonly used in spatiotemporal models [28]. In our application, it turned out that the interaction term was sufficient to capture the spatial-temporal variation, so we dropped the marginal temporal term. We used the posterior median and 90% credible intervals (CrIs) for inference. The 90% CrI was preferred over the commonly used 95% CrI, as the former was better suited for prediction purposes to avoid overly wide uncertainty. All analyses were performed in RStudio [30] and R package integrated nested Laplace approximation [31].

    To assess the predictive performance of the models, we used a one-year ahead approach: assuming we observed data from 2005 to year t, we predicted the opioid-related mortality rates and counts for year t+1. We started by considering the observed data from 2005 to 2015, and the prediction was carried out for 2016. Then we considered observed data from 2005 to 2016 and predicted counts and rates were calculated for 2017. The one-year ahead prediction procedures were carried out for each year from 2016 to 2019, and the predictions were obtained for each ZCTA, along with the prediction uncertainty (posterior 90% CrI). We used the following metrics to compare the predictive performance across different models: mean absolute error (MAE), root mean square error (RMSE), and the rank difference (RD). MAE was defined as the average absolute difference between the predictive and the observed values; that is, , for count or rate Yit. RMSE was defined as the average squared difference between the predictive and the observed values; that is, . RMSE has the advantage that it was on the same scale as the outcome variable and was thus easy to interpret. For example, if the prediction was for the ZCTA-level mortality count (ie, how many people would die from fatal opioid-related overdose), then an RMSE of 10 would roughly mean that on average the prediction was off by 10 mortality cases. Smaller MAE and RMSE would indicate better performance. The RD was motivated by the correct classification of opioid mortality by group membership, where we divided the ZCTAs into quintiles Q1 to Q5, based on the observed fatal opioid-related overdose rates and counts, and then compared to the quintiles for the predicted rates and counts. The RD was then calculated by the percentage of times when the classification was correct. Higher values of RD indicated better performance.

    Ethical Considerations

    This study was reviewed by the Health Sciences institutional review board at Tufts University and was designated as exempt from ethics approval (IRB reference number: 13288).


    The number of fatal opioid-related overdoses among people aged 20 years or older increased from 557 in 2005 to 1912 in 2019, corresponding to an increase in rate from 7.95 in 2005 to 34.4 in 2019 per 100,000 population. Out of a total of 16,377 opioid-related fatalities between 2005 and 2019, overall, 91% occurred in urban ZCTAs. The observed rates at the ZCTA level were highly variable (Figure S3 in Multimedia Appendix 1), ranging from 0 to 1316 per 100,000 population. The extreme values in the observed rates tended to occur in areas with small populations. The high instability of the observed rates, commonly known as the small area estimation problem, poses a special challenge in developing accurate prediction models, which our Bayesian spatiotemporal models helped address.

    Results of the predictive performance among the proposed candidate models, using the one-year ahead approach, are presented in Table 1 (for the ZCTA-level fatal opioid-related overdose count) and Table 2 (for the ZCTA-level fatal opioid-related overdose rate). Prediction errors were summarized for all ZCTAs and stratified by rural and urban ZCTAs. Overall, model 3 with the inseparable spatiotemporal interaction term of AR(1) and CAR structures performed best for fatal opioid-related overdose predictions, as indicated by smaller MAE and RMSE values. The RMSE showed that the predicted opioid-related death counts, on average, were off only by 2 counts per area. As expected, predictions were better for urban areas than for rural areas, since most of the deaths occurred in urban areas. Addition of demographic and socioeconomic variables generally improved the prediction performance, particularly for model 1; however, the improvement was attenuated in models 2 and 3 where we assumed inseparable spatiotemporal models. The smaller improvement in model 3 was likely because the CAR model already captured the spatial patterns tied to the demographic and socioeconomic factors and hence served as a surrogate of the demographic and socioeconomic variables for predictions. The best performing model using the rank difference showed a less clear pattern, as the values were very similar across the different models. This was particularly true when focusing on the rural ZCTAs, likely due to their small predicted rates or counts and hence very narrow ranges for each quintile.

    Table 1. Predictive performance assessment of the ZIP Code Tabulation Area (ZCTA) level opioid-related overdose death count, using the root mean and rank difference, and one-year ahead prediction starting with 2016a. The smallest root mean square error, mean absolute difference, and highest rank difference for each row are depicted in italics, indicating the best performing models.
    Count yearZCTA typeModel 1Model 2Model 3


    		Base ModelAdd SDOHbBase ModelAdd SDOHBase ModelAdd SDOH
    Root mean square error

    		2016All2.562.632.982.852.352.40

    		2017All2.061.982.292.202.041.94

    		2018All2.532.472.722.562.472.44

    		2019All2.382.332.572.432.192.18

    		2016Urban2.902.983.383.222.652.70

    		2017Urban2.342.252.622.502.302.20

    		2018Urban2.832.753.042.842.762.71

    		2019Urban2.652.592.872.702.462.45

    		2016Rural1.261.321.421.431.201.28

    		2017Rural0.980.971.031.031.030.98

    		2018Rural1.451.451.621.621.431.48

    		2019Rural1.421.411.461.471.151.16
    Mean absolute difference

    		2016All0.920.981.431.320.670.74

    		2017All0.010.020.610.500.260.20

    		2018All0.550.560.820.730.320.37

    		2019All0.170.190.580.510.000.06

    		2016Urban1.151.211.771.600.850.91

    		2017Urban0.040.060.790.630.290.22

    		2018Urban0.610.600.910.780.300.35

    		2019Urban0.130.140.640.520.030.04

    		2016Rural0.300.360.520.550.220.28

    		2017Rural0.140.080.130.160.190.13

    		2018Rural0.400.440.580.600.380.43

    		2019Rural0.290.330.430.480.060.12
    Rank difference

    		2016All0.500.510.500.490.520.53

    		2017All0.500.500.490.480.510.52

    		2018All0.500.510.500.520.490.52

    		2019All0.510.510.490.500.500.51

    		2016Urban0.540.540.550.540.570.55

    		2017Urban0.530.530.510.520.520.52

    		2018Urban0.540.550.550.560.540.54

    		2019Urban0.580.600.540.570.610.62

    		2016Rural0.400.390.400.380.370.35

    		2017Rural0.430.420.390.410.410.40

    		2018Rural0.380.370.340.360.370.37

    		2019Rural0.340.370.330.340.370.38

    aFor example, using data from 2005 to 2015, we predicted the fatal opioid-related overdose count in year 2016 for each of the 537 ZIP Code Tabulation Areas (ZCTAs) in Massachusetts, and compared them to the observed data to assess performance. Similarly, for 2017, we used data from 2006 to 2016, and predicted for 2017. The prediction error was assessed for counts for all ZCTAs and stratified by rural and urban status. Base model refers to Bayesian dynamic spatiotemporal models without any covariates; model results with added covariates of area-level contextual factors are included in the column “Add SDOH.”

    bSDOH: Social determinants of health.

    Since our 1-year prediction assessments suggested that the best performing model was model 3 with the included demographic and socioeconomic variables, we used this model for 2020 and 2021 predictions based on the opioid-related mortality data for 2005-2019. We assumed that ZCTA-level population size and demographic and socioeconomic variables had the same values in 2020 and 2021 as those in 2019. We carried out the prediction for each ZCTA and at the state level. Figure 1 presents the fitted temporal trends of opioid-related mortality in Massachusetts between 2005 and 2019, with predictions carried out for 2020 and 2021 (the gray shaded area), both as the rate per 100,000 population (panel A) and total counts (panel B). Each line represents an individual ZCTA, color-coded by urban or rural status. Figure 2 presents the maps of the predicted opioid-related mortality rates and counts for 2020 and 2021. At the ZCTA-level, the predicted opioid-related mortality rate ranged from 6.11 to 162 per 100,000 population for 2020, and 6.05 to 158 per 100,000 population in 2021. Note that because of the “smoothing” effect in Bayesian models, the estimated or predicted rate could be very close to 0 but would never be exactly 0. Therefore, for ZCTAs that may have zero observed deaths, the predicted rate would not be zero but would be skewed toward the average rate. In addition, any ZCTAs that may have low observed rates but were “surrounded” by ZCTAs with high rates would have higher predicted rates, reflecting how the spatial models borrow information from neighboring ZCTAs. Because of the large variation in population size by ZCTA, the spatial patterns seen in the maps in Figure 2 for predicted rate and counts are quite different, where the latter are generally centered around highly populated areas. The predicted count ranges from close to 0 (after rounding up to integers) to 30 (90% CrI 21-41) for 2020, and 29 (90% CrI 19-44) for 2021. We also identified the top 5 ZCTAs with the highest predicted counts, by urban or rural classifications, and presented the prediction results in Figure 3. The urban ZCTAs were located within cities with high fatal opioid-related overdose risks: Lawrence, Lynn, Quincy, Brockton, and New Bedford (Figure 3C). The rural ZCTAs were located in municipalities that were known to have high risks for fatal opioid-related overdoses, including Pittsfield, North Adams, Greenfield, Westfield, and Billerica (Figure 3D).

    At the state level, the predicted opioid-related mortality rate for the population aged 20 years and older was 35.79 per 100,000 population (90% CrI 29.4-43.4) for 2020, and 35.81 (90% CrI 26.9-46.0) for 2021. For urban areas, the predicted rate was 36.0 (90% CrI 29.2-43.7) in 2020, and 35.5 (90% CrI 26.8-46.1) in 2021. For rural areas, the predicted rate was 33.9 (90% CrI 26.6-43.3) in 2020, and 34.6 (90% CrI 25.2-53.4) in 2021. The CrIs for the prediction were wider in 2021 than in 2020 and in rural areas than in urban areas. This was expected as the further ahead we attempted to predict outcomes, the less certainty we would have. Our prediction suggested that the opioid-related overdose death counts for 2020 would be 1887 (90% CrI 1549-2285) for the whole state, with 1682 (90% CrI 1364-2041) in urban ZCTAs, and 203 (90% CrI 160-260) in rural ZCTAs, and those for 2021 would be 1888 (90% CrI 1419-2426) for the whole state with 1656 (90% CrI 1253-2149) in urban ZCTAs, and 208 (90% CrI 151-320) in rural ZCTAs.

    Table 2. Predictive performance assessment of the ZIP Code Tabulation Area (ZCTA) level opioid-related overdose death rate, using the root mean and rank difference, and one-year ahead prediction starting with 2016a. The smallest root mean square error, mean absolute difference, and highest rank difference for each row are depicted in italics, indicating the best performing models.
    Rate yearZCTA typeModel 1Model 2Model 3


    		Base ModelAdd SDOHbBase ModelAdd SDOHBase ModelAdd SDOH
    Root mean square error

    		2016All41.1840.2742.8441.9440.2838.95

    		2017All29.2130.1829.0130.2428.6329.82

    		2018All40.5440.2242.6341.8340.4740.65

    		2019All74.8075.6475.7576.0373.2675.39

    		2016Urban44.5343.3946.5245.4543.5741.65

    		2017Urban30.3131.7630.7532.3029.6331.13

    		2018Urban38.3037.6740.8239.7138.1738.24

    		2019Urban83.4884.4284.4884.7481.5984.13

    		2016Rural30.5030.4530.9830.7429.7930.62

    		2017Rural26.0825.5323.7623.9125.8026.02

    		2018Rural46.0046.3347.1047.0346.0446.47

    		2019Rural44.0544.4944.9045.3143.9944.38
    Mean absolute difference

    		2016All7.297.2713.7312.444.914.95

    		2017All3.704.014.312.815.875.83

    		2018All5.895.3610.599.264.104.05

    		2019All4.474.0110.078.831.601.57

    		2016Urban10.029.4116.3714.417.257.01

    		2017Urban2.113.095.693.565.095.24

    		2018Urban5.724.6110.058.153.012.83

    		2019Urban4.603.6410.028.242.071.91

    		2016Rural0.001.576.687.201.310.55

    		2017Rural7.946.440.640.847.937.40

    		2018Rural6.347.3912.0412.236.987.29

    		2019Rural4.135.0210.2210.400.350.68
    Rank difference

    		2016All0.300.340.310.330.310.32

    		2017All0.320.360.320.370.340.35

    		2018All0.290.310.310.300.310.34

    		2019All0.330.340.340.350.340.32

    		2016Urban0.320.380.310.370.330.38

    		2017Urban0.340.380.350.350.370.38

    		2018Urban0.320.340.300.330.330.34

    		2019Urban0.360.360.360.370.380.36

    		2016Rural0.280.240.300.240.290.26

    		2017Rural0.190.210.300.270.220.23

    		2018Rural0.210.260.240.240.260.26

    		2019Rural0.280.250.250.260.270.28

    aFor example, using data from 2005 to 2015, we predicted the fatal opioid-related overdose rate in 2016 for each of the 537 ZIP Code Tabulation Areas (ZCTAs) in Massachusetts, and compared the predicted rates to the observed rates to assess prediction performance. Similarly, for 2017, we used data from 2006 to 2016 and predicted for 2017. The prediction error was assessed for rates for all ZCTAs and stratified by rural and urban status. Base model refers to Bayesian dynamic spatiotemporal models without any covariates; model results with added covariates of area-level contextual factors are included in the column labeled “Add SDOH.”

    bSDOH: Social determinants of health.

    ‎
    Figure 1. Fitted temporal trends for opioid-related overdose (OD) mortality in Massachusetts between 2005 and 2019, with predictions for 2020 and 2021 (in gray shade). Each line represents a ZIP Code Tabulation Area (ZCTA) color-coded by its urban or rural status. Opioid-related mortality data were obtained from the Massachusetts Registry for Vital Records and Statistics, and predictions were made on the fatal opioid-related overdose rates per 100,000 population (panel A) and count (panel B). These results are from the Bayesian dynamic spatiotemporal Model 3 with ZCTA level demographic and socioeconomic variables.
    ‎
    Figure 2. Maps of the predicted ZIP Code Tabulation Area (ZCTA)–level fatal opioid-related overdose (OD) rates (per 100,000 population) and counts for 2020 (A and B) and 2021 (C and D). Data were obtained from Massachusetts Registry for Vital Records and Statistics for 2005 to 2019 and used to predict for 2020 and 2021. These predictions were obtained from the proposed Bayesian dynamic spatiotemporal Model 3 with ZCTA level demographic and socioeconomic variables. In each panel, the embedded histograms present the distribution of the predicted fatal opioid-related overdose rates or counts for that year.
    ‎
    Figure 3. Five selected ZIP Code Tabulation Areas (ZCTAs) with the highest predicted fatal opioid-related overdose (OD) rates or counts for urban (A, C) and rural (B, D) areas. The vertical bars in panels A and B present the 90% posterior credible intervals (CrIs) of the predictions for years 2020 and 2021. Maps (C, D) present the corresponding locations of the identified ZCTAs. Predictions were performed using the proposed Bayesian dynamic spatiotemporal model 3 with ZCTA level demographic and socioeconomic variables and fatal opioid-related overdose data from 2005 to 2019 in Massachusetts.

    In this analysis, we developed and compared several Bayesian spatiotemporal dynamic models for predicting small area opioid-related mortality. The prediction performance was evaluated with data from the MDPH RVRS for 2005 through 2019 using the one-year ahead approach, with and without ZCTA-level demographic and socioeconomic variables. Using data from 2005 through 2019 and the best performing model, we also predicted the fatal opioid-related death rates and counts for 2020 and 2021, respectively, along with uncertainty assessments. We also identified the ZCTAs—by urban and rural status—that had the highest predicted opioid-related mortality, both by rates and by counts.

    Our prediction showed that, if there was no interruptive change to the trends observed for 2005-2019, we would observe a stabling trend of the fatal opioid-related overdose in Massachusetts for 2020 and 2021. The stabling trend would be applicable to the entire state (ie, the state-level total fatal opioid-related overdose) as well as for most ZCTAs in Massachusetts. Our prediction also identified the ZCTAs deviating from this stabilizing trend, with continually increasing rates above the state average. Identifying ZCTAs with the predicted high risks allows for the possibility of preemptive and geo-targeted public health interventions. The prediction models presented here allowed for a more granular depiction of existing and expected trends in opioid-related overdose deaths over spatially granular units, which differed from other existing models developed for predication at larger geographical scales (eg, state or county). Such predictions are instrumental for state and local public health departments’ planning to identify and address potential service gaps for deploying harm reduction and treatment interventions. Local data are largely limited both in type and quantity; drug seizure data, for example, are often not available at more granular levels than the state, limiting input information in predictive modeling. However, up-to-date data at the local scale are instrumental to developing the prediction models and engaging communities in designing and implementing data-driven responses to reduce opioid-related harms [32].

    At the time of this analysis, we only obtained the fatal opioid-related overdose data for Massachusetts from 2005 to 2019 and partial data in 2020 (incomplete data with only the first 6 months). We designed our analysis and prediction evaluation on the basis of the data available to us, but we lacked the official statistics for 2020 and 2021 to assess the prediction accuracy for those years. The latest Massachusetts opioid-related overdose data brief released in December 2022 [33] reported the confirmed opioid-related overdose deaths among Massachusetts residents for 2020 and 2021. Although the data brief reported the total number and rate including all ages, which was different from the adult population we focused on, it was clear that we underpredicted the total number of fatal opioid-related overdose in Massachusetts for 2020 and 2021. The data brief also noted a 9% increase in the opioid-related overdose death rate in Massachusetts for 2021 over 2020. However, since the brief only reported the total opioid overdose death counts at the state level and not by ZCTAs, we were not able to assess the underprediction at a more granular spatial level. The sharp increase in the fatal opioid-related overdose was likely due to the substantive impact of the COVID-19 pandemic [8,34], which disrupted any established trends prior to 2020 that informed our prediction modeling. The time lag in the drug overdose database had been identified as a main barrier to developing fatal overdose predictive models, as reported by Borquez and Martin [35]. Our future work, in collaboration with the MDPH, plans to use the statewide Public Health Data warehouse [36] to obtain timelier and rich data sources with local information in order to improve predictive performance in small areas.

    It is important to consider the limitations inherent in the highlighted analyses, and recommendations that should be considered when developing future opioid-related overdose prediction models. First, the spatial units of ZCTAs used in our study may not be the ideal spatial unit for prediction, as the population size within ZCTAs vary substantially, compared to other spatial units of analysis (eg, census tracts) where, by design, the population sizes are much more homogenous. In addition, ZCTAs may not be sufficiently granular, from a spatial perspective, to identify local hot spots. However, they are extensively used in spatial analysis as they are most readily available in many aggregated data sources (eg, surveillance or insurance claims databases). ZCTAs provide useful geospatial information for analyses while often also satisfying data privacy concerns. Second, a better understanding of the contributing factors to local opioid overdose trends and patterns, with reliable measures representing such factors, would clearly improve prediction power and accuracy. For example, toxicology data would have helped us to include the appearance of fentanyl in the local drug supply, to aid the prediction in the shifting “waves” in opioid-related fatal overdose [37]. Third, we used the ZCTA of the decedents’ residences, rather than the injury addresses or the locations where the fatal opioid overdoses were recorded. Injury data generally have a high level of missingness (~50% in MA), and the recorded death location, if not at the decedent’s place of residence, is often recorded at a hospital, even though the injury (ie, overdose event) typically occurred elsewhere. Finally, and perhaps most challenging, is to incorporate the impact of emerging phenomena such as the COVID-19 pandemic into predictions of opioid-related overdose trends, which was not possible for our analyses given the data availability at the time of analyses. This task requires the real-time data inputs and requires a joint effort among researchers and practitioners from multiple agencies, institutions, and sectors. We have seen a lot of progress made on this front during the COVID-19 pandemic, and we hope to see more progress in the drug overdose research in the future.

    Despite the abovementioned data source and methodological limitations, our models showed promise in providing reasonable 1-year forecasts of opioid-related mortality in MA with geographic granularity using existing data, as the short-term point estimates for the number of overdoses tended to be close to the true value. Our Bayesian spatiotemporal models further demonstrated the advantages of incorporating inseparable spatiotemporal dependencies over the simpler regression models without such dependence. Since the assumed spatial dependency structure captured the spatial patterns tied to many demographic and socioeconomic factors, such models do not rely on the knowledge of the future measures of these factors in predicting opioid-related mortality. Prediction is a challenging problem in general, and though many models have been developed in various contexts, it is almost impossible to find one single model or approach that universally performs best [38]. Although our analysis could not investigate all possible predictive models for fatal opioid-related overdoses, we provided a novel approach to forecasting overdose events for small geographic areas. Compared to other predictive approaches, Bayesian models provide a natural framework where the prediction can be conveniently included in the model fitting process, by treating predictions as missing values. Our findings demonstrated the utility of sophisticated Bayesian spatiotemporal dynamic models in supporting state and local opioid surveillance and the ability to provide prediction at a granular geographic level, offering a unique opportunity for preemptive public health and policy interventions, replacing reactionary public health responses. Echoing Borquez and Martin [35], we encourage the field to consider a modeling consortium for opioid-related prediction models, where different modeling techniques could be ensembled.

    Acknowledgments

    This research was partially supported by funding from the Centers for Disease Control and Prevention (Assessing High Risk Opioid Prescribers and Fatal and Non-Fatal Opioid Overdose grant INTF7311HH2 500224100 awarded to the principal investigator TJS), and the National Institute Drug Abuse (grant R01DA054267 awarded to CB and TJS).

    Authors' Contributions

    CB and TJS conceived and designed the analysis. TJS contributed to the acquisition of the data. CB and KZ contributed to data processing and curation, conducted the data analysis, and drafted the initial manuscript. WL, MRL, DB, and TJS contributed to critical revision of the manuscript. All authors contributed to the interpretation of the findings, and reviewed and approved the final manuscript.

    Conflicts of Interest

    None declared.

    ‎
    Multimedia Appendix 1

    Trends of opioid-related deaths in Massachusetts.

    PDF File (Adobe PDF File), 9073 KB
    1. Ciccarone D. Fentanyl in the US heroin supply: A rapidly changing risk environment. Int J Drug Policy 2017 Aug;46:107-111 [FREE Full text] [CrossRef] [Medline]
    2. Ciccarone D. The triple wave epidemic: Supply and demand drivers of the US opioid overdose crisis. International Journal of Drug Policy 2019 Sep;71:183-188 [FREE Full text] [CrossRef] [Medline]
    3. Understanding drug overdoses and deaths. Centers for Disease Control and Prevention.   URL: https://www.cdc.gov/drugoverdose/epidemic/index.html [accessed 2023-01-19]
    4. Somerville NJ, O'Donnell J, Gladden RM, Zibbell JE, Green TC, Younkin M, et al. Characteristics of Fentanyl Overdose - Massachusetts, 2014-2016. MMWR Morb Mortal Wkly Rep 2017 Apr 14;66(14):382-386 [FREE Full text] [CrossRef] [Medline]
    5. Scholl L, Seth P, Kariisa M, Wilson N, Baldwin G. Drug and Opioid-Involved Overdose Deaths - United States, 2013-2017. MMWR Morb Mortal Wkly Rep 2018 Jan 04;67(5152):1419-1427 [FREE Full text] [CrossRef] [Medline]
    6. Ahmad FB, Cisewski JA, Rossen LM, Sutton P. Provisional drug overdose death counts. National Center for Health Statistics. 2023.   URL: https://www.cdc.gov/nchs/nvss/vsrr/drug-overdose-data.htm [accessed 2023-01-19]
    7. Watts M. Overdose deaths accelerating COVID-19 pandemic, expanded prevention efforts needed. Florida Department of Health in Paso County. 2021.   URL: http://pasco.floridahealth.gov/newsroom/2021/02/overdose-deaths-accelerating-​during-covid-19.html [accessed 2023-01-19]
    8. AMA issue brief: drug overdose epidemic worsened during COVID pandemic. Alliant Health Solutions. 2021.   URL: https:/​/www.​alliantquality.org/​media_library/​ama-issue-brief-drug-overdose-epidemic-worsened-during-covid-pandemic/​ [accessed 2023-01-19]
    9. Data Brief: Opioid-Related Overdose Deaths among Massachusetts Residents 2022 December. Current opioid statistics | Mass.gov. 2022 Dec.   URL: https://www.mass.gov/lists/current-opioid-statistics#updated-data-%E2%80%93-as-of-​december-2022- [accessed 2023-01-18]
    10. Data brief: opioid-related overdose deaths among Massachusetts residents (2020-11). Massachusetts Department of Public Health. 2020.   URL: https://archives.lib.state.ma.us/handle/2452/837612 [accessed 2023-01-19]
    11. Number of unintentional opioid-related deaths by county, MA residents: 2000-2015. Massachusetts Department of Public Health. 2016.   URL: http:/​/www.​mass.gov/​eohhs/​docs/​dph/​quality/​drugcontrol/​county-level-pmp/​overdose-deaths-by-​county-including-map-may-2016.​pdf [accessed 2023-01-19]
    12. Bharel M, Bernson D, Averbach A. Using Data to Guide Action in Response to the Public Health Crisis of Opioid Overdoses. NEJM Catalyst 2020 Sep;1(5). [CrossRef]
    13. Marks C, Carrasco-Escobar G, Carrasco-Hernández R, Johnson D, Ciccarone D, Strathdee SA, et al. Methodological approaches for the prediction of opioid use-related epidemics in the United States: a narrative review and cross-disciplinary call to action. Transl Res 2021 Aug;234:88-113 [FREE Full text] [CrossRef] [Medline]
    14. Bozorgi P, Porter DE, Eberth JM, Eidson JP, Karami A. The leading neighborhood-level predictors of drug overdose: A mixed machine learning and spatial approach. Drug Alcohol Depend 2021 Dec 01;229(Pt B):109143. [CrossRef] [Medline]
    15. Schell RC, Allen B, Goedel WC, Hallowell BD, Scagos R, Li Y, et al. Identifying predictors of opioid overdose death at a neighborhood level with machine learning. Am J Epidemiol 2022 Feb 19;191(3):526-533 [FREE Full text] [CrossRef] [Medline]
    16. Abell-Hart K, Rashidian S, Teng D, Rosenthal RN, Wang F. Where Opioid Overdose Patients Live Far From Treatment: Geospatial Analysis of Underserved Populations in New York State. JMIR Public Health Surveill 2022 Apr 12;8(4):e32133 [FREE Full text] [CrossRef] [Medline]
    17. Basak A, Cadena J, Marathe A, Vullikanti A. Detection of Spatiotemporal Prescription Opioid Hot Spots With Network Scan Statistics: Multistate Analysis. JMIR Public Health Surveill 2019 Jun 17;5(2):e12110 [FREE Full text] [CrossRef] [Medline]
    18. Lyle Cooper R, Thompson J, Edgerton R, Watson J, MacMaster SA, Kalliny M, et al. Modeling dynamics of fatal opioid overdose by state and across time. Prev Med Rep 2020 Dec;20:101184 [FREE Full text] [CrossRef] [Medline]
    19. Campo DS, Gussler JW, Sue A, Skums P, Khudyakov Y. Accurate spatiotemporal mapping of drug overdose deaths by machine learning of drug-related web-searches. PLoS One 2020 Dec 7;15(12):e0243622 [FREE Full text] [CrossRef] [Medline]
    20. Sumetsky N, Mair C, Wheeler-Martin K, Cerda M, Waller LA, Ponicki WR, et al. Predicting the Future Course of Opioid Overdose Mortality: An Example From Two US States. Epidemiology 2021 Jan;32(1):61-69 [FREE Full text] [CrossRef] [Medline]
    21. Bauer C, Champagne-Langabeer T, Bakos-Block C, Zhang K, Persse D, Langabeer JR. Patterns and risk factors of opioid-suspected EMS overdose in Houston metropolitan area, 2015-2019: A Bayesian spatiotemporal analysis. PLoS One 2021 Mar 11;16(3):e0247050 [FREE Full text] [CrossRef] [Medline]
    22. Li ZR, Xie E, Crawford FW, Warren JL, McConnell K, Copple JT, et al. Suspected heroin-related overdoses incidents in Cincinnati, Ohio: A spatiotemporal analysis. PLoS Med 2019 Nov 12;16(11):e1002956 [FREE Full text] [CrossRef] [Medline]
    23. Jalal H, Buchanich JM, Roberts MS, Balmert LC, Zhang K, Burke DS. Changing dynamics of the drug overdose epidemic in the United States from 1979 through 2016. Science 2018 Sep 21;361(6408) [FREE Full text] [CrossRef] [Medline]
    24. ZCTA level fatal opioid overdose rate in Massachusetts (MA).   URL: https://spatiotemporal-data-science.shinyapps.io/MA_Interactive_Map/ [accessed 2023-01-19]
    25. American Community Survey (ACS). United States Census Bureau.   URL: https://www.census.gov/programs-surveys/acs [accessed 2023-01-19]
    26. Knorr-Held L, Besag J. Modelling risk from a disease in time and space. Stat Med 1998 Sep 30;17(18):2045-2060. [CrossRef] [Medline]
    27. Bauer C, Wakefield J, Rue H, Self S, Feng Z, Wang Y. Bayesian penalized spline models for the analysis of spatio-temporal count data. Stat Med 2016 May 20;35(11):1848-1865 [FREE Full text] [CrossRef] [Medline]
    28. Rue H, Held L. Gaussian Markov Random Fields: Theory and Applications. London: Chapman & Hall; 2005.
    29. Besag J. Spatial Interaction and the Statistical Analysis of Lattice Systems. Journal of the Royal Statistical Society: Series B (Methodological) 2018 Dec 05;36(2):192-225. [CrossRef]
    30. RStudio is an integrated development environment (IDE) for R. GitHub.   URL: https://github.com/rstudio/rstudio [accessed 2023-01-19]
    31. Lindgren F, Rue H. Bayesian spatial modelling with R-INLA. J Stat Softw 2015;63(19):1-25. [CrossRef]
    32. Sprague Martinez L, Rapkin BD, Young A, Freisthler B, Glasgow L, Hunt T, et al. Community engagement to implement evidence-based practices in the HEALing communities study. Drug Alcohol Depend 2020 Dec 01;217:108326 [FREE Full text] [CrossRef] [Medline]
    33. Data Brief: Opioid-Related Overdose Deaths among Massachusetts Residents 2022 December. Massachusetts Department of Public Health. 2022 Dec.   URL: https://www.mass.gov/lists/current-opioid-statistics#updated-data-%E2%80%93-as-of-​december-2022- [accessed 2023-01-19]
    34. Slavova S, Rock P, Bush HM, Quesinberry D, Walsh SL. Signal of increased opioid overdose during COVID-19 from emergency medical services data. Drug Alcohol Depend 2020 Sep 01;214:108176 [FREE Full text] [CrossRef] [Medline]
    35. Borquez A, Martin NK. Fatal overdose: Predicting to prevent. Int J Drug Policy 2022 Jun;104:103677 [FREE Full text] [CrossRef] [Medline]
    36. Data Brief: Chapter 55 – 2017 Legislative Report. Public Health Data Warehouse (PHD): Publications. 2017.   URL: https://www.mass.gov/lists/public-health-data-warehouse-phd-publications [accessed 2022-01-19]
    37. U.S. Congress. Senate. Subcommittee on Africa and Global Health Policy of the Committee on Foreign Relations. Zimbabwe After the Elections: Hearing before the Subcommittee on Africa and Global Health Policy of the Committee on Foreign Relations. In: The SHAFR Guide Online. Leiden: BRILL; 2022.
    38. COVID-19 Forecasts: Deaths. Centers for Disease Control and Prevention.   URL: https://www.cdc.gov/coronavirus/2019-ncov/science/forecasting/forecasting-us.html [accessed 2023-01-18]

    ‎
    CAR: conditional autoregressive model
    CrI: credible interval
    MAE: mean absolute error
    MDPH: Massachusetts Department of Public Health
    RD: rank difference
    RMSE: root mean square error
    RVRS: Registry of Vital Records and Statistics
    ZCTA: ZIP Code Tabulation Area

    Edited by A Mavragani, H Bradley; submitted 26.07.22; peer-reviewed by Z Li, WHJ Lo-Ciganic, S Linder; comments to author 24.10.22; revised version received 14.12.22; accepted 26.12.22; published 10.02.23

    Copyright

    ©Cici Bauer, Kehe Zhang, Wenjun Li, Dana Bernson, Olaf Dammann, Marc R LaRochelle, Thomas J Stopka. Originally published in JMIR Public Health and Surveillance (https://publichealth.jmir.org), 10.02.2023.

    This is an open-access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work, first published in JMIR Public Health and Surveillance, is properly cited. The complete bibliographic information, a link to the original publication on https://publichealth.jmir.org, as well as this copyright and license information must be included.