Patterns and risk factors of opioid-suspected EMS overdose in Houston metropolitan area, 2015-2019: A Bayesian spatiotemporal analysis

Background Opioid-related overdose deaths are the top accidental cause of death in the United States, and development of regional strategies to address this epidemic should begin with a better understanding of where and when overdoses are occurring. Methods and findings In this study, we relied on emergency medical services data to investigate the geographical and temporal patterns in opioid-suspected overdose incidents in one of the largest and most ethnically diverse metropolitan areas (Houston Texas). Using a cross sectional design and Bayesian spatiotemporal models, we identified zip code areas with excessive opioid-suspected incidents, and assessed how the incidence risks were associated with zip code level socioeconomic characteristics. Our analysis suggested that opioid-suspected overdose incidents were particularly high in multiple zip codes, primarily south and central within the city. Zip codes with high percentage of renters had higher overdose relative risk (RR = 1.03; 95% CI: [1.01, 1.04]), while crowded housing and larger proportion of white citizens had lower relative risks (RR = 0.9; 95% CI: [0.84, 0.96], RR = 0.97, 95% CI: [0.95, 0.99], respectively). Conclusions Our analysis illustrated the utility of Bayesian spatiotemporal models in assisting the development of targeted community strategies for local prevention and harm reduction efforts.


Methods and findings
In this study, we relied on emergency medical services data to investigate the geographical and temporal patterns in opioid-suspected overdose incidents in one of the largest and most ethnically diverse metropolitan areas (Houston Texas). Using a cross sectional design and Bayesian spatiotemporal models, we identified zip code areas with excessive opioid-suspected incidents, and assessed how the incidence risks were associated with zip code level socioeconomic characteristics. Our analysis suggested that opioid-suspected overdose incidents were particularly high in multiple zip codes, primarily south and central within the city. Zip codes with high percentage of renters had higher overdose relative risk (RR = 1.03; 95% CI: [1.01, 1.04]), while crowded housing and larger proportion of white citizens had lower relative risks (RR = 0.9; 95% CI: [0.84, 0.96], RR = 0.97, 95% CI: [0.95, 0.99], respectively).

Conclusions
Our analysis illustrated the utility of Bayesian spatiotemporal models in assisting the development of targeted community strategies for local prevention and harm reduction efforts.

Introduction
Although slightly decreasing from 2017 to 2018, opioid-related overdose remains a leading cause of injury-related mortality in the US, with nearly 70% of drug overdoses involving opioids [1]. According to the US Centers for Disease Control and Prevention (CDC), while deaths rates involving all opioids decreased by 2%, rates involving synthetic opioids, such as fentanyl, saw a 10% increase [1]. Furthermore, recent geospatial analysis on opioid-related overdose has shown a geographical shift from rural areas to large metropolitan areas [2]. Houston, Texas is the fourth largest and one of the most ethnically diverse counties in the US. The 2018 census projections place the total population of the city at 2.3 million [3] and the Houston metropolitan area at 6.9 million people with a growth rate of 2.02% [4]. Like most of the country, Texas has seen a significant increase in the number of overdose deaths involving synthetic opioids. However, deaths involving heroin in Texas have more than doubled and prescription opioids have remained stable [5].
Monitoring of real-time drug overdose patterns is an essential tool in combatting the epidemic; yet, federal and state surveillance systems have considerable limitations. National drug surveillance systems which rely on a top-down approach, are slow to report patterns, and can even lag behind changing overdose trends [6]. Additionally, state and local level EMS data collection systems are not standardized. EMS agencies differ drastically in resources and local EMS agencies with limited resources may be reluctant to invest in more sophisticated data collection software. Due to the variability in collection practices, researchers often rely on field notes and data captured by dispatch systems. Similarly, data from the Survey on Drug Use and Health (SAMHSA) is limited to self-report data and may underrepresent certain populations, such as low-income individuals, ethnic minorities [7], homeless, non-English speaking and technology-impaired individuals [8]. Reliance solely on emergency department (ED) reporting of opioid-related overdose events may also underestimate incidents, since overdose survivors often refuse ambulance transport [9]. Finally, the National Vital Statistics System relies on local coroner and emergency department reported causes of mortality, which may be inconsistent by locality and underreported [10].
Local population-level data, when collected and analyzed in a timely manner, can fill in the gaps left by other surveillance programs. Emergency medical services (EMS) and first responders are critical parts of the emergency care system in the US and the first phase of emergency care. There are more than 20 million EMS transports each year, and emergency 9-1-1 services offer immediate access to an operator who can provide basic life support coaching until help arrives on scene [11]. In most urban and suburban areas, EMS is a component of local fire departments (FDs) [12] and embedded into local communities where they are strategically positioned to engage with residents. Community FDs often serve as a nerve center for public health and safety education, from home safety [13] to emergency first aid [14]. Because of their close ties to the community, local EMS are the first to respond to medical emergencies and therefore the first to notice new trends, for instance, changes in opioid overdose patterns. Early recognition of these hot spots is crucial for timely and targeted intervention and investigation over time, shed light on changes in opioid use and demographics, and projections of future hot spots. The use of Bayesian spatiotemporal models, which have recently gained popularity in research on opioid-related overdose and mortality [15], could be used to provide insight into the spatiotemporal risk of opioid-related incidents and hence guide subsequent policymaking and intervention design. Model based risk estimates offer advantages over the observed risks in that the model-based approaches examine a phenomenon that exists in a particular place and point in time by observing variables that are shared by geographical locations in addition to individual characteristics. The observed EMS data often present a lot of noise, particularly in less populated areas. Model-based approaches allow us to tease out the noise and reveal the true underlying pattern, and to identify "hot spots" of disease incidence and forecast risks. In addition, the model-based approaches allow the investigation of the potential contributing factors to the "hotspots", that can inform policy and public health intervention, and aid the resource allocation and distribution of state and federal funds to combat the opioid crisis.

