Introduction

Over the last decade, the USA has experienced a dramatic rise in opioid abuse. Opioid overdose deaths have increased from 21,088 in 2010 to 49,860 in 2019 [1]. The opioid crisis has affected both rural and urban areas but is most acute in regions experiencing high levels of poverty and marginalization, where residents have low access to healthcare. Cities that are sites of great socio-spatial inequalities have experienced high rates of opioid overdose deaths (OODs) [2]. However, studies show that OODs are not uniformly distributed within cities, rather they tend to be spatially clustered in neighborhoods experiencing concentrated poverty, marginalization, and deprivation [3–8]. Further, studies report that the current overdose crisis is fueled by numerous demographic and socioeconomic factors [9, 10], where socio-spatial inequities influence substance use disorder. As Sadler and Furr Holden [11] note, built environmental characteristics cause OODs to cluster by way of a deprivation amplification effect. An effective response will require an integrated, data-guided approach that involves healthcare professionals, policymakers, the justice system, and a diverse array of community partners and organizations [12]. The spatial analysis, modeling, and mapping capabilities of geographic information systems (GIS) offer a powerful approach to study the opioid crisis [13–15]. This approach has enabled researchers to examine the effectiveness of local policy responses [11, 16].

As the opioid crisis worsens, a comprehensive and granular framework to guide initiatives at the community level is desperately needed. Since the influence of key factors varies with geographic scales and across communities, relationships that are evident at more global scales (e.g., state) may not appear at the local scale (e.g., neighborhood) [17]. As a result, resources and services may be misdirected and opportunities to identify effective community-level solutions may be missed. Thus, the selection of an appropriate scale of analysis is fundamental to meaningful geographic inquiries into public health crises. Policy responses must also be sensitive to such variations across scales. This is especially true in ethnically diverse and racially segregated cities that show pronounced racial disparities in socioeconomic and health conditions. While the general relationship between socioeconomic status and opioid overdose has been well-documented [17–30], there is a clear need for multiscale geospatial analysis and modeling of relationships and their spatial variations across communities.

This study undertakes a multiscale geospatial analysis of the opioid crisis in Milwaukee County, WI, where 418 overdose deaths were reported in 2019 [31]. The death count has risen to 537 in 2020 [32], resulting in an incidence rate of more than 67.5 per 100,000—well over twice than that of the rest of Wisconsin and one of the highest in the USA. As Wisconsin’s largest and ethnically most diverse urban metropolitan area, Milwaukee provides a unique opportunity to study the intersection among socioeconomic factors, drug overdoses and the deleterious effects of racial segregation [32]. Moreover, the rise in overdose deaths among African Americans while rates of overdose are declining in the White population [33], a trend that is also evident in Milwaukee County [34], suggests that our current response to the opioid crisis has disproportionate benefit for White communities.

This paper makes significant contributions to the study of the opioid crisis. First, we are attentive to the significance of geographic scale and recognize that a complex, diverse urban environment must be examined through a fine-grained multi-scalar approach [35]. Therefore, to explicate Milwaukee’s opioid crisis, we examine it at the county-wide scale, municipal scale, census tract scale, and neighborhood scale. We employ multiscale geographically weighted regression (MGWR) for modeling, an innovative methodological approach that explains how geospatial patterns vary across scales, enabling us to examine the differential influence of factors globally, regionally, and locally [35, 36]. Such a multi-scalar approach is relatively new to opioid studies, which tend to focus on a single geographical scale (state level or county level) as the unit of analysis. Second, to further examine known and discover new OOD determinants, we compiled a comprehensive dataset of 225 candidate variables, taking into account the various community level factors that shape OODs. We were particularly attentive to the influence of racial segregation and structural inequalities on OODs. Third, to analyze the opioid crisis in the context of racial segregation and associated structural inequalities, we divide Milwaukee County into three regions based on the racial composition of census tracts. Consequently, we examine spatial variations in OODs and policy responses across communities specified by racial/ethnic identities. Our findings show that policy interventions had differential impacts on communities based on race/ethnicity. The study therefore provides new empirical and methodological insights that can guide future research and community-targeted practices, policies, and interventions.

Study Area, Materials, and Methods

