A large-scale retrospective study of opioid poisoning in New York State with implications for targeted interventions

Opioid overdose related deaths have increased dramatically in recent years. Combating the opioid epidemic requires better understanding of the epidemiology of opioid poisoning (OP). To discover trends and patterns of opioid poisoning and the demographic and regional disparities, we analyzed large scale patient visits data in New York State (NYS). Demographic, spatial, temporal and correlation analyses were performed for all OP patients extracted from the claims data in the New York Statewide Planning and Research Cooperative System (SPARCS) from 2010 to 2016, along with Decennial US Census and American Community Survey zip code level data. 58,481 patients with at least one OP diagnosis and a valid NYS zip code address were included. Main outcome and measures include OP patient counts and rates per 100,000 population, patient level factors (gender, age, race and ethnicity, residential zip code), and zip code level social demographic factors. The results showed that the OP rate increased by 364.6%, and by 741.5% for the age group > 65 years. There were wide disparities among groups by race and ethnicity on rates and age distributions of OP. Heroin and non-heroin based OP rates demonstrated distinct temporal trends as well as major geospatial variation. The findings highlighted strong demographic disparity of OP patients, evolving patterns and substantial geospatial variation.

monitor the opioid epidemic at the community or patient level. This depends on critical knowledge to answer many questions. For example, which regions or communities have most serious opioid overdose problem and what are the latest trends? What population groups have a higher overdose risk?
As recommended by U.S. Department of Health and Human Services, better data and reporting can help us better understand the crisis 17 . With increased accessibility of health data driven by open data initiatives, largescale patient level data analysis provides an opportunity for a data-driven approach to identity patterns of the opioid epidemic and to discover potential causes of opioid-related deaths through studying diagnoses from hospital visits at large scale. The NYS Statewide Planning and Research Cooperative System (SPARCS) 18 collects patient level details on patient characteristics, diagnoses and treatments, services, and charges for each inpatient stay and outpatient visit for emergency department, ambulatory surgery, and outpatient services. It also includes locations of patients (street addresses). Such unprecedented access to large numbers of patient records linked with spatial data allows researchers to explore opioid poisoning in NYS with significant improvement of accuracy and coverage compared to prior studies on non-medical opioid use at a national, state or urban level 19,20 . However, current state-wide opioid reports based on SPARCS are only limited to county level.
Data driven studies from this approach could provide major impact for improving public health studies. Our work on statewide data analytics offers a critical step for community level interventions to prevent or reduce non-medical opioid use and related opioid overdose by identifying high risk communities within which to apply more specific countermeasures. Analyzing opioid data at fine spatial resolution such as zip code can also reveal community level distribution patterns of demographic, regional and temporal trends, which can provide crucial knowledge for local governments, residents and health systems and service professionals seeking effective solutions to the opioid crisis 21 .
Aiming to better manage the opioid-related public health crisis, we investigated demographic disparities, temporal and geographic patterns of hospital, emergency, ambulatory surgery and outpatient visits (2010-2016) for NYS resulting in at least one OP discharge diagnosis in SPARCS. We performed the following studies to assess: (1) demographics to find disparities of OP rates between age, gender, and race/ethnicity groups; (2) temporal trends of OP rates and any disparities between demographic groups; (3) heroin vs. non-heroin poisoning to understand the differences between two major components of the opioid epidemic; (4) geographic trends of OP rates across NYS at zip code level; and (5) correlations between demographic and socio-economic factors that may be linked to the higher risk of OP.