Study design
We used a retrospective cross sectional study design of prehospital data between January 1, 2015 and December 31, 2019, and mapped suspected opioid-overdose events to zip codes in the Houston metropolitan area. The study population was approximately 24% white, and 48% Hispanic. We assessed the geographical and temporal patterns in opioid-suspected overdose incidents. We developed Bayesian spatiotemporal models to identify zip code with excessive opioid-suspected incidents (i.e., hot spots), and assessed how the incidence risks were associated with zip code level socioeconomic characteristics.

Data sources
We relied on de-identified incident location data extracted from the patient care record system of the Houston Fire Department database, which was provided under a business associates agreement between the University of Texas Health Science Center at Houston and the fire department. We defined the opioid-suspected overdose as an incident which involves administration of naloxone, an overdose reversal medication that is effective for opioids [16]. For the spatiotemporal analysis, we included a total of 2630 EMS calls on opioid-suspected overdose incidents in 84 zip codes concentrated in the densely populated inner-loop area in Houston, TX (S1 Fig

Demographic and socioeconomic variables
Zip code demographic and socioeconomic status (SES) variables were obtained from the US Census Bureau American Community Survey (ACS) 5-year estimates. ACS data were only available for years between 2015 and 2018, and so we used the data from 2018 for 2019. These variables included total population, race/ethnicity (% of White and % of Hispanic), employment (% unemployed), poverty (% living under poverty), education level (% with bachelor degree or higher), income (in dollars), insurance (% of population uninsured), occupation (% of population working as blue collar) and housing (% of renters and % of population living with crowded housing). The complete ACS variables and their summary statistics are shown in Table 1, for the 84 zip codes in Houston metropolitan areas included in this analysis. From 2015 to 2018, we observed a decrease from 5.54% to 4.12% in the average percentage of unemployment, a decrease from 19.1% to 17.3% in the average percentage of living with poverty, a decrease from 25.8% to 22.5% in the average percentage of population with no health insurance, and an increase from $29,700 to $32,700 in average per capita income. However, substantial spatial variation of these demographic and SES variables was observed at these zip codes, see

Bayesian spatiotemporal models
We followed the general inseparable Bayesian spatiotemporal modeling framework by first assuming the appropriate distribution of the observed data, then assign the structures to the spatial and temporal components. Specifically, let Y it denoted the observed opioid-suspect EMS overdose counts in year t and zip code i, we assumed a Poisson distribution with relative SIRs have a straightforward interpretation, where value greater than 1 indicated more observed counts than the expected, and hence could be defined as high incidence areas or "hot spots." SIR with a value less than 1 indicated observed counts fewer than the expected, and therefore "cold spots." SIR of value 1 indicated equal observed and the expected counts. Second, we could compare observed SIR and model based relative risk (RR), as RR provided a smoothed version of SIR and was generally considered more accurate. In the Poisson model above, we further decomposed the relative risk into the following components where α was the overall intercept, and x it is the vector of zip code level covariates with coefficient vector β. We included a random effect φ i to account for the spatial correlation, and a random effect δ it to account for the spatiotemporal interaction. In Bayesian spatiotemporal models, there random effects were assumed to have different distributions as priors. Here we considered several popular spatial models for random effects φ i , and compared model performance using three model selection criteria: deviance information criterion (DIC), widely applicable Bayesian information criterion (WAIC) and conditional predictive ordinate (CPO). In all three criteria, models with the smaller values of these measures were considered better [17]. The spatial models we considered here included the Besag, York, and Molliè (BYM) model [18], the Besag's proper spatial model [19] and Leroux model [20]. Details of the models can be found in their original papers but we provided a brief summary here. The BYM model considered the random effect φ i the summation of a spatially structured term u i and a spatially unstructured term v i . The spatially structured term was assumed to have an intrinsic autoregressive (ICAR) model, where the spatial dependency between area i and j was encoded in the adjacency matrix W. The entries in W were assigned value 1 if two areas shared the same geographical boundary, and 0 otherwise. The spatial unstructured term was assumed to be independently and identically distributed. The ICAR model was also called "Besag improper" model as its joint distribution was improper (i.e., the precision matrix was singular). The Besag proper model overcame this issue by constructing a non-singular precision matrix. The third model considered here was the Leroux model where the conditional distribution of the spatial random effect was The parameter ρ quantified the degree of spatial correlation of φ i , with ρ = 0 corresponding to independence and ρ = 1 representing strong spatial correlation throughout the region. Residual variation not explained by spatial correlation was captured by the variance parameter σ 2 . We used Integrated Nested Laplace Approximation (INLA) for Bayesian inference as an alternative to the Markov Chain Monte Carlo (MCMC), for its fast computation [21]. All analyses were performed in R using package INLA (R Studio, Boston MA).

Results
EMS calls for opioid-suspected overdose incidence increased from 449 in 2015 to 666 in 2019 (Fig 2). We observed a slight decrease in 2017 (n = 473) compared to 2016 (n = 487), with a significant increase in calls in 2018 (n = 555) and again in 2019. In addition, we also investigated the temporal trend of incidents by month, day of the week and time of the day. Opioid-suspected overdose incidents occurred more frequent during the summer time between May and September, on Saturdays, and around the time between 7pm and 11pm. The observed SIRs of opioid-suspected overdose EMS incidents are shown in Fig 3, where we clearly observed the spatial patterns. Darker shades indicated areas with higher risks, primarily in the south and central area of Houston. We also observed substantial changes in the spatial variation over time, which provided strong evidence to include spatiotemporal interaction in the Bayesian spatiotemporal models. Without the interaction term, we essentially assume the same spatial pattern and hotspots every year. With the interaction term in the model, we capture the changing spatial pattern over time and allow different hotspots every year. The analysis results from three Bayesian spatiotemporal models with different spatial structures were presented in Table 2. All model selection criteria of DIC, WAIC and CPO pointed the BYM model as the best fitting model, and therefore we made inference based on this model. Of all the SES variables investigated, we found that zip code neighborhoods with higher percentage of renters were associated with a small increase in the likelihood of higher opioid-suspected EMS relative risk (RR = 1.03; 95% CI: [1.01, 1.04]), while neighborhoods with more white population were associated with a small decease in the likelihood (RR = 0.97; 95% CI: [0.95, 0.99]). Neighborhoods with crowded housing, defined as more than one person occupying a room, appeared to have lower relative risk of opioid-suspected EMS incidents (RR = 0.9; 95% CI: [0.84, 0.96]) No other SES variables were found statistically significant. Model-based zip code level RR from the BYM model was presented in Fig 3. The spatial patterns were largely similar compared to the SIR, with elevated risks in the south part of downtown areas throughout the years, and

Discussion
Our study found that the number of EMS calls for opioid-suspected overdose significantly increased in the Houston metropolitan area between 2015 and 2019, and this trend is continuing into 2020 with 530 suspected overdose calls from January to June; a 16% increase over the same time last year. Across the nation, opioid overdose rates have seen a 54% increase in major metropolitan areas from 2016 to 2017 [1]. Originally fueled by prescription opioids, recent rises in overdoses are now driven by heroin and fentanyl, which is making the urban drug market a hotbed of overdose mortality. Applying spatiotemporal modeling to overdose incidence data, as demonstrated in this analysis can help communities struggling with overdoses, particularly those in rural areas, forecast overdose trends and develop a targeted approach to early intervention and prevention efforts. Additionally, temporal trends can be tracked and help inform staffing practices of EMS providers. Our analysis revealed more overdoses occur in the summer months, between May and September and between 7pm and 11pm. This information, combined with geospatial trends can assist municipalities with hiring and staffing practices, ensuring enough adequately trained personnel are on duty during those high-risk times.
Community strategies for understanding their own temporal and spatial patterns could be useful in positioning resources and developing neighborhood programs. The use of real-time surveillance data, from hospitals or EMS, offers specific advantages over reliance solely on death data. For example, discrepancies in reporting by state coroners was estimated in one study to that nearly 70,000 opioid-related deaths may have gone unreported since 1999 [22]. Because national statistics on drug intoxication deaths rely on data derived from state death certificates, inconsistencies in reporting from state to state can have serious implications since codes used to report causes of death are applied differently from state to state. Additionally, rural counties may lack resources for investigating drug overdose deaths and may over-use a general "presence of unspecified drugs, medicaments, and biological substances" code for reporting opioid-related overdose deaths. Similarly, overreliance on ED overdose data can similarly under-estimate overdose events, since many overdoses do not result in an ED visit. Not every 911-call results in transport to the ED [16] as survivors can refuse transport. Research has revealed a significant portion of people who inject drugs refuse EMS transport, citing fear of harassment and discrimination, potential arrest, and anticipated financial costs [23,24]. Similarly, Good Samaritan laws have increased bystander administration of naloxone to suspected opioid overdose victims and although recommended, not every resuscitation is followed by a call to 911 [25,26]. Although it is difficult to track, several studies have reported large percentages of overdose survivors do not get treated by EMS [27]. An estimated 25 to 50 nonfatal overdose occur for every overdose death [28] and EMS call surveillance can provide the data necessary to fill in the data gaps in identifying overdose trends.
Some communities have implemented programs that utilize EMS incident data for targeted interventions. A local Fire Department/EMS and law enforcement collaboration, in Columbus, Ohio provides free transportation to treatment facilities for overdose survivors who decline transport to the ED [29]. The REACT program also provides harm-reduction and education outreach to communities that are hard-hit by the opioid crisis [29]. In San Francisco, the DOPE program also focuses on harm reduction and conducts community outreach, using an EMS early warning system that monitors call data for any surges in overdose events [30]. Analysis of DOPE has found the program sites tend to be located in densely populated minority communities with significant economic disparity, and who experience more overdose mortality. Furthermore, the program has demonstrated marked success in its efforts to reduce opioid-related overdose deaths [30]. A program in Houston, Texas, the Houston Emergency Opioid Response System (HEROES) uses information from EMS runs involving naloxone resuscitation to deploy a mobile response team including a peer recovery coach and licensed paramedic to engage overdose survivors into treatment [31]. Geospatial analysis of the most current data on suspected opioid-overdose calls can inform community programs on trending hot spots but can also help target specific populations that are experiencing increased overdose events by including certain demographic characteristics in the analysis. Interventions that target special populations, such as youth, females of childbearing age, and elder populations can use geospatial analysis to narrow their target areas. Spatiotemporal modeling provides advantages over monitoring real-time EMS call data. First, it helps visualize a complex and dynamic problem that fluctuates across time and space. Second, geospatial analysis can actualize changes across time spans, from days to years and can identify patterns to assist with causation research. Spatiotemporal models can help inform policy makers with budget allocation and staffing concerns. Bayesian spatiotemporal modeling approaches can examine existing trends and apply them to forecast patterns, which is useful for prevention efforts. Altogether, Bayesian spatiotemporal modeling can help EMS providers, researchers, and lawmakers refine programs and operate in a constantly changing environment.
Future policies and funding should consider incorporating geospatial modeling of first responder and other data sources within each community, to provide a more comprehensive response to the current opioid epidemic.

Limitations
This study has several limitations. First, we focused here on one community and one major EMS provider, primarily since national EMS data at the location level do not exist. Houston, Texas is the fourth largest and most diverse city in the United States [3]; however, it may not be representative of every major metropolitan area in the US that is experiencing an influx of opioid-related overdoses, nor do our results apply to rural areas which may have different demographics and risk factors. Additionally, as with most EMS data, not all overdoses included here were necessarily opioid-related, although they were all suspected opioid overdoses indicated by the use of naloxone. EMS systems obtained data are subject to distinct limitations of convenience sampling and missing data can cause calculation errors even with imputation [32]. Finally, EMS data is not completely representative of all overdoses in a community, since patients may not always call 911 for emergency assistance.

Conclusions
We found that several popular Bayesian spatiotemporal models were useful in identifying the socioeconomic covariates that were potentially associated with the high EMS incidence risk. We noted clear patterns emerging geographically, as well as over time. Areas with high rates of renters had higher overdose incidence risk, while crowded housing and larger proportion of white citizens had lower rates. Bayesian models could be useful to develop targeted community strategies for local prevention and harm reduction efforts. Langabeer.