Located along the western coast of Lake Michigan, Milwaukee is a densely populated county. It has an estimated population of 945,726 and demographically breaks down to 58.9% White, 26.3% Black, and 15.6% Hispanic or Latino origin (of any race) [37]. Milwaukee County is not only the most populous county in Wisconsin but also the most racially and economically segregated [38, 39]. The City of Milwaukee, with an estimated population of 590,157, is located within Milwaukee County. Although crime, poverty, and other social factors often overlap, mapping these factors in Milwaukee shows dramatic divisions that align with the city’s history of racial segregation [40, 41]. Fatal overdose rates do not necessarily follow socioeconomic lines but instead may depend on racial composition, urbanicity, and family organization within a community [42]. Thus, to explicate the opioid crisis in Milwaukee County, a comprehensive dataset that captured neighborhood-level socioeconomic and demographic factors was established.

Data Sources and Methods

Our study drew data from a wide range of community resources and relied on three basic data categories. The first source of data captured the precise location of fatal OODs. The 2019 dataset was provided by the Milwaukee County Medical Examiner office, which investigates all suspicious deaths occurring in Milwaukee County. The second category of data was health-related, with several subtypes. Public health measures were provided by the Wisconsin Department of Health Services. Emergency medical call for service (EMS) data were provided by the Milwaukee Fire Department. Opioid availability data were derived from The Automation of Reports and Consolidated Orders System (ARCOS) dataset from the Drug Enforcement Administration (DEA) [43]. The third category of data was demographic- and socioeconomic-related. We collected the data from US Census Bureau’s 2018 American Community Survey 5-Years Estimates of census tract level data (Data.census.gov). Crime data were provided electronically by the Milwaukee County Police Department for 2019. To consider the effect of digital divide (disparity in access to computer/internet across households), we calculated the “digital divide index” for all census tracts in Wisconsin, using Gallardo’s [44] proposed formula.

Overall, we compiled a dataset of 225 variables, which we identified as explanatory (see Table S1 for full list of variables). All variables were collected, processed, and joined to the administrative boundary shapefile of census tracts collected from the TIGER/Line database (www.census.gov), using ArcGIS Desktop 10.7. Through statistical analysis, we identified the 12 most significant explanatory variables (out of 225) for spatial modeling.

Spatial modeling enables us to statistically investigate the geographic relationships/correlations between explanatory variables and disease outbreak [45–47]. Conventional “global” regression modeling assumes that relationships are constant across a study area; i.e., relationships do not change across space [48]. However, many spatial processes differ with the geographical context and scale. We therefore use MGWR, as it is sensitive to the effects of scale and can more accurately capture spatial heterogeneity, minimize overfitting, mitigate concurvity, and reduce bias in the parameter estimates [35, 36, 49, 50]. To our knowledge, this paper provides the first application of the MGWR framework for modeling OOD determinants. Following [51, 52], composite maps were prepared to visualize the parameter estimates and their statistical significance. The MGWR maps are presented for each variable to investigate the parameter estimate spatial heterogeneity.

Results

As metropolitan counties struggle to respond to the opioid crisis, it has become evident that new strategies are needed to better define the problem and its underlying causes [17, 53]. A significant challenge to addressing the crisis is that influential factors can vary significantly across neighborhoods/communities. As a result, many interventions and policies lack universal effectiveness. This is particularly problematic in cities such as Milwaukee. Indeed, Milwaukee is one of the most ethnically diverse cities in the USA and, at the same time, it is one of the most racially segregated [37, 38]. Milwaukee is also among the US metropolitan areas experiencing the greatest increase in the number of overdose deaths [31]. Thus, Milwaukee County provides a unique opportunity to study the intersection among socioeconomic factors, factors related to race and discrimination, and the opioid epidemic.

GIS-based mapping of data is a powerful approach that permits inference of complex interactions among variables based on their temporal-spatial relationships. Geospatial approaches have been previously used to examine factors that influence the opioid crisis [11, 14, 15, 16]. However, prior approaches have been limited as they are based on assumption that all modeled processes operate at the same spatial scale. This creates interpretational challenges and limits the ability to make accurate inferences about relationships that could guide community-targeted policies and interventions [36]. MGWR allows a more flexible exploration of the relationships at varying spatial scales. To our knowledge, this is the first published multi-scalar investigation of the opioid crisis.

