Skip to main content
Download PDF

Risk factors for persistent fatal opioid-involved overdose clusters in Massachusetts 2011–2021: a spatial statistical analysis with socio-economic, accessibility, and prescription factors

BMC Public Health volume 24, Article number: 1893 (2024) Cite this article

Abstract

Background

Fatal opioid-involved overdose rates increased precipitously from 5.0 per 100,000 population to 33.5 in Massachusetts between 1999 and 2022.

Methods

We used spatial rate smoothing techniques to identify persistent opioid overdose-involved fatality clusters at the ZIP Code Tabulation Area (ZCTA) level. Rate smoothing techniques were employed to identify locations of high fatal opioid overdose rates where population counts were low. In Massachusetts, this included areas with both sparse data and low population density. We used Local Indicators of Spatial Association (LISA) cluster analyses with the raw incidence rates, and the Empirical Bayes smoothed rates to identify clusters from 2011 to 2021. We also estimated Empirical Bayes LISA cluster estimates to identify clusters during the same period. We constructed measures of the socio-built environment and potentially inappropriate prescribing using principal components analysis. The resulting measures were used as covariates in Conditional Autoregressive Bayesian models that acknowledge spatial autocorrelation to predict both, if a ZCTA was part of an opioid-involved cluster for fatal overdose rates, as well as the number of times that it was part of a cluster of high incidence rates.

Results

LISA clusters for smoothed data were able to identify whether a ZCTA was part of a opioid involved fatality incidence cluster earlier in the study period, when compared to LISA clusters based on raw rates. PCA helped in identifying unique socio-environmental factors, such as minoritized populations and poverty, potentially inappropriate prescribing, access to amenities, and rurality by combining socioeconomic, built environment and prescription variables that were highly correlated with each other. In all models except for those that used raw rates to estimate whether a ZCTA was part of a high fatality cluster, opioid overdose fatality clusters in Massachusetts had high percentages of Black and Hispanic residents, and households experiencing poverty. The models that were fitted on Empirical Bayes LISA identified this phenomenon earlier in the study period than the raw rate LISA. However, all the models identified minoritized populations and poverty as significant factors in predicting the persistence of a ZCTA being part of a high opioid overdose cluster during this time period.

Conclusion

Conducting spatially robust analyses may help inform policies to identify community-level risks for opioid-involved overdose deaths sooner than depending on raw incidence rates alone. The results can help inform policy makers and planners about locations of persistent risk.

Peer Review reports

Introduction

The opioid overdose crisis has persisted as one of the most significant public health challenges of the past two decades in the US. Between 1999 and 2019, the national age-adjusted opioid-involved overdose rate increased from 2.9 per 100,000 population to 15.5 per 100,000 [1]. Rates continued to increase steeply during the COVID-19 pandemic [2]. The rate climbed to 21.4 per 100,000 during 2020 and then to 24.7 per 100,000 in 2021 [3]. Commonly understood drivers of the overdose crisis include drug supply and demand-side pressures. Supply-side drivers have evolved over multiple waves, from prescription opioids, to heroin, to illicitly manufactured fentanyl [4, 5]. Demand-side drivers include deindustrialization and concentrated poverty, pain arising from work-related injuries, [5] income inequality, [6] and added stress, isolation, and economic disadvantage connected to the COVID-19 pandemic [7, 8].

More studies have begun to consider the association between drug use, opioid-related mortality, and the built environment, defined by Ezell and colleagues as “the purposeful creation and spatial arrangement of housing, sidewalks, roadways, retail and institutional buildings, public transit, and green spaces.” [9,10,11] Research has already suggested that the built environment influences health and health behaviors [12], including substance use [9, 10, 13] and opioid-related mortality [11]. For example, analgesic opioid-involved overdose fatalities were found to be more likely to occur in “fragmented” neighborhoods than in higher-income neighborhoods in New York City [14]. Chichester et al. found that bus stops and public schools were associated with increased risk of opioid overdose in rural areas of an Alabama county and that inpatient treatment centers, transitional living facilities, express loan establishments, and liquor vendors were associated with increased opioid overdose risk in urban areas of the same county [15]. Inequality and racial and ethnic composition of neighborhoods have also been identified as correlates of increased opioid mortality [16, 17].

The connection between substance use and built environment variables (access to public restrooms, access to pharmacies, and driving distance to services, defined in our study as fast-food restaurants, gas stations, and highway exits) is also important. Prior studies found that people who inject drugs often use drugs in public restrooms [9, 18, 19]. Pharmacies represent an important access variable for several reasons. During the initial wave of the overdose crisis, pharmaceutical prescriptions, either legitimate, diverted, or potentially inappropriate, fed the opioid supply [4, 17, 20]. In addition, naloxone (a medication to reverse an opioid overdose) is available at pharmacies without a prescription, [21, 22] although this provision may vary by neighborhood socio-demographic levels [23]. For example, prescription opioid poisoning increased more in postal codes with greater pharmacy density in California in the time period between 2001 and 2010 [24]. Road access to services may mediate several factors related to substance use, such as access to meeting places to buy illicit substances, as well as access to harm reduction services including syringe services programs [9].

Over the past two decades, Massachusetts’ opioid-involved mortality rate has often been higher than the national rate and, at times, twice as high [25]. The most recent data released by the Massachusetts Department of Public Health indicated that opioid-involved mortality rate peaked at 33.5 per 100,000 population in 2022 [26]. Studies focusing on the overdose crisis in Massachusetts have explored potential intersections between the built environment and opioid overdoses. A spatial analysis of potentially inappropriate opioid prescribing (PIP) identified several overdose and PIP clusters, but did not find a significant overlap between the two [20]. Other studies identified a rural county in Massachusetts with both good access to harm reduction services but persistently high overdose rates, [27] and that a majority of overdoses in the state occurred at home between 2015 and 2017 [28]. A recent study of opioid-related deaths in Massachusetts incorporated psychosocial, economic, built environment, and health-related variables using multilevel mixed-effects regression models, found that none of the built environment variables had a statistically significant association with opioid-related mortality [29].