Methods
This study was approved by the Stony Brook University Institutional Review Board and the Office of Quality and Patient Safety, Department of Health of NYS, and all methods were performed in accordance with the relevant guidelines and regulation. Informed consent was not needed as the study had no contact with participants, which was determined by the Stony Brook University Institutional Review Board. The data were obtained from a NYS administrative database. All study analyses were completed between 2017 and 2019.
Study populations and data sources. In this retrospective cohort study, we used data from SPARCS, an all-payer health claims database for NYS, required for all Article 28 licensed facilities 18 . We identified OP patients based on ICD codes (primary and secondary diagnosis codes, ICD-9 from January 1, 2010 to September 30, 2015, and ICD-10 from October 1, 2015 to December 31, 2016). We will continue update our analyses when new SPARCS data beyond 2016 become available. Only NYS patients (with a valid NYS 5-digital zip code) with at least one OP related diagnosis were included in analyses. The identified data set included patient level information, such as demographics, diagnoses, treatments, services, addresses of patients, among others. The racial and ethnic categories utilized in the analysis (in Table 1) is consistent with SPARCS methodology, which included non-Hispanic Whites, African-American, Hispanic, non-Hispanic Asian categories.
Statistical and spatial analyses. All historical records per individual were aggregated using an encrypted unique patient identifier and ordered by the corresponding hospital visiting dates. For analyzing spatial patterns, the patients were also aggregated by patient residential zip code, rather than facility zip code.
For temporal trends, if a patient had multiple OP-related clinical contacts in a year, only the first recorded episode was included. In an analysis combining 7 years' data, each unique patient was only considered once at their index OP event even if the patient had had multiple OP episodes over more than one year.
OP rates overall and within each sociodemographic category were estimated and compared across 7 years using logistic regression. To evaluate whether OP rates increase over the years, we performed contrast F tests based on the logistic regression. For correlation studies at zip code level, we used TIGER/Line with Selected Demographic and Economic Data from the US Census Bureau 24 . To evaluate the correlations between sociodemographic factors with OP rates at zip code level, multiple linear regression models were fit for OP rates as the dependent variable. Zip code level factors, i.e. median household income, race and ethnicity, gender, % of high school education, house units and population density (persons per square mile) were included as covariates. Regression coefficients were standardized in order to compare the effects of factors on different scales.
All maps were generated using ArcGIS Desktop 10.5 (Esri, Redlands CA). Statistical analyses were performed using SAS version 9.4 (SAS Institute, Inc., Cary, NC).