To further examine expected and discover unexpected OOD determinants, a comprehensive dataset of 225 candidate variables was first compiled. The number of variables was narrowed down through a two-step approach. First, we identified the variables that exhibited high collinearity based on their Pearson’s correlation coefficients and global variance inflation factors (VIFs) when evaluated against one another [54]. Second, the most uncorrelated factors were selected as stepwise mixed regression inputs to choose a subset of factors by removing non-significant explanatory variables. Based on these analyses, we derived 12 explanatory variables: “naloxone availability,” “crime rate,” “inequality of household income,” “disability rate,” “resided in the same unit for 20 to 30 years,” “total dispensed opioid 2006–2014,” “resided in the same unit for 4 to 8 years,” “houses built before 1950,” “hospital and clinic accessibility,” “renters spent 30% of income on rent,” “households with young children,” and “population over 25 with college degree.”

Global Model

To provide context for our MGWR results, we conducted global ordinary least squares (OLS) regression on the 12 selected variables. Results for the 12 selected OOD explanatory variables in Milwaukee County’s 296 census tracts, with an R2 of 0.371, are summarized in Table 1.

Table 1 Results from the ordinary least squares regression, regression coefficient, and t-statistics were provided to examine the significance of variables, and VIFs (variance inflation factors) were provided to investigate their dependence

The global model produced a relatively moderate R2, indicating that about 40% of the variation in OODs is accounted for by the selected variables. Multicollinearity was not an issue, as the VIFs for each explanatory variable were under 10 [54]. Based on a standard t-value threshold of 1.96 for a 95% confidence level, each of the variables, except educational attainment (population over 25 years old with a college degree), was statistically significant.

MGWR Results

The global OLS model assumes that relationships are constant across the study area. We employed MGWR to address scalar variations, deter unexplainable high levels of spatial heterogeneity, and control multicollinearity and concurvity in the local data subsets. Calibrating an MGWR model produces a vector of optimal bandwidths that describes the spatial scale at which each process in the model varies [51]. MGWR was applied to the same set of explanatory variables used in the global OLS model. The R2 increased to 0.560 in the MGWR model from 0.371 in the global OLS model; i.e., MGWR performed better than the global OLS model. AIC (Akaike information criterion) is an estimator of prediction error and an indicator of relative quality of statistical models. In this case, the AIC decreased 55 units, from 730 in the global model to 675 in the MGWR model, suggesting less prediction error and higher model quality in MGWR [55].

In our model, we used the number of neighboring census tracts for each census tract (observation) as the bandwidths. In Table 2, bandwidths related to each explanatory variable are listed (the global model assumes a bandwidth of infinity). Eight relationships affected OOD rates at a global scale, with relatively large bandwidths indicating that almost all the data are included in each local subset. Two associations occurred at a regional scale with bandwidths implying several hundred nearest census tract neighbors (less than 200); “naloxone availability” varied locally, with a relatively small bandwidth (less than 50).

Table 2 Results from the MGWR; using bandwidth (neighboring census tracts) and effective parameters, the scale at which different processes operate is determined

Our MGWR model identified eight key factors that influenced OODs in Milwaukee County in 2019 at a global (i.e., county-wide) scale. Of these, “disability,” “neighborhood instability,” and “total dispensed opioid drugs (2006–2014)” are the main variables that had positive associations with OODs across the county. Overall, these findings are consistent with the results from prior research, suggesting that higher rates of overdose are associated with impoverished socioeconomic conditions [11, 56–58], physical pain, and over-prescription of opioid medications for pain [59, 60]. The factors “educational attainment,” “access to healthcare,” “households with children,” and “percentage of renters spending 30% of their income on rent” all had negative associations with overdose deaths. As the first three factors are all indicators of community stability, their negative relationships with overdose risk are not surprising. By contrast, the negative association between the “percentage of renters spending 30% of their income on rent” (an indicator of financial hardship) and OODs was unexpected and will require further investigation.

Two factors displayed regional spatial variation (Fig. 1). “Crime rate” influenced OODs in central, southern, and western parts of Milwaukee County but had little influence in northeastern Milwaukee County. “Inequality of household income” was associated with OODs in central and southern Milwaukee County, but not in other regions. Such variations are better understood with finer grained analysis, when examined in the context of inequalities shaped by racial segregation. To examine these spatial variabilities, we conducted regional MGWR modeling at the census tract levels (Table 4).

Fig. 1
figure 1