While these studies across the US and within Massachusetts have contributed importantly to our understanding of the relationships between substance use and the built environment, many often fail to use spatial statistical methods that properly account for excess zero counts for overdose outcomes and spatial autocorrelation, two of the most common confounding factors in spatial epidemiology [30]. With spatial analysis of rare events, such as overdose deaths, the distribution of death counts is often “zero-inflated” (i.e., a large number of locations have no deaths while a few locations have many). Several methods have been used to address this issue, including use of: small area estimation techniques; [31] a rate smoothing technique to examine characteristics of prescription opioid poisoning; smoothing rates in the ZIP Code Tabulation Area (ZCTA) with small populations; and Empirical Bayes smoothing. [32,33,34,35] These methods include data from adjacent spatial units that enable refined estimates for all locations, but are of elevated importance in locations that have low population density such as rural areas in western Massachusetts. Two models that acknowledge spatial autocorrelation and account for excess zeros include zero-inflated Poisson regression models, [36, 37] such as the Besag-York-Mollie model, [38] and a Bayesian hierarchical space–time misalignment Poisson model [24, 39].

The problem of excess zeros masking potentially non-zero true counts and rates exists in Massachusetts analyses because the geographical distribution of opioid overdose deaths in Massachusetts spans urban, suburban, and sparsely populated rural areas that are likely to have very low or zero counts. To date, we are only aware of a few studies of opioid-involved overdose deaths using Massachusetts data that have employed techniques that acknowledge zero-inflation and spatial autocorrelation at smaller spatial units (such as the ZCTA level), [40] but these techniques could greatly help in enhancing our understanding of geographic variation in overdose rates across Massachusetts.

The goal of our study was to identify fatal overdose clusters as well as socioeconomic and built environment factors associated with opioid-related overdose rates in Massachusetts from 2011 to 2021 while acknowledging spatial autocorrelation. We aimed to: (1) address the issue of zero-inflated opioid-involved overdose rates by comparing three spatial methods (raw rates, Empirical Bayes, and Empirical Bayes Spatial); (2) utilize the raw and Empirical Bayes rate as well as Empirical Bayes Local Indicator of Spatial Autocorrelation (LISA) to identify statistically significant fatal opioid overdose clusters in Massachusetts between 2011 and 2021; (3) derive socio-environmental and pharmacological variables using PCA in those clusters; and (4) model community-level factors associated with overdose rates using Conditional Autoregressive Empirical Bayes logistic regression models to predict if a ZCTA was part of a cluster and Conditional Autoregressive Empirical Bayes zero inflated Poisson regression models to predict how many times it was part of a cluster during the time period while acknowledging spatial autocorrelation of the outcome (Fig. 1).

Data

Data sources

We obtained fatal opioid-related overdose data by address of residence for decedents between 2011 and 2021 from the Massachusetts Registry of Vital Records and Statistics (RVRS) [41]. These data included sociodemographic characteristics of decedents, including sex, race, ethnicity, and age at the time of death. We obtained opioid prescription data from the Massachusetts Prescription Monitoring Program (MA PMP), aggregated at the zip code or ZCTA level [42]. We obtained address level data for services (gas stations, fast food restaurants, pharmacies), and “access to infrastructure” measures (highway exits, major roads) from Data Axle and from MassGIS (Massachusetts Bureau of Geographic Information) [43, 44]. We obtained pharmacy addresses from the Massachusetts Board of Registration in Pharmacy [22]. Finally, we compiled population level sociodemographic data at the ZCTA level, from the US Census Bureau’s American Community Survey (ACS), relying on 2011 to 2015 and the 2017 to 2021 5-year estimates for people aged ≥ 10 years for calculating rates [45]. We used ACS 2011–2015 data for the years 2011–2016 and 2017–2021 data for the years 2017–2021.

Outcome

The two outcomes were whether a ZCTA was part of a cluster or not and the number of times a ZCTA was within a fatal overdose cluster based on annual opioid-related overdose rates between 2011 and 2021. We calculated the LISA statistic for each year between 2011 and 2021 for the raw rates, and for Empirical Bayes smoothed rates. We also calculated the Empirical Bayes LISA for the same time period, as this adjusts for variation in the underlying population. [46, 47] We did not use the spatial EB rates for calculating clusters. This is because spatial EB rates were useful to identify locations at risk, but they also resulted in over-smoothing of the rates.

Covariates

PIP measures were aggregated by ZCTA. We defined PIP measures, established through our previous research, for the years 2011–2017 at the ZCTA scale: high Morphine Milligram Equivalents or MME ( > = 90), co-prescribing of benzodiazepines and opioids, poly-pharmacy opioid prescription fills, multiple provider episodes (i.e., doctor shopping), >=3 cash purchases of opioid prescriptions and opioid prescriptions without a pain diagnosis. We were not able to obtain PIP data for the years 2018–2021.

Socioeconomic measures

We selected socioeconomic measures from the ACS at the ZCTA-level based on the literature and our previous research [20, 48, 49]. Poverty as a percentage of households living below the poverty threshold was defined by the five-year ACS estimates. We also included median age, and population percentages by race and ethnicity for non-Hispanic White, Black, Asian, American Indian and Alaskan Native, and Hispanic communities.

Built environment measures

We compiled built environment measures based on the literature and our previous research [20, 48, 49]. Specifically, we selected gas stations and fast-food restaurant locations as a proxy for access to public restrooms, locations where overdoses often occur [18, 19, 43, 50, 51]. We used pharmacy addresses to analyze the spatial distribution of access to sources of over-the-counter naloxone [22]. We compiled highway exits and major roads from a Massachusetts GIS agency (MassGIS) as a proxy for access to services. For each ZCTA we calculated the distance to the nearest exit, major road, pharmacy, fast food restaurant, and gas station from the centroid of the ZCTA. We also calculated the average gas station density for each ZCTA.

Methods

Figure 1 provides an overview of the data and methods used in this paper.

Death rate mapping

Using ArcGIS Pro 3.1 (ESRI, Redlands, CA) we created maps of raw and Empirical Bayes and Empirical Bayes Spatial smoothed opioid overdose death rates at the ZCTA level across the state.

Spatial smoothing techniques

Smoothing techniques work by detrending the overdose rate within target polygons (ZCTAs) and by using data from neighboring polygons, thus allowing for the calculation of a local average overdose measure that is less susceptible to variation due to outlier values. We employed two spatial smoothing techniques to derive stable estimates for spatial measures (fatal overdose rates). These included: (1) Empirical Bayes Method; and, (2) Empirical Bayes Spatial method [52,53,54].