Results
Among 146,598,009 encounters of 26,413,181 patients visiting at any facility (inpatient, outpatient, emergency room, ambulatory surgery) in SPARCS, 58,481 patients with an OP-related diagnosis on 72,102 hospital visits between 2010 and 2016 were identified in NYS through the residential zip codes. The overall OP rate of NYS was 296.2 per 100,000 persons over the seven years ( Table 1). The counts and rates of OP patients steadily increased during 2010-2016 within each demographic group (all P < 0.0001) and for NYS overall (P < 0.0001). Each increase overall and within each demographic group demonstrated significant linear trends (all P < 0.0001). The comparison between different gender, age, race and ethnicity groups also suggested demographic disparities.
Demographic-based analysis. There were strong disparities among racial/ethnical groups. Figure 1A showed a significantly higher rate of 352.1 per 100,000 among non-Hispanic whites with a significantly lower rate of 31.9 among non-Hispanic Asians (overall P < 0.0001). Among different age groups (categorical: 0 to < 85 in 5-year increments, and ≥ 85) for all OP patients in 2010-2016, there was a high peak among young adults (ages in 20-30) and two local peaks among middle ages (ages in 50-60) and younger children (ages < 5) (Fig. 1B). The age distributions showed significantly different patterns of OP rates between racial/ethnical groups (overall P < 0.0001). Non-Hispanic white OP patients were more concentrated among young adults (20-35 age groups) (Fig. 1C) whereas the Non-Hispanic African-American and Non-Hispanic Asian OP patients had a higher concentration among 50-70 age groups (Fig. 1D) and above 75 age groups (Fig. 1F) respectively. There was no obvious peak for the age distributions of Hispanic OP patients with the relatively highest rates among the 45-60 age groups (Fig. 1E). The OP rates were not uniformly distributed over age groups within each racial/ ethnical group (all p values of goodness of fit test < 0.0001, Figs. 1C, 1D, 1E and 1F).
Temporal trends of demographic characteristics. The OP rates for all demographic groups increased dramatically after 2014 (Fig. 2). Across years 2010 to 2016, males always had a higher rate than females (all P < 0.0001) ( Fig. 2A). Non-Hispanic whites and non-Hispanic African-Americans increased at a higher rate than Hispanics and non-Hispanic Asians (  www.nature.com/scientificreports/ patients with heroin poisoning and those with non-heroin poisoning (Kolmogorov-Smirnov test P < 0.0001) (Fig. 3B). Young adults within 20-35 age range had the highest OP rates attributable to heroin. Older OP patients were more likely to overdose on non-heroin opioids. Heroin and Non-heroin based OP groups showed distinct temporal trends (Fig. 3C).
Year 2015 saw a dramatic increase in the Non-heroin poisoning rate and by 2016, the rate was approximately 4 times that of OP attributable to Heroin (P < 0.0001). The rate of heroin poisoning was accelerating in males from 2010 to 2016 (Fig. 3D) and was significantly higher than that of females overall all years (P < 0.0001). The rate of Non-heroin poisoning among females began to surpass that of males after 2015 (P for 2015 = 0.51, P for 2016 < 0.0001, Fig. 3E). The OP rates for Heroin dramatically and consistently accelerated among non-Hispanic whites from 2010 to 2016 (Fig. 3F) and is significantly higher than that of non-Hispanic African American and Hispanic (P < 0.0001). The Non-heroin OP rate began to increase after 2014 among all race and ethnicity groups (all P < 0.0001) (Fig. 3G). For Heroin poisoning, the rate among young adults (25-34) had the most significant increase (Fig. 3H). For non-heroin poisoning, the rate of patients among elders (55-64 and 65 +) had most significant increase (both P < 0.0001) (Fig. 3I). Figure 4 presents the maps for both OP rates (Fig. 4A) and counts (Fig. 4B) in NYS in 2016. While major cities (New York, Buffalo, Rochester, Albany) had elevated counts of OP patients (Fig. 4B), New York City (NYC, including the Boroughs of Manhattan, Bronx, Staten Island, Queens and Brooklyn) had lower OP rates. Long Island (Nassau and Suffolk Counties) had both comparatively high counts and rates. Heroin and Non-heroin based OP mappings demonstrated www.nature.com/scientificreports/ distinct spatial patterns. Most big cities had higher counts of both heroin and non-heroin OP patients (Fig. 4D,  4F). Higher Non-heroin OP rates (Fig. 4E) were not limited to big cities and appeared across the whole state. An animation video (Appendix Video 1) in the Supplemental Content demonstrates the evolving spatial-temporal patterns of OP rates in NYS.

Sociodemographic factors associated with OP.
We studied potential correlations between OP rates and demographic and social-economic factors by analyzing a group of demographic and socio-economic factors at the zip code level. Among 8 factors, no high association was found (Supplemental Content, Appendix Table 1). The median household income had a significant but weak negative correlation relationship with overall, heroin and non-heroin OP rates (standardized regression coefficients are − 0.235, − 0.096 and − 0.252 respectively, both P < 0.0001), suggesting zip codes with higher OP incidence rate had a lower median household income. Population density (person per square mile) also had a significant but weak negative correlation relationship with heroin OP rates (standardized regression coefficient − 0.139, P < 0.0001), suggesting zip codes with higher OP incidence rate had a lower population density.

Discussion
This study provided a large-scale patient level analysis of demographic disparities, spatial patterns, temporal patterns and sociodemographic factors for OP in NYS. With state-wide data including address information, we had sufficient sample size to conduct patient level evaluation at a finer, zip code level of geographic resolution. This study should provide essential new knowledge for local governments and healthcare institutions for design of more focused and targeted interventions. An important finding was that between 2010 and 2016, the overall rate of OP in NYS increased for patients > 65 years by 741.5% compared to the overall 364.6% increase in the general OP rate. This jump is due to the OP rate of 20 per 100,000 in 2010, a relatively low base rate compared to other groups 15-64 years in the same year. Whereas the other racial/ethnic groups had increases of 2.9 to 4 times, Asian, non-Hispanic patients saw a seven-fold increase in the OP rate from 2010 to 2016, with a peak in patients greater than 75. Given that Asian, non-Hispanic patients have a comparatively low prevalence of heroin-associated OP, it appears that patterns of opioid prescribing may be involved in the increased rates of OP in this group from 2010 to 2016. As such, the study reveals an important finding that can be further explored regarding identifying potential factors such as medical comorbidity (e.g., disorders increasing vulnerability to respiratory depression), differential sensitivity to opioid-induced respiratory depression, and prescribing styles of physicians who treat Asian, non-Hispanic patients. However, we also demonstrate that for the overall group > 65 the acceleration in the OP rate is greater than other age groups 15-64 for both Heroin and Non-heroin OP.
Coinciding with the implementation of widespread mitigation measures for the COVID-19 pandemic, there is a pressing need for opiate overdose reversal interventions and effective alternatives to pain management, which involves the need for increasing the capacity of primary care and health care center, including expanding overdose prevention education and naloxone provision [25][26][27] . Examples of evidence-based targeting interventions include the Hub and Spoke models and medication assisted therapy (MAT). Specialized opioid treatment and comprehensive system wide programs such as the Hub and Spoke have been developed in Vermont. Opioid use disorders (OUD) treatment centers are the hub and designated clinics that prescribe buprenorphine are the spokes 28 . Clinic-based MAT and the Project Extension for Community Healthcare Outcomes (ECHO) virtual clinics have successfully utilized an integrated team approach to support opioid management nationwide.
While there have been successes in urban health care centers in the implementation of these programs, there are regional challenges for rural communities and hospitals where OP rates have been steadily increasing 29,30 . These challenges include training personnel to deliver such programs and the shortage of trained and waivered buprenorphine prescribers, which has been exacerbated during the COVID-19 pandemic. An integral part of MAT is the education of patients and community members on opioid overdose prevention, stabilization on opioid antagonist or agonist medication, as well as naloxone distribution in high need communities with the highest rates of heroin related deaths and opioid overdoses reported in this study 31 . Limitations of this study. There are some limitations inherent to the SPARCS all-payer database from which our analysis derives. Most importantly, there is the potential for systematic over-and under-counting the number of unique individuals who had an episode of opioid poisoning during the look-back period from 2010 to 2016. In 2015, the SPARCS diagnosis codes transitioned from ICD-9 to ICD-10 and might account for a portion of the increased rates observed afterwards 32 . In particular, the well documented rise in rates of OP related to synthetic opioids such as fentanyl was not captured directly prior to October 2015 as fentanyl and its analogues were not routinely screened for in clinical settings, and the SPARCS identifiers only had a non-specific code of "Poisoning; Other" (Supplemental Content, Appendix Table 2). Additionally, due to the coding variations between ICD-9 and ICD-10, in which unintentional overdose data was not consistently available, the data does not capture patients who intentionally overdosed after misusing opioids. Given the subsequent availability of this data after ICD-10, this data will be available and reported in future studies.
SPARCS also lacks integration with vital records, thus for example, the records of clinical encounters of OP patients who did not survive prior to arrival at a hospital were not included and may have contributed to a systematic undercount. We intended to further specify the study findings in future work pending incorporating a dataset that has been requested from NYS containing death records linked to SPARCS. For individuals who get treated outside of the SPARCS healthcare provider network, OP events may also have been underreported.
The uniqueness of SPARCS data in NYS presented us with an opportunity to advance methodology not quite available in other States. The findings concerning NYS may not generalize to other States with both large www.nature.com/scientificreports/ metropolitan and urban centers, as well as rural settings. ZIP codes are more likely attributed to roads, post offices, and other facilities within the U.S. postal system, which limited the association between OP incidences and the census related demographic or socioeconomic information. SPARCS has not been able to release latest data for research use due to logistics issues and the delay from COVID-19. We will extend our study once SPARCS is ready to release the latest dataset.

Conclusions and relevance.
Increasingly availability of open health data provides unique opportunity for a more precise understanding opioid epidemic at large scale. This study demonstrates major disparities among demographic groups and geospatial regions, and evolving patterns of the opioid epidemic in NYS.