Composite maps of MGWR parameter estimate surfaces for crime rate and inequality of household income which tend to show regional patterns of spatial heterogeneity

Surprisingly, “naloxone availability” parameter estimates displayed significant local spatial variation (Fig. 2). Timely administration of naloxone can prevent OODs, and naloxone has been used by emergency medical personnel to pharmacologically reverse overdoses for decades [61]. Since 1996, community-based programs have provided naloxone to opioid users, their communities, and service/care providers to reverse potentially fatal opioid overdoses [62]. As expected, “naloxone availability” had a strong negative association with the OOD rate in affluent, predominantly White, and suburban census tracts.

Fig. 2
figure 2

Composite map of MGWR parameter estimate surfaces for naloxone availability which tend to show local patterns of spatial heterogeneity

In contrast, census tracts with high population density or tracts with predominantly African American or Hispanic population displayed different results. Here, “naloxone availability” parameter estimates are either not significantly different from zero or, in Hispanic neighborhoods, have a positive association with OODs. In other words, naloxone availability is not having a negative impact on OODs in predominantly African American and Hispanic communities. The reasons for the observed spatial variations are unclear and understanding them will be critical for guiding neighborhood-level policies and practices. Potential contributing factors include the types of opioids or drug combinations used, hesitancy to engage agencies due to fear of incarceration, lack of adequate educational programs, drug use in isolation versus with peers, the absence of effective overdose follow-up programs, and barriers to naloxone access (versus availability).

Racial segregation and its associated inequalities have such profound effects on neighborhood wellbeing in Milwaukee County that it is perceivable that residents are living in three different worlds. For example, the overall Milwaukee County crime rate in 2019 was 42.9 per 1000 residents. In predominantly White communities, it was 18.1 per 1000 residents. However, the crime rate in predominantly African American communities was 83.9 per 1000 residents, while in predominantly Hispanic communities, it was 52.4 per 1000 residents. Thus, in a hyper-segregated metropolitan area such as Milwaukee, analyzing the opioid crisis with respect to segregation and structural inequalities is vital to unpacking the different layers of the epidemic. Accordingly, many associations with overdose deaths appeared to vary with the racial/ethnic composition of communities. To explicate these relationships, we divided Milwaukee County into three regions: census tracts in which White (non-Hispanic) residents are the majority, tracts in which African American residents are the majority, and tracts in which Hispanic residents are the majority.

Three independent Pearson correlation coefficient and VIF evaluations were conducted to detect multicollinearity in each region. To select a subset of variables by removing non-significant explanatory variables, stepwise forward regression analyses were conducted. Different explanatory variables were selected for each region to develop a robust model to investigate variability of the dependent variable “opioid drug overdose death rate” in three different regions of Milwaukee County. In Table 3, the impactful variables identified by stepwise regression analyses for each region are presented.

Table 3 Stepwise regression selected explanatory variables included in the modeling of opioid overdose death in the study for each region

By dividing Milwaukee County into three regions according to segregated neighborhoods and conducting separate models, we observed significant improvement in the accuracy of geospatial modeling, based on comparisons of R2 values and AICs [55]. The R2 for the predominantly Hispanic community increased from 0.969 in the global MGWR model to 0.995 in the regional MGWR model, and the AIC decreased to −28 in the regional scale MGWR model from 2 in the global MGWR model. The R2 for the predominantly Black community increased to 0.718 in the MGWR model from 0.437 in the global model, and the AIC decreased to 219 in the regional MGWR model from 252 in the global model. Similarly, the R2 for the predominantly White community increased to 0.746 in the MGWR model from 0.480 in the global model and the AIC decreased to 335 in the MGWR model from 397 in the global model.

The results in Table 4 show variations across census tracts based on their racial/ethnic compositions. An examination of spatial variance in the parameter estimations of regional MGWR supports the possibility of differential influence of interventions/policies across communities defined by racial/ethnic compositions. These findings are addressed in the “Discussions.”

Table 4 MGWR results: Regional determinant of opioid overdose death among different races/ethnicities

Discussions