The Empirical Bayes smoothing for rate calculations relies on calculating a weighted average of the raw rate for each ZCTA and the state average, with weights proportional to the underlying population at risk. Therefore, ZCTAs with a small population at risk will tend to have their rates adjusted considerably, while rates for ZCTAs with larger populations at risk will not change much. The second method was Empirical Bayes Spatial smoothing, which is like Empirical Bayes smoothing, except that the reference rate is computed for a group of neighboring ZCTAs that share a boundary with each individual ZCTA, rather than taking the same overall reference rate for all ZCTAs. We used the R package rgeoda to calculate Empirical Bayes and Empirical Bayes Spatial smoothed rates [53, 55]. We then used the LISA statistic to identify and map clusters (i.e., ZCTAs with high overdose rates surrounded by ZCTAs with high rates or high-high (HH) locations and ZCTAs with high rates surrounded by low rates or high-low (HL) locations for (i) the raw rates, (ii) the Empirical Bayes smoothed rate and (iii) using the Empirical Bayes LISA method [47, 53]. We did not calculate LISA for Empirical Bayes Spatial rates to avoid over-smoothing the data. Clusters were considered statistically significant if the local pseudo-p value p < 0.05 as estimated by rgeoda. LISA were also calculated using the rgeoda package. We used the queen’s contiguity criterion to define all the spatial weights matrix calculations.

Principal components analysis (PCA)

To avoid multicollinearity amongst the PIP, access, and socioeconomic variables due to significant cross correlation, we used PCA with a Varimax rotation to extract new variables that summarized the covariance in the PIP, socioeconomic, and access variables. PCA is a dimension reduction technique that is commonly used to reduce the number of variables used in analyses while preserving the information from them. We used a cut point of 0.25 for the included factor loadings as recommended in other studies [56]. The psych and dplyr packages in R were used for PCA and data management [57].

Statistical modeling

We fit three Bayesian logistic regression models where the dichotomous outcome characterized whether a ZCTA was in a HH or HL cluster or not (0/1), as identified by (i) the raw overdose rate LISA, (ii) the Empirical Bayes rate LISA and (iii) the Empirical Bayes-corrected LISA developed by Assuncao and Reis [47]. The components derived from the PCA were the explanatory variables for all three models. Models fit with all the variables (without PCA) resulted in very high VIF values due to multicollinearity. We used Conditional Autoregressive (CAR) Bayes regression models to fit both the logistic regression model focused on whether the ZCTA was in a cluster for individual years between 2011 and 2021, as well as a zero-inflated Poisson regression model for the number of times the ZCTA was located within a LISA cluster, as identified by the smoothing methods in the 11 years between 2011 and 2021. Zero-inflated Poisson regression models are fitted when the outcome is a count with excess zeros. CAR models, which acknowledged the spatial dependency of the outcome within a Bayesian framework, were used for both the logistic as well as the ZIP regression models. The package CARBayes was used to fit the models in R.

Results

The use of rate smoothing techniques facilitated the identification of high opioid overdose rate ZCTAs surrounding Pittsfield, Springfield, Fall River, and New Bedford, Massachusetts. Several of these locations had high overdose fatality rates in 2021 based on the smoothed data that were not easily identifiable in maps of the raw rates for the same year (Fig. 2). Figure 3 shows that the impact varied over time. For example, in 2011 the smoothed rates enabled identification of clusters in the Boston area, in addition to the area south of Worcester, Fall River, New Bedford, and Cape Cod (Fig. 3). Likewise, in 2014, using the raw rate alone resulted in fewer and more scattered clusters − 25 than the Empirical Bayes methods, which identified 45 and 53 ZCTAs as being part of LISA HH or HL clusters. The Empirical Bayes LISA and LISA of the Empirical Bayes rate methods also depicted more geographic variability, facilitating identification of clusters near Wareham in the southeast and Haverhill in the northeast. In 2017, the numbers of clusters identified based on the raw rates alone was comparable to those from the other methods − 36 versus 28 from the Empirical Bayes LISA method. The mapped results suggested that overdose deaths persisted in Worcester, Fall River, and Boston. By 2021, the number of ZCTAs identified within high overdose clusters using raw rates was much lower: 13. However, using Empirical Bayes LISA, the number of ZCTAs in clusters was 38. Figure 4 highlights locations of persistent HH and HL overdose clusters, with the most notable ongoing clusters in Worcester, New Bedford, Fall River, and Wareham. Many of these were not identified by the raw rate LISA, which only showed clusters in western Massachusetts in 2021 (Fig. 3).

We compared sociodemographic and built environment access variables in the statistically significant ZCTA LISA clusters, which illuminated important differences in sociodemographic factors. In ZCTAs that were always identified as being part of HH or HL LISA clusters, the mean percentages of the population living in poverty, and residents who were non-Hispanic Black or Hispanic, were higher than in ZCTAs that never appeared in fatal overdose clusters (Table 1). In terms of built environment variables, ZCTAs that were located within persistent clusters were closer to major roads and highway exits.

The results of the PCA suggested that the four components we used in this analysis explained about 62% of the variance in the data (Table 2). The loadings on the first PCA component suggested that it was a measure of PIP because it had a positive loading (indicating positive correlation) with all PIP measures. This component had high scores in the central parts of the state along Interstate 90 and along Interstate 95 in the south (Fig. 5). The second component was correlated with poverty, non-Hispanic Black, American Indian and Alaska Native (AIAN) and Hispanic population percentages and we refer to it as the “minority and poverty” component. Locations that had high values for the minority and poverty component were found across Massachusetts in locations like Worcester, Springfield, New Bedford, Fall River, and in the southern area of Boston, within the Dorchester neighborhood (Fig. 5).

The third component, which we termed “distance to amenities,” was related to being relatively far from fast food restaurants, pharmacies, and gas stations. The ZCTAs with high values for the distance from amenities component were located not only in high poverty locations, such as urban neighborhoods in Boston, Worcester, New Bedford, and Springfield, but also in the more rural areas of Western Massachusetts, and on Cape Cod and the islands. Some high-income suburbs close to Boston like Carlisle, and Grafton (which is near Worcester) also had moderately high values for this component suggesting that it was not correlated with income (Fig. 5).

The fourth component, which we referred to as the “rurality” component, loaded positively on median age, and distance to a highway and major roads. As can be seen in Fig. 5, high scores for the rurality component were identified in western Massachusetts and Cape Cod, while low scores were observed in urban areas within and near Boston, Springfield, New Bedford, and Worcester.

The results of the CAR Bayesian logistic regressions for the outcome focused on whether a ZCTA was in a HH or HL cluster (derived from three different LISA statistics) for each year from 2011 to 2021 are shown in Table 3. We only included coefficients with statistically significant credible intervals from the CAR Bayesian regressions in Table 3 for brevity. Notably, the Empirical Bayes LISA, which was the most stable indicator of high overdose clusters observed in the maps, was the earliest to pick up any trends in the models over time. The LISA on Empirical Bayes smoothed rates also had similar results. LISA estimates of high overdose clusters from raw rates eventually also showed similarly significant credible intervals, but only after the trends had been observed in the other two model estimates.

The results of the CAR Bayes logistic model, predicting whether ZCTAs were located within a high overdose rate cluster using the Empirical Bayes LISA, suggested that for most years between 2011 and 2021, fatal opioid-involved overdose cluster ZCTAs were more likely to have a high minority and poverty component. The credible intervals for the minority and poverty component were lowest early in the study period in 2012 [2.13–4.47] and rose in 2021 [4.1–51.9] but they were always positive and significant. The PIP component was not significant in any of the models that were fitted. The distance to amenities and rural components, when significant, had a protective effect in that ZCTAs with high scores on these variables meaning that locations that were far from amenities were less likely to be in high overdose clusters. The odds ratio on the distance to amenities component varied from a higher range [0.29–0.89] in 2018 to a smaller range in in 2021 of [0.001–0.06]. The 95% credible intervals for the rurality component odds ratio varied from a lower range more recently in 2019 [0.05–0.67] to a wider range [0.07–0.99] in 2014.

The CAR Bayesian models predicting the number of times a ZCTA was in a high overdose rate cluster (Table 4) always had a positive and significant value within the credible interval for the minority and poverty component for the three models. This suggests that ZCTAs with high scores on this component in Worcester, Springfield, New Bedford, Brockton, Boston, and Fall River were more likely to have been a persistent overdose cluster over the 11-year period. The negative coefficients for the distance to amenities and rurality suggested that ZCTAs scoring high on distance to amenities and rurality were less likely to be located within persistent high overdose rate clusters (Table 4). The PIP component and the rurality component were not significant in any of the models.

Discussion

In this study, we identified overdose clusters using rate smoothing techniques and found persistent and significant high incidence clusters in Massachusetts. We estimated PCA components to account for multicollinearity while incorporating opioid prescription, socio-environmental and built environment variables; and we acknowledged spatial autocorrelation in our models in the presence of spatial dependency in both the location of high overdose clusters and their persistence over the study period. Using PCA, we identified four unique components in Massachusetts that described ZCTA characteristics: PIP, Minority and Poverty, Distance to amenities, and Rurality. In the logistic regression models predicting if a ZCTA was part of a cluster at high risk of opioid-related overdose, we found that the minority and poverty factor was a significant predictor of a high overdose cluster in most years. In predicting the count of the number of times a ZCTA was part of a cluster, we found that the minority and poverty factor was again a significant predictor. Some recent studies have also noted the value of constructing variables that measure social disadvantage through the use of PCA [58, 59].

Zero inflation of disease rates is common in spatial epidemiological studies, which has resulted in several techniques to address this issue [31,32,33,34,35,36,37,38,39, 60]. We chose Empirical Bayes, and Empirical Bayes spatial smoothing methods, which reduced variance in the opioid fatality rates and thus resulting in better identification of high incidence clusters in ZCTAs that had relatively low population density. These rate smoothing techniques offer considerable promise in identifying regional clusters of persistent risk for fatal overdoses in Massachusetts earlier in the overdose epidemic than raw rates. Our application of the Empirical Bayes methods supported the use of this technique as a method to detect patterns at smaller spatial scales in which there was a mix of rural and urban areas [34]. Another recent study in Illinois noted the usefulness of EB clusters in a temporal study that could “serve as a proxy for pervasive risk” and “emerging risk” [61].

Our results may also indicate that minority and poverty presence in a ZCTA is significant in predicting high incidence clusters in Massachusetts. Several recent studies have also noted the role of poverty and race in predicting opioid overdose incidence in other states [16, 62,63,64,65,66]. The protective effect suggested by distance to amenities and rurality in the models suggests that measures of built environment and access to opioid supply may need to be incorporated into future analyses [16, 65]. It should also be noted that the negative association with rurality was no longer significant after 2019 which may be worth investigating in future work. A recent Massachusetts Department of Public Health report also noted this resurgence of opioid overdose prevalence in the most rural areas of the state in recent years [26].

The mapped LISA cluster estimates demonstrate that high overdose rate clusters that persisted over the years were comprised of higher percentages of non-Hispanic Black and AIAN populations, as well as Hispanic populations, and higher percentages of people living in poverty [66]. In 2011–2015, the overdose epidemic in Massachusetts and elsewhere was characterized in the media as affecting people who were middle income and non-Hispanic White as part of the “deaths of despair” narrative [67,68,69,70,71]. The role of poverty in opioid-involved overdoses in minoritized people who use drugs should be surveilled carefully as opioid overdoses continue to rise in these groups [72, 73, 74, 75].

Our findings have several limitations. For privacy reasons, the prescribing data (PIP) in this study were only available at the ZCTA scale and, therefore, our results are subject to ecological fallacy. Models that are at the individual level may show different relationships than our findings, which were aggregated to the ZCTA. Furthermore, we could not obtain PIP data for 2018–2021 and the PCA is based only on 2011–2017 data. However, the data for the earlier time periods were spatially clustered so it is likely that the later time periods will also show similar spatial and statistical correlations. The statistical associations that we identified may only be applicable to Massachusetts. Several laws and regulations were enacted in the state during the opioid overdose epidemic that may not be generalizable to other states. In 2019, for instance, the Governor of Massachusetts enacted a law that allowed no more than a seven-day supply of prescribed opioids (with certain exceptions), which, along with other policies, resulted in a notable decline in the number of opioid prescriptions dispensed [76].

Conclusion

We employed a range of spatial analytical and statistical modeling approaches to better understand the opioid overdose crisis in Massachusetts and improve targeting of public health interventions. Smoothing methods allowed us to derive more stable estimates for opioid overdose rates and enhanced our understanding of the spatial clustering patterns related to fatal opioid-involved overdoses. We recommend using Empirical Bayes LISA when estimating clusters and Empirical Bayes smoothing for deriving rates especially for identifying locations at risk, particularly when tracking places that have a variety of population densities over space and time. These methods were able to identify clusters earlier in the epidemic than raw rates were able to in low population density locations. PCA facilitated a better understanding of unique community-level spatial and socioeconomic variables that may explain opioid-involved incidence. By using components that combined several socioeconomic, built environment, and prescription variables in the regression models, which also acknowledged spatial autocorrelation, we were able to characterize the significant contributing factors to opioid-related deaths in Massachusetts. Future research should investigate minority and poverty variables and their potential proxies. Many factors impact the cascade of opioid use, opioid use disorder, and opioid-involved overdose mortality, and by identifying factors at the beginning of the cascade, as we have done here, we can inform policies to intervene early in the cascade and prevent opioid-involved deaths.

Table 1 Mean (SD) of demographic and built environment variables by number of years in which a ZIP code tabulation area (ZCTA) was an empirical Bayes (EB) LISA HH (high surrounded by high) or HL (high surrounded by low)
Full size table
Table 2 Principal Components Analysis using a Varimax rotation: loadings of the opioid prescription, socioeconomic, and infrastructure variables for the first four components
Full size table
Table 3 Significant coefficients within the 95% credible interval (for CAR bayesian logistic model) associated with a ZCTA was a LISA (Local Indicators of Spatial Association) HH (high surrounded by high) or HL (high surrounded by low), Massachusetts, 2011–2021
Full size table
Table 4 CAR bayesian zero inflated Poisson regression coefficients predicting the number of times the ZCTA was a LISA (Local Indicators of Spatial Association) HH (high surrounded by high) or HL (high surrounded by low), Massachusetts, 2011–2021
Full size table
Fig. 1
figure 1

Flow diagram visualizing data and methods described in this paper. Note: PIP is potentially inappropriate opioid prescribing. ACS represents data from the US Census Bureau’s American Community Survey. Data Axle was used for built environmental variables such as gas stations and pharmacies; MassGIS is the Massachusetts State GIS Data provider; LISA (Local Indicators of Spatial Association) statistics generate clusters

Full size image
Fig. 2
figure 2

Fatal opioid overdose rate quintiles in 2021 in Massachusetts by ZIP Code Tabulation Area (ZCTA): 2(a) Raw rate; 2(b) Empirical Bayes (EB) smoothed rates; 2(c) Spatial EB smoothed rates 2(d) Reference map. This comparison of maps displaying raw and smoothed overdose rates highlights that identification of local high incidence clusters is easier in the EB and EB spatial maps than in the raw rate map’s patchwork of high and low values

Full size image
Fig. 3
figure 3

Fatal opioid overdose LISA clusters and outliers (p < 0.05) for 2011, 2014, 2017 and 2021: 3(a) Raw, 3(b) Empirical Bayes (EB) LISA, 3(c) EB smoothed death rates LISA by ZIP Code Tabulation Area (ZCTA), Massachusetts. LISA (Local Indicators of Spatial Association) analyses generate clusters of HH (high ZCTAs surrounded by ZCTAs with high rates), HL (high surrounded by low rates), LL (low surrounded by low rates), LH (low surrounded by high rates) or not significant

Full size image
Fig. 4
figure 4

Number of times (i.e., years) the ZCTA was a HH or HL cluster between 2011–2021 using Empirical Bayes LISA (Local Indicators of Spatial Association)

Full size image
Fig. 5
figure 5

Principal Components Analysis scores mapped in deciles for Massachusetts by ZIP Code Tabulation Area (ZCTA): (a) RC1: potentially inappropriate prescribing (PIP); (b) RC2: Minority and poverty; (c) RC3: Distance to amenities (d) RC4: Rurality

Full size image

Data availability

The datasets used and/or analyzed during the current study are available from the corresponding author or the Massachusetts Department of Public Health upon request.

References

  1. Hedegaard H. Drug Overdose Deaths in the United States, 1999–2019. 2020;(394):8.

  2. National Center for Health Statistics. Provisional Drug Overdose Death Counts. Centers for Disease Control and Prevention. 2023 [cited 2023 Dec 7]. https://www.cdc.gov/nchs/nvss/vsrr/drug-overdose-data.htm.

  3. Spencer MR, Minino AM, Warner. Margaret. Drug Overdose Deaths in the United States, 2001–2021. NCHS Data Brief, no 457 [Internet]. Hyattsville, MD: National Center for Health Statistics; 2022. https://doi.org/10.15620/cdc:122556.

  4. Ciccarone D. The triple wave epidemic: Supply and demand drivers of the US opioid overdose crisis. Int J Drug Policy [Internet]. 2019 Feb 2 [cited 2019 Oct 31]; http://www.sciencedirect.com/science/article/pii/S0955395919300180.

  5. Cerdá M, Krawczyk N, Hamilton L, Rudolph KE, Friedman SR, Keyes KM. A critical review of the social and behavioral contributions to the overdose epidemic. Annu Rev Public Health. 2021;42(1):95–114.

    Article  PubMed  PubMed Central  Google Scholar 

  6. Thombs RP, Thombs DL, Jorgenson AK, Harris Braswell T. What is driving the drug overdose epidemic in the United States? J Health Soc Behav. 2020;61(3):275–89.

    Article  PubMed  Google Scholar 

  7. 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;214:108176.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Walters SM, Seal DW, Stopka TJ, Murphy ME, Jenkins WD. COVID-19 and people who use drugs–A commentary. Health Behav Policy Rev. 2020;7(5):489–97.

    Article  PubMed  PubMed Central  Google Scholar 

  9. Ezell JM, Ompad DC, Walters S. How urban and rural built environments influence the health attitudes and behaviors of people who use drugs. Health Place. 2021;69:102578.

    Article  PubMed  PubMed Central  Google Scholar 

  10. Li Y, Miller HJ, Root ED, Hyder A, Liu D. Understanding the role of urban social and physical environment in opioid overdose events using found geospatial data. Health Place. 2022;75:102792.

    Article  PubMed  Google Scholar 

  11. Hembree C, Galea S, Ahern J, Tracy M, Markham Piper T, Miller J, et al. The urban built environment and overdose mortality in New York City neighborhoods. Health Place. 2005;11(2):147–56.

    Article  CAS  PubMed  Google Scholar 

  12. Northridge ME, Sclar ED, Biswas P. Sorting out the connections between the built environment and health: a conceptual framework for navigating pathways and planning healthy cities. J Urban Health. 2003;80(4):556–68.

    Article  PubMed  PubMed Central  Google Scholar 

  13. Deering KN, Rusch M, Amram O, Chettiar J, Nguyen P, Feng CX, et al. Piloting a ‘spatial isolation’ index: the built environment and sexual and drug use risks to sex workers. Int J Drug Policy. 2014;25(3):533–42.

    Article  PubMed  Google Scholar 

  14. Cerdá M, Ransome Y, Keyes KM, Koenen KC, Tardiff K, Vlahov D, et al. Revisiting the role of the Urban Environment in Substance Use: the case of analgesic overdose fatalities. Am J Public Health. 2013;103(12):2252–60.

    Article  PubMed  PubMed Central  Google Scholar 

  15. Chichester K, Drawve G, Giménez-Santana A, Sisson M, McCleskey B, Dye DW, et al. Pharmacies and features of the built environment associated with opioid overdose: a geospatial comparison of rural and urban regions in Alabama, USA. Int J Drug Policy. 2020;79:102736.

    Article  PubMed  Google Scholar 

  16. Johnson LT, Shreve T. The ecology of overdose mortality in Philadelphia. Health Place. 2020;66:102430.

    Article  PubMed  Google Scholar 

  17. Sadler RC, Furr-Holden D. The epidemiology of opioid overdose in Flint and Genesee County, Michigan: implications for public health practice and intervention. Drug Alcohol Depend. 2019;204:107560.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Fozouni L, Buchheit B, Walley AY, Testa M, Chatterjee A. Public restrooms and the opioid epidemic. Subst Abuse. 2020;41(4):432–6.

    Article  Google Scholar 

  19. Sutter A. Curtis, Matt, Frost, Taeko. Public drug use in eight U.S. cities: Health risks and other factors associated with place of drug use. Int J Drug Policy. 2019;64:62–9.

    Article  PubMed  Google Scholar 

  20. Stopka TJ, Amaravadi H, Kaplan AR, Hoh R, Bernson D, Chui KKH, et al. Opioid overdose deaths and potentially inappropriate opioid prescribing practices (PIP): a spatial epidemiological study. Int J Drug Policy. 2019;68:37–45.

    Article  PubMed  PubMed Central  Google Scholar 

  21. Smart R, Pardo B, Davis CS. Systematic review of the emerging literature on the effectiveness of naloxone access laws in the United States. Addiction. 2021;116(1):6–17.

    Article  PubMed  Google Scholar 

  22. Massachusetts Bureau of Substance Addiction Services. Getting Naloxone from a Pharmacy [Internet]. https://www.mass.gov/service-details/getting-naloxone-from-a-pharmacy.

  23. Egan KL, Foster SE, Knudsen AN, Lee JGL. Naloxone availability in Retail pharmacies and Neighborhood inequities in Access. Am J Prev Med. 2020;58(5):699–702.

    Article  PubMed  Google Scholar 

  24. Cerdá M, Gaidus A, Keyes KM, Ponicki W, Martins S, Galea S, et al. Prescription opioid poisoning across urban and rural areas: identifying vulnerable groups and geographic areas: geography of prescription opioid poisoning. Addiction. 2017;112(1):103–12.

    Article  PubMed  Google Scholar 

  25. Kaiser Family Foundation. Opioid Overdose Death Rates and All Drug Overdose Death Rates per 100,000 Population (Age-Adjusted) [Internet]. KFF. 2022 [cited 2022 Nov 22]. https://www.kff.org/other/state-indicator/opioid-overdose-death-rates/.

  26. Massachusetts Department of Public Health. Data Brief: Opioid-Related Overdose Deaths Among Massachusetts Residents, June 2024 [Internet]. 2024 Jun. https://www.mass.gov/doc/opioid-related-overdose-deaths-among-ma-residents-june-2024-0/download.

  27. Stopka TJ, Jacque E, Kelso P, Guhn-Knight H, Nolte K, Hoskinson R, et al. The opioid epidemic in rural northern New England: an approach to epidemiologic, policy, and legal surveillance. Prev Med. 2019;128:105740.

    Article  PubMed  PubMed Central  Google Scholar 

  28. Pustz J, Srinivasan S, Larochelle MR, Walley AY, Stopka TJ. Relationships between places of residence, injury, and death: spatial and statistical analysis of fatal opioid overdoses across Massachusetts. Spat Spatio-Temporal Epidemiol. 2022;100541.

  29. Flores MW, Cook BL, Mullin B, Halperin-Goldstein G, Nathan A, Tenso K, et al. Associations between neighborhood-level factors and opioid-related mortality: a multi-level analysis using death certificate data. Addiction. 2020;115(10):1878–89.

    Article  PubMed  PubMed Central  Google Scholar 

  30. Wu CC, Chu YH, Shete S, Chen CH. Spatially varying effects of measured confounding variables on disease risk. Int J Health Geogr. 2021;20(1):45.

    Article  PubMed  PubMed Central  Google Scholar 

  31. Rossen LM, Khan D, Warner M. Trends and Geographic patterns in drug-poisoning death rates in the U.S., 1999–2009. Am J Prev Med. 2013;45(6):e19–25.

    Article  PubMed  PubMed Central  Google Scholar 

  32. Rossen LM, Khan D, Warner M. Hot spots in mortality from drug poisoning in the United States, 2007–2009. Health Place. 2014;26:14–20.

    Article  PubMed  Google Scholar 

  33. Sharareh N, Hess R, White S, Dunn A, Singer PM, Cochran J. A vulnerability assessment for the HCV infections associated with injection drug use. Prev Med. 2020;134:106040.

    Article  PubMed  Google Scholar 

  34. Hester L, Shi X, Morden N. Characterizing the geographic variation and risk factors of fatal prescription opioid poisoning in New Hampshire, 2003–2007. Ann GIS. 2012;18(2):99–108.

    Article  Google Scholar 

  35. McLuckie C, Pho MT, Ellis K, Navon L, Walblay K, Jenkins WD, et al. Identifying areas with Disproportionate Local Health Department Services Relative to opioid overdose, HIV and Hepatitis C diagnosis rates: a study of Rural Illinois. Int J Environ Res Public Health. 2019;16(6):989.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Hernandez A, Branscum AJ, Li J, MacKinnon NJ, Hincapie AL, Cuadros DF. Epidemiological and geospatial profile of the prescription opioid crisis in Ohio, United States. Sci Rep. 2020;10(1):4341.

    Article  PubMed  PubMed Central  Google Scholar 

  37. Cramb SM, Baade PD, White NM, Ryan LM, Mengersen KL. Inferring lung cancer risk factor patterns through joint bayesian spatio-temporal analysis. Cancer Epidemiol. 2015;39(3):430–9.

    Article  PubMed  Google Scholar 

  38. Besag J, York J, Mollié A. Bayesian image restoration, with two applications in spatial statistics. Ann Inst Stat Math. 1991;43(1):1–20.

    Article  Google Scholar 

  39. Zhu L, Waller LA, Ma J. Spatial-temporal disease mapping of illicit drug abuse or dependence in the presence of misaligned ZIP codes. GeoJournal. 2013;78(3):463–74.

    Article  PubMed  Google Scholar 

  40. Bauer C, Zhang K, Li W, Bernson D, Dammann O, LaRochelle MR, et al. Small area forecasting of opioid-related mortality: bayesian Spatiotemporal Dynamic modeling Approach. JMIR Public Health Surveill. 2023;9:e41450.

    Article  PubMed  PubMed Central  Google Scholar 

  41. Massachusetts Registry of Vital Records and Statistics. Opioid Overdose Death Data [Internet]. 2019 [cited 2020 Aug 16]. https://www.mass.gov/orgs/registry-of-vital-records-and-statistics.

  42. Massachusetts Prescription Monitoring Program. Opioid prescribing data [Internet]. 2019. https://www.mass.gov/orgs/prescription-monitoring-program.

  43. Data Axle Reference Solutions. U.S. Businesses [Internet]. 2019. http://www.referenceusa.com/Home/Home.

  44. Massachusetts Bureau of Geographic Information. MassGIS. [cited 2020 Nov 30]. MassGIS. https://www.mass.gov/orgs/massgis-bureau-of-geographic-information.

  45. U.S. Census Bureau. American Community Survey (ACS) 5-year estimates, 2012–2016 [Internet]. 2017. Available from: www.census.gov.

  46. Ord JK, Getis A. Local spatial Autocorrelation statistics: Distributional issues and an application. Geogr Anal. 1995;27(4):286–306.

    Article  Google Scholar 

  47. Assunção RM, Reis IA. A new proposal to adjust Moran’s I for Population Density. Stat Med. 1999;18.

  48. Rowe C, Santos GM, Vittinghoff E, Wheeler E, Davidson P, Coffin PO. Neighborhood-Level and Spatial Characteristics Associated with Lay Naloxone reversal events and opioid overdose deaths. J Urban Health. 2016;93(1):117–30.

    Article  PubMed  PubMed Central  Google Scholar 

  49. Chichester K, Drawve G, Sisson M, McCleskey B, Dye DW, Cropsey K. Examining the neighborhood-level socioeconomic characteristics associated with fatal overdose by type of drug involved and overdose setting. Addict Behav. 2020;111:106555.

    Article  PubMed  Google Scholar 

  50. Buchheit BM, Crable EL, Lipson SK, Drainoni ML, Walley AY. Opening the door to somebody who has a chance. – the experiences and perceptions of public safety personnel towards a public restroom overdose prevention alarm system. Int J Drug Policy. 2021;88:103038.

    Article  PubMed  Google Scholar 

  51. Wolfson-Stofko B, Bennett AS, Elliott L, Curtis R. Drug use in business bathrooms: an exploratory study of manager encounters in New York City. Int J Drug Policy. 2017;39:69–77.

    Article  PubMed  Google Scholar 

  52. Shi X, Duell E, Demidenko E, Onega T, Wilson B, Hoftiezer D. A polygon-based locally-weighted-average method for smoothing Disease Rates of small units. Epidemiology. 2007;18(5):523–8.

    Article  PubMed  Google Scholar 

  53. Anselin L. Exploring spatial data with GeoDa TM: a workbook. Chicago, IL: GeoDa; 2010.

    Google Scholar 

  54. Cressie N. Statistics for Spatial Data. In: Statistics for Spatial Data [Internet]. John Wiley & Sons, Ltd; 1993 [cited 2021 Nov 23]. (Wiley Series in Probability and Statistics). https://onlinelibrary.wiley.com/doi/abs/10.1002/9781119115151.ch1.

  55. Anselin L, Lozano N, Koschinsky J. Rate transformations and smoothing. Urbana. 2006;51.

  56. Costello AB, Osborne J. Best practices in exploratory factor analysis: four recommendations for getting the most from your analysis. Explor Factor Anal. 2005;10(7):11.

    Google Scholar 

  57. Revelle W. Personality Project: The psych package for personality and psychological research [Internet]. 2023. https://personality-project.org/r/psych/.

  58. Schendl A, Park G, Xu Z. The spatial prevalence and associated factors of opioid overdose mortality in Milwaukee County, Wisconsin (2003–2018). Spat Spatio-Temporal Epidemiol. 2022;43:100535.

    Article  Google Scholar 

  59. Yang TC, Shoff C, Choi SWE, Sun F. Multiscale dimensions of county-level disparities in opioid use disorder rates among older Medicare beneficiaries. Front Public Health. 2022;10:993507.

    Article  PubMed  PubMed Central  Google Scholar 

  60. Carter JG, Mohler G, Ray B. Spatial concentration of opioid overdose deaths in Indianapolis: an application of the Law of Crime Concentration at Place to a Public Health Epidemic. J Contemp Crim Justice. 2019;35(2):161–85.

    Article  Google Scholar 

  61. Kolak MA, Chen YT, Joyce S, Ellis K, Defever K, McLuckie C, et al. Rural risk environments, opioid-related overdose, and infectious diseases: a multidimensional, spatial perspective. Int J Drug Policy. 2020;85:102727.

    Article  PubMed  PubMed Central  Google Scholar 

  62. Fink DS, Keyes KM, Branas C, Cerdá M, Gruenwald P, Hasin D. Understanding the differential effect of local socio-economic conditions on the relation between prescription opioid supply and drug overdose deaths in US counties. Addiction. 2023;118(6):1072–82.

    Article  PubMed  PubMed Central  Google Scholar 

  63. Marshall JR, Gassner SF, Anderson CL, Cooper RJ, Lotfipour S, Chakravarthy B. Socioeconomic and geographical disparities in prescription and illicit opioid-related overdose deaths in Orange County, California, from 2010–2014. Subst Abuse. 2019;40(1):80–6.

    Article  Google Scholar 

  64. Kurani S, McCoy RG, Inselman J, Jeffery MM, Chawla S, Finney Rutten LJ, et al. Place, poverty and prescriptions: a cross-sectional study using Area Deprivation Index to assess opioid use and drug-poisoning mortality in the USA from 2012 to 2017. BMJ Open. 2020;10(5):e035376.

    Article  PubMed  PubMed Central  Google Scholar 

  65. Stopka TJ, Larochelle MR, Li X, Bernson D, Li W, Ackerson LK, et al. Opioid-related mortality: dynamic temporal and spatial trends by drug type and demographic subpopulations, Massachusetts, 2005–2021. Drug Alcohol Depend. 2023;246:109836.

    Article  PubMed  PubMed Central  Google Scholar 

  66. Bauer C, Hassan GH, Bayly R, Cordes J, Bernson D, Woods C, et al. Trends in Fatal Opioid-related overdose in American Indian and Alaska native communities, 1999–2021. Am J Prev Med. 2024;66(6):927–35.

    Article  PubMed  Google Scholar 

  67. Monnat SM. The contributions of socioeconomic and opioid supply factors to U.S. drug mortality rates: urban-rural and within-rural differences. J Rural Stud. 2019;68:319–35.

    Article  Google Scholar 

  68. Massachusetts Department of Public Health. Data Brief: Opioid-Related Overdose Deaths among Massachusetts Residents, 2013–2021 [Internet]. 2023 Dec. https://www.mass.gov/doc/opioid-related-overdose-deaths-among-ma-residents-december-2023/download.

  69. Plunk AD, Grucza RA, Peglow SL. It is past time to think more inclusively about deaths of despair. Am J Bioeth. 2018;18(10):29–31.

    Article  PubMed  Google Scholar 

  70. Peters DJ, Monnat SM, Hochstetler AL, Berg MT. The Opioid Hydra: understanding overdose mortality epidemics and syndemics across the rural-urban Continuum. Rural Sociol. 2020;85(3):589–622.

    Article  PubMed  Google Scholar 

  71. Case A, Deaton A. Rising morbidity and mortality in midlife among white non-hispanic americans in the 21st century. Proc Natl Acad Sci. 2015;112(49):15078–83.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  72. Venkataramani AS, Bair EF, O’Brien RL, Tsai AC. Association between Automotive Assembly Plant Closures and Opioid Overdose Mortality in the United States: a difference-in-differences analysis. JAMA Intern Med. 2020;180(2):254.

    Article  PubMed  Google Scholar 

  73. Langabeer JR, Chambers KA, Cardenas-Turanzas M, Champagne-Langabeer T. County-level factors underlying opioid mortality in the United States. Subst Abuse. 2022;43(1):76–82.

    Article  CAS  Google Scholar 

  74. Larochelle MR, Slavova S, Root ED, Feaster DJ, Ward PJ, Selk SC et al. Disparities in Opioid Overdose Death trends by Race/Ethnicity, 2018–2019, from the HEALing communities Study. Am J Public Health. 2021;e1–4.

  75. Furr-Holden D, Milam AJ, Wang L, Sadler R. African americans now outpace whites in opioid-involved overdose deaths: a comparison of temporal trends from 1999 to 2018. Addiction. 2021;116(3):677–83.

    Article  PubMed  Google Scholar 

  76. Supply limitations for. opiate prescriptions; exception for palliative care [Internet]. M.G.L., 19D 2019. https://malegislature.gov/Laws/GeneralLaws/PartI/TitleXV/Chapter94C/Section19D.

Download references

Acknowledgements

The authors are grateful to the reviewers and editors for their careful reading of this article and their constructive feedback on the content as well as their editorial suggestions on language. We are especially grateful for their comments on the statistical and spatial analysis methods.

Funding

This study was part of a larger project, Assessing High Risk Opioid Prescribers and Fatal and Non-Fatal Opioid Overdose, funded by the CDC and MDPH. Grant ID: INTF7311HH2 500224100 (PI: Stopka).

Author information

Authors and Affiliations

  1. Department of Urban and Environmental Policy and Planning, Tufts University, Medford, MA, USA

    Sumeeta Srinivasan

  2. Department of Public Health and Community Medicine, Tufts University School of Medicine, Boston, MA, USA

    Jennifer Pustz & Thomas J. Stopka

  3. Institute for Health Metrics and Evaluation, Seattle, WA, USA

    Elizabeth Marsh

  4. Prescription Monitoring Program, Massachusetts Department of Public Health, Boston, MA, USA

    Leonard D. Young

Authors
  1. Sumeeta Srinivasan
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Jennifer Pustz
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Elizabeth Marsh
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Leonard D. Young
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Thomas J. Stopka
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

SS and TJS designed the study. SS, LY and TJS reviewed the data. SS conducted the statistical analysis. JP and SS prepared the figures. SS, JP and EM collated the literature. All authors interpreted the analyses and wrote the manuscript.

Corresponding author

Correspondence to Sumeeta Srinivasan.

Ethics declarations

Ethics approval and consent to participate

This study, focused on secondary analysis of existing data, was reviewed by the Tufts University Health Sciences IRB and designated as not human subjects research.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Srinivasan, S., Pustz, J., Marsh, E. et al. Risk factors for persistent fatal opioid-involved overdose clusters in Massachusetts 2011–2021: a spatial statistical analysis with socio-economic, accessibility, and prescription factors. BMC Public Health 24, 1893 (2024). https://doi.org/10.1186/s12889-024-19399-5

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12889-024-19399-5

Keywords

BMC Public Health

ISSN: 1471-2458

Contact us