As seen in Table 4, the relationships of some factors with overdose deaths were similar across tracts. For example, communities with higher crime rates or those comprised predominantly of middle-aged populations had higher rates of OODs, regardless of racial/ethnic composition. In other cases, relationships varied across census tracts. Positive associations between “inequality of household income” and overdose deaths were evident in the predominantly African American and White communities, but not in the Hispanic community. Prescription opioid availability (historically dispensed oxycodone) and the prevalence of disability had strong positive associations in both Hispanic and White communities, but not within the African American community. In contrast, the percentage of families with young children had strong negative associations with overdose death rates in the Hispanic and White communities, but not in the African American community. Living in a stable neighborhood with more access to healthcare was associated with lower OOD rates in the African American and Hispanic communities, but not in the White communities. Marital status (divorce or death of a spouse), population density, and alcohol availability were positively associated with OODs in the African American majority, but not among the Hispanic or White communities. College enrollment was negatively associated with OODs, and a history of physical pain was positively associated with OODs in the Hispanic American community but not in the predominantly African Americans or White communities. Healthcare service accessibility influenced OOD rates in the predominantly White community, but not in the African American or Hispanic majority communities. The precise reasons for geographic variation in these relationships are unclear and require further investigation.

Perhaps the most surprising observation was that naloxone availability negatively influenced OODs in the White majority, but not the African American or Hispanic majority communities. In fact, in areas of Milwaukee that are predominantly Hispanic, there was a positive association between naloxone availability and OODs, potentially reflecting community efforts to increase distribution, despite ineffectiveness in reducing overdose risk. The NARCAN® Direct Program, supported by the Wisconsin Department of Health Services, is an integral component of the state’s response to the opioid crisis and provides Narcan at no cost to community agencies to distribute with the goal of reducing OODs. While it is possible that community members are acquiring naloxone from other sources, these findings suggest that, despite naloxone availability, there are obstacles that are preventing its effective use in African American and Hispanic communities. These may include the willingness or ability of community members to engage the agencies through which naloxone is being provided and/or lack of awareness of the availability of naloxone and the benefits of its use. Alternatively, other factors, such as the types/potency of opioids or drug combinations involving non-opioids being used, could vary across communities.

Another surprising observation was related to the association between incarceration rates and OODs. Unexpectedly, incarceration rates influenced overdose deaths in the White majority community but not within the African American or Hispanic communities. African Americans are no more likely than White Americans to use illicit drugs but are 6–10 times more likely to be incarcerated for drug offenses [63]. While the threat of incarceration may be a deterrent to drug use and the risk for OODs is negligible during incarceration, overall overdose rates are higher in former inmates [64]. Our findings raise questions about whether or not more aggressive policing and criminalization of drug possession in communities of color are effective measures for preventing overdose.

OODs are declining in White Americans but are on the rise in communities of color, particularly among African Americans [33]. This trend is also evident in Milwaukee County [31]. It raises concerns that the policies, resources, and interventions that are being designed and implemented are disproportionately benefiting White communities, thereby creating inequity and widening health disparities. Understanding how factors differentially influence OODs across communities should guide targeted policies and approaches aimed at prevention, intervention, and harm reduction. While some of the observations may specifically apply to communities in Milwaukee, many of the relationships are likely evident in communities in other metropolitan areas.

There are a number of research gaps and opportunities for future work that build on the present findings. First, our data collection ended just prior to the onset of the COVID-19 pandemic. Evidence suggests that overdose deaths have increased during the pandemic [65, 66], and the impact of the pandemic on the overdose crisis likely also varies across communities. Second, while the list of analyzed factors was extensive, many factors of likely importance (e.g., types of opioids or other drugs that contributed to relapse, types of community programs and treatment approaches, mental health indicators) were not included. Third, we georeferenced overdoses based on the location of the overdose but not the victims’ addresses of residence. Finally, there is a need to extend our analyses to regions beyond Milwaukee County so that we can identify and address factors that are differentially influential in suburban and rural communities. In addition to expanding our MGWR analysis, future work will also involve community engagement and qualitative assessments so that we can further understand these relationships prior to sharing data with policymakers and community organizations to support the targeted implementation of policies and initiatives.

Conclusion

In conclusion, we report that the factors that influence OODs in Milwaukee County vary across the diverse communities in this highly segregated metropolitan area. The observed geographic variation in relationships includes the impact of naloxone availability and incarceration rates on overdose deaths with pronounced differences between White communities and communities of color. Understanding community-level factors that contribute to overdose risk should guide targeted community-level solutions. Overall, our findings demonstrate the value of precision epidemiology using MGWR analysis for defining and guiding responses to public health challenges.