Predicting Heroin Overdose in Mesa
Sylvia Lam and Claudia Low
13 December, 2024
Introduction
This report aims to predict heroin overdose risk in Mesa, Arizona. Heroin is an illegal, dangerous and addictive opioid drug that has pervaded society. Heroin abuse is a leading cause of preventable deaths in the United States, necessitating the implementation of actionable measures to address this issue. The City of Mesa ranks as the second-highest city for drug use and is among the top cities for heroin usage nationwide (American Addiction Centers, 2024). Addressing this crisis requires targeted interventions to allocate limited public health resources effectively.
To address the heroin overdose issue, we propose NARVEND, a geospatial dashboard designed to support decision-making for public health officials in Mesa to identify areas most in need of Narcan nasal spray vending machines. Narcan nasal spray can reverse the life-threatening effects of opioid overdose, but timing is critical, and it must be administered as soon as overdose symptoms are observed. This is because heroin is a short-acting opioid that can result in life-threatening respiratory depression, leading to unconsciousness, brain damage, or death quickly. Placing these vending machines in locations that have high risk of heroin overdose but lack nearby Narcan supply will expand access to Narcan for both drug users and non-drug users in the vicinity to carry around with them or easily access in emergencies. Given budget and resource constraints, we recommend deploying seven vending machines to dispense free Narcan nasal spray, the the life-saving medication.
In order to determine the optimal placement of these vending machines, we developed a geospatial risk prediction model trained on reported heroin overdoses in Mesa in 2020, incorporating various environmental and socio-economic factors. The predictions of the model form the risk basemap within the NARVEND dashboard. By overlaying current resource locations, the dashboard highlights underserved areas with the highest overdose risks but limited access to existing resources. These identified gaps inform data-driven recommendations for the most suitable placement of life-saving Narvan nasal spray vending machines, ensuring that resources are allocated where they are needed most.
The presentation pitch to the Mesa public health officials on this project can be viewed here.
This analysis can be replicated across other cities facing opioid overdose issues, and provides critical value in identifying underserved communities most in need of accessible Narcan supply.
# Setup and load in required packages
knitr::opts_chunk$set(echo = TRUE, warning=FALSE)
library(tidyverse)
library(sf)
library(RSocrata)
library(viridis)
library(spatstat.explore)
library(raster)
library(spdep)
library(FNN)
library(grid)
library(gridExtra)
library(knitr)
library(kableExtra)
library(tidycensus)
library(classInt)
library(tmap)
library(dplyr)
library(osmdata)
options(scipen = 999)
root.dir = "https://raw.githubusercontent.com/urbanSpatial/Public-Policy-Analytics-Landing/master/DATA/"
source("https://raw.githubusercontent.com/urbanSpatial/Public-Policy-Analytics-Landing/master/functions.r")# Load in opioid data
opioid_mesa <- read.csv("Fire_and_Medical_Opioid_Overdose_Incidents_20241121.csv") %>%
filter(!is.na(Latitude)) %>%
st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326) %>%
st_transform("ESRI:102249")
# Load in spatial data
## Maricopa Boundary
maricopaBoundary <- st_read("maricopaCounty.geojson") %>%
filter(NAME == "MARICOPA") %>%
st_transform("ESRI:102249")
## Mesa Boundary
mesaBoundary <- st_read("mesaBoundary.geojson") %>%
st_transform("ESRI:102249")
## Airport Boundary
airportBoundary <- st_read("airport.geojson") %>%
filter(OBJECTID == 6711) %>%
st_transform("ESRI:102249")
# Load in census data
census_data <-
get_acs(geography = "tract", variables = c("B01003_001E","B19013_001E", "B02001_002E", "B01001_002E",
"B17001_002E", "B15002_015E", "B15002_032E", "B15011_001E", "B23025_001E",
"B23025_002E"),
year=2022, state="AZ", county = "Maricopa County", geometry=TRUE, output="wide") %>%
st_transform('ESRI:102249') %>%
rename(TotalPop = B01003_001E,
MedHHInc = B19013_001E,
Whites = B02001_002E,
Male = B01001_002E,
BelowPovertyLine = B17001_002E,
MaleBach = B15002_015E,
FemaleBach = B15002_032E,
TotalBach = B15011_001E,
LaborForceTotal = B23025_001E,
InLaborForce = B23025_002E,
) %>%
mutate(pctWhite = (Whites / TotalPop),
pctMale = (Male / TotalPop),
pctBelowPovertyLine = (BelowPovertyLine/TotalPop),
pctBach = ((MaleBach + FemaleBach) / TotalBach),
pctEmployed = (InLaborForce / LaborForceTotal)
) %>%
dplyr::select(-NAME, -starts_with("B"))
# Load in risk factor variables
## Parks [https://data.mesaaz.gov/Parks-Recreation-and-Community-Facilities/Parks-Locations-And-Amenities/djym-pkpp/about_data]
parks <- read.csv("parks.csv") %>%
dplyr::select(Latitude, Longitude) %>%
st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326) %>%
st_transform("ESRI:102249") %>%
mutate(Legend = "Parks")
## Homeless [https://citydata.mesaaz.gov/Community-Services/Unsheltered-Point-In-Time-PIT-Count-2022-Details-M/efjd-c5mi/about_data]
homeless <- read.csv("homeless.csv") %>%
dplyr::select(Latitude, Longitude) %>%
na.omit %>%
st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326) %>%
st_transform("ESRI:102249") %>%
mutate(Legend = "Homeless")
## Crime [https://data.mesaaz.gov/Police/Police-Incidents/39rt-2rfj/about_data]
crime_incidents <- read.csv("crime_incidents.csv") %>%
na.omit %>%
st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326) %>%
st_transform("ESRI:102249")
### All crime incidents
all_incidents <- read.csv("crime_incidents.csv") %>%
dplyr::select(Latitude, Longitude, Crime.Type) %>%
na.omit %>%
st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326) %>%
st_transform("ESRI:102249") %>%
mutate(Legend = "Crime_Incidents") %>%
dplyr::select(-Crime.Type)
### Assault incidents
assault <- crime_incidents %>%
dplyr::filter(grepl("ASSAULT", Crime.Type, ignore.case = TRUE)) %>%
dplyr::select(geometry) %>%
mutate(Legend = "Assault")
### Burglary incidents
burglary <- crime_incidents %>%
dplyr::filter(grepl("BURGLARY", Crime.Type, ignore.case = TRUE)) %>%
dplyr::select(geometry) %>%
mutate(Legend = "Burglary")
## Drug incidents
drug <- crime_incidents %>%
dplyr::filter(grepl("DRUG", Crime.Type, ignore.case = TRUE)) %>%
dplyr::select(geometry) %>%
mutate(Legend = "Drug")Exploratory Analysis
Temporal and Spatial Distribution of Heroin Overdose
The outcome variable in the geospatial risk model is counts of reported heroin overdose, provided by the Fire and Medical Emergency Department of Mesa. The number of opioid overdose cases from 2018 to 2023 is presented in Figure 1. There was a surge in overdose cases in Mesa in 2020 and 2021, likely attributed to the increased social isolation and economic stress from the Covid-19 pandemic which increased drug usage and risk of overdose. The 2020 overdose case data highlights a period of crisis, revealing the most vulnerable populations, and was therefore utilized in the model. This geospatial risk model aims to map the latent risk of heroin overdose, the underlying likelihood of overdose across space. This is unlikely to be the same as the reported overdose density but rather, is inferred from chosen risk factors and observed spatial patterns.
The spatial distribution of heroin overdose points in 2020 are first plotted across the map of Mesa in Figure 2. It is immediately noticeable and important to highlight that these points are not the exact locations of the overdose incidents, as the dataset has snapped the points to certain standard locations by approximately 1/3 mile each in order to protect anonymity. One notices a concentration of overdose incidents in the west central area.
# Bar chart of number of opioid overdose cases by year
opioid_mesa %>%
st_drop_geometry() %>%
na.omit() %>%
dplyr::select(c(Year, Incident_ID, PD.Save)) %>%
count(Year) %>%
ggplot() +
geom_bar(aes(x = Year, y = n, fill = as.factor(Year)),
position = "dodge", stat = "identity", fill = "#8a9ab1") +
scale_x_continuous(breaks = seq(min(opioid_mesa$Year), max(opioid_mesa$Year), by = 1)) +
labs(
title = "Opioid Overdose Cases by Year",
y = "Number of Opioid Overdose Cases",
x = "Year",
caption = "Figure 1"
) +
theme_minimal() +
theme(
axis.text.x = element_text(angle = 45, vjust = 0.5),
panel.grid = element_blank(),
plot.caption = element_text(hjust = 0)
)
# Remove the airport from Mesa boundary
mesaBoundaryNoAirport <- st_difference(mesaBoundary, airportBoundary)
# Filter opioid cases to only 2020 within Mesa City boundary
opioid_mesa2020 <- opioid_mesa %>%
filter(Year == 2020)
opioid_mesa2020 <- opioid_mesa2020 %>%
filter(sapply(st_intersects(opioid_mesa2020, mesaBoundaryNoAirport), length) > 0)
# Plot out spatial distribution of overdose points in 2020 across Mesa
ggplot() +
geom_sf(data = mesaBoundaryNoAirport, fill = "#333333", colour = "black", lwd = 0) +
geom_sf(data = opioid_mesa2020) +
geom_sf(data = opioid_mesa2020, colour="white", size=0.5, alpha = 0.7, show.legend = "point") +
labs(title= "Opioid Overdose cases in Mesa, 2020",
caption = "Figure 2") +
theme_minimal() +
theme(panel.grid = element_blank(),
axis.text = element_blank(),
axis.title = element_blank(),
axis.ticks = element_blank(),
plot.caption = element_text(hjust = 0))
Creating fishnet
A fishnet map overlaying Mesa was then created to observe the density distribution using counts of heroin overdose incidents per grid cell.
Hexagonal rather than square grids are used as the centroids of the grid are used in some parts of feature engineering. Hexagonal grids ensure equal distances from a cell’s centroid to its neighbors, unlike square grids where diagonal neighbors are farther away, leading to uneven distances. Additionally, hexagonal grids minimize edge effects by enabling smoother transitions at the boundaries, reducing distortions for cells near the edges.
Due to the lack of a precise location in the dataset of heroin overdose, the fishnet grid overlaying the city has a cell size of 1000 meters, which is approximately 2/3 of a mile. This cell size was chosen as considering the snapping of points to the nearest 1/3 mile standardized locations, a grid cell that is 2/3 mile wide will ensure the genuine location of the overdose is still within the same grid cell. Although larger cell size reduces granularity in risk predictions, it still ensures that overdose events are captured within a manageable spatial resolution while minimizing location misclassification errors, allowing for meaningful risk predictions within such data limitations.
It should be noted that the Phoenix-Mesa Gateway Airport was removed from the map as part of the analyses as the census variable values for this tract were extreme outliers. There are no residents in the airport and security in the airport is also much tighter, thus heroin overdose cases are highly unlikely.
# Create a fishnet grid with a cell size of 1000m
fishnet <-
st_make_grid(maricopaBoundary,
cellsize = 1000,
square = FALSE) %>%
.[maricopaBoundary] %>%
st_sf() %>%
mutate(uniqueID = 1:n())
# Crop the fishnet to Mesa boundary
mesaFishnet_old <- st_intersection(fishnet, mesaBoundaryNoAirport)
mesaFishnet <- mesaFishnet_old %>%
mutate(uniqueID = 1:n())
# Aggregate overdose points to the fishnet
## Add a value of 1 to each overdose case, sum them with aggregate
overdose_net <-
dplyr::select(opioid_mesa2020) %>%
mutate(countOverdoseCase = 1) %>%
aggregate(., mesaFishnet, sum) %>%
mutate(countOverdoseCase = replace_na(countOverdoseCase, 0),
uniqueID = 1:n(),
cvID = sample(round(nrow(mesaFishnet) / 24),
size=nrow(mesaFishnet), replace = TRUE))
# Map number of overdose points per fishnet cell
ggplot() +
geom_sf(data = overdose_net, aes(fill = countOverdoseCase), color = NA) +
scale_fill_viridis(name = "Count of Overdose Cases") +
labs(title = "Count of Overdose Case for the fishnet",
caption = "Figure 3") +
theme_minimal() +
theme(panel.grid = element_blank(),
axis.text = element_blank(),
axis.title = element_blank(),
axis.ticks = element_blank(),
plot.caption = element_text(hjust = 0))
Risk factors
A range of risk factors were selected to predict for heroin overdose in the risk prediction model.
Crime data was taken from the Mesa Police Department. All crimes were selected as one predictor variable, as well as specific crimes of assault, burglary and drug-related crime. Drug-related crime include arrests for possessing drug-related paraphernalia and drug use. Heroin overdose would likely also correlate with areas of violent and property crime such as assault and burglary. Drug use in general likely highly correlate with heroin overdose as these are locations where there may be more frequent substance usage.
Other environmental factors considered include the presence of parks, measured by the number of parks within each grid cell. Public parks are often associated with increased drug-related activity due to their accessibility, and limited surveillance. Spatial data on parks was provided by the Parks, Recreation and Community Facilities Department of Mesa. Proximity to homeless individuals was also considered a risk factor due to the association between homelessness and increased drug use. Homeless hotspots often coincide with areas of higher substance use due to limited access to stable housing, healthcare, and support services, thus this risk factor helps account for spatial concentrations of vulnerable populations, which may be useful in the prediction of the spatial distribution of heroin overdoses. The dataset from the Mesa data hub is a point-in-time count of unsheltered individuals on a particular morning in January 2022, but considering the difficulty of having a true spatial distribution of homeless individuals in a city, this dataset is the closest proxy.
Census data from the American Community Survey 2022 5-year estimates were next chosen to account for socio-demographic variables of the census tract. Substance abuse is generally correlated with areas with lower income, employment levels and education levels. The variables included were percent white, percent male, percent below poverty line, percent with Bachelor’s degree, percent employed and median household income. The area of each census tract within each grid cell was calculated, with the corresponding census variable values weighted based on the proportion of the tract’s area intersecting the hexagonal cell. These weighted values were then aggregated for each fishnet grid cell.
# Calculate the area of each census tract within a grid cell and find the proportion of the grid cell it covers
mesaFishnet <- st_transform(mesaFishnet, crs = st_crs(census_data))
census_data$area_census_tract <- st_area(census_data)
fishnet_intersection <- st_intersection(mesaFishnet, census_data)
fishnet_intersection$area_intersection <- st_area(fishnet_intersection)
fishnet_intersection <- fishnet_intersection %>%
mutate(area_proportion = as.numeric(area_intersection) / as.numeric(area_census_tract))
# Weight the census variable by the area proportion
fishnet_intersection <- fishnet_intersection %>%
mutate(weighted_pctWhite = pctWhite * area_proportion,
weighted_pctMale = pctMale * area_proportion,
weighted_pctBelowPovertyLine = pctBelowPovertyLine * area_proportion,
weighted_pctBach = pctBach * area_proportion,
weighted_pctEmployed = pctEmployed * area_proportion,
weighted_MedHHInc = MedHHInc * area_proportion)
# Aggregate weighted variables for each fishnet grid cell
fishnet_census <- fishnet_intersection %>%
group_by(uniqueID) %>%
summarize(
pctWhite = sum(weighted_pctWhite, na.rm = TRUE) / sum(area_proportion, na.rm = TRUE),
pctMale = sum(weighted_pctMale, na.rm = TRUE) / sum(area_proportion, na.rm = TRUE),
pctBelowPovertyLine = sum(weighted_pctBelowPovertyLine, na.rm = TRUE) / sum(area_proportion, na.rm = TRUE),
pctBach = sum(weighted_pctBach, na.rm = TRUE) / sum(area_proportion, na.rm = TRUE),
pctEmployed = sum(weighted_pctEmployed, na.rm = TRUE) / sum(area_proportion, na.rm = TRUE),
MedHHInc = sum(weighted_MedHHInc, na.rm = TRUE) / sum(area_proportion, na.rm = TRUE)
)The risk factors were plotted as part of an exploratory analysis. One can observe how the different types of crimes all concentrate in a similar western area of the city as the overdose incidents. Unsheltered individuals also concentrated in that same area. There is also an overall higher percentage of individuals below poverty and with lower median household income in west Mesa.
# Create a dataset that combines all the risk factor variables and sync it to the fishnet
vars_net_neigh <- parks %>%
rbind(parks, homeless, all_incidents, assault, burglary, drug) %>%
st_join(mesaFishnet, join=st_within) %>%
st_drop_geometry() %>%
group_by(uniqueID, Legend) %>%
summarize(count = n()) %>%
left_join(mesaFishnet, ., by = "uniqueID") %>%
spread(Legend, count, fill=0) %>%
dplyr::select(-`<NA>`) %>%
ungroup()
censusonly <- st_drop_geometry(fishnet_census)
vars_net <- left_join(vars_net_neigh, censusonly, by = "uniqueID") %>%
dplyr::select(-CITY, -OBJECTID.1, -Zoning, -Description, -Acres, -Previous_Zoning, -RelatedCase, -Ordinance, -Link, -OBJECTID)
vars_net.long <-
gather(vars_net, Variable, value, -geometry, -uniqueID)
vars_net.long <- vars_net.long %>%
mutate(Variable = recode(Variable,
"pctWhite" = "Percent White",
"pctMale" = "Percent Male",
"pctBelowPovertyLine" = "Percent Below Poverty Line",
"pctBach" = "Percent with Bachelor's Degree",
"pctEmployed" = "Percent Employed",
"MedHHInc" = "Median Household Income"))
# Plot risk factor values by fishnet
vars <- unique(vars_net.long$Variable)
mapList <- list()
for(i in vars){
mapList[[i]] <-
ggplot() +
geom_sf(data = filter(vars_net.long, Variable == i), aes(fill=value), colour=NA) +
scale_fill_viridis(name = "", option = "rocket") +
labs(title=i) +
theme_void() +
theme(plot.title = element_text(hjust = 0.5))
}
mapList <- lapply(mapList, ggplotGrob)
empty_space <- textGrob("", gp = gpar(col = "transparent"))
grid.layout <- grid.arrange(
empty_space,
arrangeGrob(grobs = mapList, ncol = 3),
top = textGrob("Risk Factors by Fishnet",
gp = gpar(fontsize = 18)),
bottom = textGrob("Figure 4", x = 0, hjust = 0, gp = gpar(fontsize = 10)),
heights = unit(c(0.05,0.8), "npc")
)
final_net <-
left_join(overdose_net, st_drop_geometry(vars_net), by="uniqueID")
# Join in centroids of fishnets to polygon for census tract neighborhoods
final_net <-
st_centroid(final_net) %>%
st_join(dplyr::select(census_data, GEOID), by = "uniqueID") %>%
st_drop_geometry() %>%
left_join(dplyr::select(final_net, geometry, uniqueID)) %>%
st_sf() %>%
na.omit()Correlation Analysis
A correlation plot for each risk factor and the outcome variable of heroin overdose can be seen in Figure 5. It is observed that assault incidents, burglary incidents, and all crime incident counts had a strong correlation to the outcome variable ranging from 0.72 to 0.77.
Median household income and percent below poverty had a moderate correlation with overdose at -0.36 and 0.44 respectively, but the other census variables capturing gender, education and employment rate had the lowest correlations of 0.1 or less, on the other hand were found to have little to no correlation with the outcome variables. Census variables may have had lower correlations as other more important environmental factors and income factors may overshadow simple demographic indicators. The method used to calculate census variable values per grid cell also assumed uniform distribution of a census variable within each census tract, and a singular census variable for a 2/3 mile by 2/3 mile grid cell may also mask more local variations, resulting in weaker correlations.
# Plot correlation of risk factor variables to overdose count
correlation.long <-
st_drop_geometry(final_net) %>%
dplyr::select(-uniqueID, -cvID, -GEOID) %>%
gather(Variable, Value, -countOverdoseCase)
correlation.long <- correlation.long %>%
mutate(Variable = recode(Variable,
"pctWhite" = "Percent White",
"pctMale" = "Percent Male",
"pctBelowPovertyLine" = "Percent Below Poverty Line",
"pctBach" = "Percent with Bachelor's Degree",
"pctEmployed" = "Percent Employed",
"MedHHInc" = "Median Household Income"))
correlation.cor <-
correlation.long %>%
group_by(Variable) %>%
summarize(correlation = cor(Value, countOverdoseCase, use = "complete.obs"))
ggplot(correlation.long, aes(Value, countOverdoseCase)) +
geom_point(size = 0.2) +
geom_text(data = correlation.cor, aes(label = paste("r =", round(correlation, 2))),
x=-Inf, y=Inf, vjust = 1.5, hjust = -.1) +
geom_smooth(method = "lm", se = FALSE, colour = "black") +
facet_wrap(~Variable, ncol = 2, scales = "free") +
labs(title = "Heroin Overdose count as a function of risk factors",
caption = "Figure 5",
x = "Correlation value",
y = "Count of Overdose Cases") +
theme_minimal()+
theme(plot.caption = element_text(hjust = 0))
Spatial Autocorrelation
The presence of spatial autocorrelation patterns of heroin overdose was then investigated, which is the extent to which overdose is correlated with itself across space. Local Moran’s I values are high in the west, indicating similar values clustering together. The low p-values indicate statistical significance. In reference back to Figure 3, this spatial cluster is a hotspot of high overdose count values. A large significant hotspot can be seen in the central west, along with a small one in the center
To capture this spatial process in the prediction model, two additional risk factors were included: whether a grid cell was part of a significant hotspot and the average nearest neighbor distance from each cell centroid to its nearest significant cluster. Incorporating these local spatial characteristics of heroin overdose can improve the accuracy of risk predictions.
# Define neighbors as queen contiguity
final_net.nb <- poly2nb(as_Spatial(final_net), queen=TRUE)
final_net.weights <- nb2listw(final_net.nb, style="W", zero.policy=TRUE)
# Calculate Local Moran's I and join results to fishnet
final_net.localMorans <-
cbind(
as.data.frame(localmoran(final_net$countOverdoseCase, final_net.weights)),
as.data.frame(final_net)) %>%
st_sf() %>%
dplyr::select(Overdose_Count = countOverdoseCase,
Local_Morans_I = Ii,
P_Value = `Pr(z != E(Ii))`) %>%
mutate(Significant_Hotspots = ifelse(P_Value <= 0.05, 1, 0)) %>%
gather(Variable, Value, -geometry)
#Plot out Overdose Count, Local Moran's I, P-Value and Significant Hotspots
title_mapping <- c(
"Overdose_Count" = "Overdose Count",
"Local_Morans_I" = "Local Moran's I",
"P_Value" = "P Value",
"Significant_Hotspots" = "Significant Hotspots"
)
vars <- unique(final_net.localMorans$Variable)
varList <- list()
for(i in vars){
varList[[i]] <-
ggplot() +
geom_sf(data = mesaBoundaryNoAirport, fill = "white") +
geom_sf(data = filter(final_net.localMorans, Variable == i),
aes(fill = Value), colour = NA) +
scale_fill_viridis(name = "") +
theme_minimal() +
labs(title = title_mapping[[i]]) +
theme(legend.position = "bottom",
panel.grid = element_blank(),
axis.text = element_blank(),
axis.title = element_blank(),
axis.ticks = element_blank())
}
plot_layout_moran <- arrangeGrob(
grobs = varList,
ncol = 4,
top = textGrob("Local Moran's I Statistics, Heroin Overdose", gp = gpar(fontsize = 14)),
bottom = textGrob("Figure 6", x = 0, hjust = 0, gp = gpar(fontsize = 8))
)
grid.newpage()
grid.draw(plot_layout_moran)
# Add spatial process variables to fishnet
final_net <-
final_net %>%
mutate(overdose.isSig =
ifelse(localmoran(final_net$countOverdoseCase,
final_net.weights)[,5] <= 0.0000001, 1, 0)) %>%
mutate(overdose.isSig.dist =
nn_function(st_coordinates(st_centroid(final_net)),
st_coordinates(st_centroid(
filter(final_net, overdose.isSig == 1))), 1))Results and Model Evaluation
Modeling and Cross Validation
The model was run with the 12 risk factors, with Leave-One-Group-Out Cross Validation by census tract then conducted. This test evaluates how well the model performs not only at a broader city-wide scale, but at more local spatial scales. Each census tract serves as a validation set once, while the remaining neighborhoods form the training set. This process is repeated, with each neighborhood taking a turn as the holdout fold in successive iterations.
# Collate all risk factors for regression
reg.ss.vars <- c("Assault", "Burglary",
"Homeless", "Parks", "pctWhite", "pctMale",
"pctBelowPovertyLine", "pctBach", "pctEmployed", "MedHHInc", "overdose.isSig", "overdose.isSig.dist")
# Leave-One-Group out Cross Validation
reg.ss.spatialCV <- crossValidate(
dataset = final_net,
id = "GEOID",
dependentVariable = "countOverdoseCase",
indVariables = reg.ss.vars) %>%
dplyr::select(cvID = GEOID, countOverdoseCase, Prediction, geometry)
# Regression summary
reg.summary <-
mutate(reg.ss.spatialCV, Error = Prediction - countOverdoseCase,
Regression = "Spatial LOGO-CV: Spatial Process") %>%
st_sf() Calculating Errors across space
Across census tracts, most have a relatively small mean absolute error (MAE), although there are several outliers with high error. The census tract with an ID of 04013421302 had the highest mean error, and it lies within the spatial cluster in the west, along with most other tracts which have mean absolute error larger than eight.
# Calculate errors by census tract
error_by_reg_and_fold <-
reg.summary %>%
group_by(Regression, cvID) %>%
summarize(Mean_Error = mean(Prediction - countOverdoseCase, na.rm = T),
MAE = mean(abs(Mean_Error), na.rm = T),
SD_MAE = mean(abs(Mean_Error), na.rm = T)) %>%
ungroup()
# Plot histogram of OOF (out of fold) errors
error_by_reg_and_fold %>%
ggplot(aes(MAE)) +
geom_histogram(bins = 30, colour="black", fill = "#8a9ab1") +
facet_wrap(~Regression) +
geom_vline(xintercept = 0) + scale_x_continuous(breaks = seq(0, 8, by = 1)) +
labs(title="Distribution of MAE", subtitle = "LOGO cross validation",
x="Mean Absolute Error", y="Count",
caption = "Figure 7") +
theme_minimal() +
theme(plot.caption = element_text(hjust = 0))
# Explore errors of specific neighborhoods
check <- error_by_reg_and_fold %>%
arrange(desc(MAE))
neighborhoods_explore<- final_net %>%
group_by(GEOID) %>%
summarize(total = sum(countOverdoseCase, na.rm = T))Overall, the mean MAE is 2.09 with a standard deviation of 3.4. While this indicates a moderate level of prediction accuracy, it is important to emphasize that the model predicts latent risk of heroin overdose rather than directly estimating observed overdose counts. This distinction means the model identifies areas that may be at risk of overdose even if no overdoses have been officially reported, reflecting the underlying spatial patterns associated with overdose risk factors.
Ultimately, evaluating model accuracy is challenging due to potential selection bias in the reported overdose data used for training. Under-reporting is an issue as heroin overdoses may go unreported due to fear of legal consequences, social stigma, or the absence of bystanders to call emergency services. It is also possible that overdoses are disproportionately reported in high-surveillance areas, even if actual risk levels are comparable elsewhere. Given that the model is trained on such biased observational data, there is a possibility that these biases are embedded in the predictions.
Nonetheless, our model for predicting latent risk remains valuable in determining optimal vending machine locations, offering a more informed approach compared to random or arbitrary decision-making.
# Present mean MAE and standard deviation MAE of regression
st_drop_geometry(error_by_reg_and_fold) %>%
group_by(Regression) %>%
summarize(Mean_MAE = round(mean(MAE), 2),
SD_MAE = round(sd(MAE), 2)) %>%
kable(caption = "MAE of Regression Model") %>%
kable_styling("striped", full_width = F) %>%
footnote(
general = "Table 1",
general_title = "",
footnote_as_chunk = TRUE,
threeparttable = TRUE
)| Regression | Mean_MAE | SD_MAE |
|---|---|---|
| Spatial LOGO-CV: Spatial Process | 2.09 | 3.4 |
| Table 1 |
# Plot errors by fishnet
error_by_reg_and_fold %>%
filter(str_detect(Regression, "LOGO")) %>%
ggplot() +
geom_sf(aes(fill = MAE)) +
facet_wrap(~Regression) +
scale_fill_viridis() +
labs(title = "Overdose errors by LOGO-CV Regression",
caption = "Figure 8") +
theme_minimal() +
theme(panel.grid = element_blank(),
axis.text = element_blank(),
axis.title = element_blank(),
axis.ticks = element_blank(),
plot.caption = element_text(hjust = 0))
Model performance across census groups
The generalizability of the model across census groups was next evaluated. It has to be acknowledged that our model is not particularly generalizable as it tends to predict less accurately for majority non-white and lower-income neighborhoods. These areas exhibit a higher mean error in the regression model. It is ultimately challenging for the model to be highly generalizable due to the inherent nature of the spatial distribution of heroin overdose, where there is mainly one large hotspot. This concentration in overdose cases creates a disparity and limits its ability to accurately predict risk in areas with different socio-demographic profiles. As the majority of overdose incidents are concentrated within this single hotspot, predicting risk in areas outside this zone becomes difficult due to the lack of representative data.
While the generalizability of the model is constrained by these circumstance, it is still a valuable tool for identifying high-risk areas and informing decision-making. The model provides a useful foundation for allocating Narcan vending machines more effectively to locations with the highest overdose risks. Future improvements over time, such as incorporating additional risk factors and continuous refinement of the model, could enhance its applicability to a broader range of neighborhoods.
# Retrieve census data to differentiate census tract neighborhoods by racial and income contexts
census_across_groups <-
get_acs(geography = "tract",
variables = c("B01003_001E","B02001_002E","B19013_001E"),
year=2022, state="AZ", county = "Maricopa County", geometry=TRUE, output = "wide") %>%
st_transform('ESRI:102249') %>%
rename(TotalPop = B01003_001E,
Whites = B02001_002E,
MedHHInc = B19013_001E) %>%
mutate(percentWhite = Whites / TotalPop,
raceContext = ifelse(percentWhite > 0.5, "Majority White", "Majority Non-White"),
incomeContext = ifelse(MedHHInc > mean(MedHHInc, na.rm = TRUE), "High Income", "Low Income"))
# Present mean error by neighborhood racial context
reg.summary %>%
filter(str_detect(Regression, "LOGO")) %>%
st_centroid() %>%
st_join(census_across_groups) %>%
na.omit() %>%
st_drop_geometry() %>%
group_by(Regression, raceContext) %>%
summarize(mean.Error = mean(Error, na.rm = T)) %>%
spread(raceContext, mean.Error) %>%
kable(caption = "Mean Error by neighborhood racial context") %>%
kable_styling("striped", full_width = F) %>%
footnote(
general = "Table 2",
general_title = "",
footnote_as_chunk = TRUE,
threeparttable = TRUE
)| Regression | Majority Non-White | Majority White |
|---|---|---|
| Spatial LOGO-CV: Spatial Process | -1.467314 | 0.1752361 |
| Table 2 |
# Present mean error by neighborhood income context
reg.summary %>%
filter(str_detect(Regression, "LOGO")) %>%
st_centroid() %>%
st_join(census_across_groups) %>%
na.omit() %>%
st_drop_geometry() %>%
group_by(Regression, incomeContext) %>%
summarize(mean.Error = mean(Error, na.rm = T)) %>%
spread(incomeContext, mean.Error) %>%
kable(caption = "Mean Error by neighborhood income context") %>%
kable_styling("striped", full_width = F) %>%
footnote(
general = "Table 3",
general_title = "",
footnote_as_chunk = TRUE,
threeparttable = TRUE
)| Regression | High Income | Low Income |
|---|---|---|
| Spatial LOGO-CV: Spatial Process | 0.0218438 | 0.196284 |
| Table 3 |
Predicting for 2021 Heroin Overdose
The model was then applied to heroin overdose cases in 2021. Figure 9 shows the reported overdose cases in 2021 exhibiting a similar spatial pattern as 2020.
# Filter to opioid cases in only 2021 within Mesa City boundary
opioid_mesa2021 <- opioid_mesa %>%
filter(Year == 2021)
opioid_mesa2021 <- opioid_mesa2021 %>%
filter(sapply(st_intersects(opioid_mesa2021, mesaBoundaryNoAirport), length) > 0)
# Plot out spatial distribution of overdose points in 2021 across Mesa
ggplot() +
geom_sf(data = mesaBoundaryNoAirport, fill = "#333333", colour = "black", lwd = 1) +
geom_sf(data = opioid_mesa2021) +
geom_sf(data = opioid_mesa2021, colour="white", size=0.2, alpha = 0.7, show.legend = "point") +
labs(title= "Opioid Overdose cases in Mesa, 2021",
caption = "Figure 9") +
theme_minimal() +
theme(panel.grid = element_blank(),
axis.text = element_blank(),
axis.title = element_blank(),
axis.ticks = element_blank(),
plot.caption = element_text(hjust = 0))
# Add 2021 cases to fishnet
opioid2021 <- opioid_mesa2021 %>%
dplyr::select(Incident_ID, geometry) %>%
na.omit() %>%
st_transform('ESRI:102249') %>%
.[mesaFishnet,]A density map of these reported cases as well as the risk predictions made by our model were then plotted. While the kernel density and risk predictions both highlight a cluster in the west, the risk predictions tend to overpredict for areas of low density of reported overdose cases. This may potentially be indicative of latent risk where there is risk of overdose despite a lack of officially reported cases.
# Set Kernel Density Estimate (KDE) parameters
overdose_ppp <- as.ppp(st_coordinates(opioid_mesa2020), W = st_bbox(final_net))
overdose_KD.500 <- spatstat.explore::density.ppp(overdose_ppp, 500)
overdose_KDE_sum <- as.data.frame(overdose_KD.500) %>%
st_as_sf(coords = c("x", "y"), crs = st_crs(final_net)) %>%
aggregate(., final_net, mean)
#Create breaks for KDE
kde_breaks <- classIntervals(overdose_KDE_sum$value,
n = 5, "fisher")
# Split into risk categories
overdose_KDE_sf <- overdose_KDE_sum %>%
mutate(label = "Kernel Density",
Risk_Category = classInt::findCols(kde_breaks),
Risk_Category = case_when(
Risk_Category == 5 ~ "5th",
Risk_Category == 4 ~ "4th",
Risk_Category == 3 ~ "3rd",
Risk_Category == 2 ~ "2nd",
Risk_Category == 1 ~ "1st")) %>%
cbind(
aggregate(
dplyr::select(opioid_mesa2021) %>% mutate(overdosecaseCount = 1), ., sum) %>%
mutate(overdosecaseCount = replace_na(overdosecaseCount, 0))) %>%
dplyr::select(label, Risk_Category, overdosecaseCount)#Create breaks for risk predictions
ml_breaks <- classIntervals(reg.ss.spatialCV$Prediction,
n = 5, "fisher")
# Split into risk categories
overdose_risk_sf <-
reg.ss.spatialCV %>%
mutate(label = "Risk Predictions",
Risk_Category =classInt::findCols(ml_breaks),
Risk_Category = case_when(
Risk_Category == 5 ~ "5th",
Risk_Category == 4 ~ "4th",
Risk_Category == 3 ~ "3rd",
Risk_Category == 2 ~ "2nd",
Risk_Category == 1 ~ "1st")) %>%
cbind(
aggregate(
dplyr::select(opioid_mesa2021) %>% mutate(overdosecaseCount = 1), ., sum) %>%
mutate(overdosecaseCount = replace_na(overdosecaseCount, 0))) %>%
dplyr::select(label,Risk_Category, overdosecaseCount)# Plot kernel density and risk predictions map for 2021 heroin overdose cases
rbind(overdose_KDE_sf, overdose_risk_sf) %>%
na.omit() %>%
gather(Variable, Value, -label, -Risk_Category, -geometry) %>%
ggplot() +
geom_sf(aes(fill = Risk_Category), colour = NA) +
#geom_sf(data = sample_n(opioid_mesa2021, 1201), size = .5, colour = "black") +
facet_wrap(~label, ) +
scale_fill_viridis(discrete = TRUE) +
labs(title="Comparison of Kernel Density and Risk Predictions",
subtitle="2021 heroin overdose",
caption = "Figure 10") +
theme_minimal() +
theme(panel.grid = element_blank(),
axis.text = element_blank(),
axis.title = element_blank(),
axis.ticks = element_blank(),
plot.caption = element_text(hjust = 0))
For the lower risk categories, our risk predictions tend to overpredict which as seen in the bar graph in Figure 11.
# Create custom color palette for bar chart
custom_colors <- c("Kernel Density" = "#587bb0",
"Risk Predictions" = "#8a9ab1")
# Plot bar chart of percentage of 2021 overdose per risk category for kernel density vs risk prediction models
rbind(overdose_KDE_sf, overdose_risk_sf) %>%
st_drop_geometry() %>%
na.omit() %>%
gather(Variable, Value, -label, -Risk_Category) %>%
group_by(label, Risk_Category) %>%
summarize(overdosecaseCount = sum(Value)) %>%
ungroup() %>%
group_by(label) %>%
mutate(Pcnt_of_test_set_crimes = overdosecaseCount / sum(overdosecaseCount)) %>%
ggplot(aes(Risk_Category,Pcnt_of_test_set_crimes)) +
geom_bar(aes(fill=label), position="dodge", stat="identity") +
scale_fill_manual(values = custom_colors, name = "Model") +
labs(title = "Risk prediction vs. Kernel density, 2021 overdose",
y = "% of 2021 Overdose (per model)",
x = "Risk Category",
caption = "Figure 11") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, vjust = 0.5),
panel.grid = element_blank(),
plot.caption = element_text(hjust = 0))
Model Application to Dashboard
Mapping current resources
We now have the risk predictions of heroin overdose based on the model above. However, this was not the only consideration when deciding on the placement of the Narcan vending machines. The location of existing resources should also be factored in so as to prioritize areas which are less accessible to current emergency services. In cases of opioid overdose, the lack of oxygen for three to five minutes can be life-threatening, making the timely availability of Narcan critical.
To optimize life-saving interventions, we identified locations within a 500-meter (approximately 1/3 mile) radius of medical facilities (such as hospitals, clinics, and pharmacies), fire stations, and police stations. These areas are likely to receive prompt emergency responses, with access to Narcan and personnel actively involved in drug-related and overdose incidents. A 1/3-mile buffer was created around all medical facilities, police stations, and fire stations. As a result, we have positioned the Narcan vending machines outside these zones, ensuring that individuals in areas with limited or delayed access to emergency services can still obtain free Narcan nasal spray when needed, potentially saving lives.
# Getting Open Street Map data for current narcan resources
## Medical Facilities
medical_mesa <- opq(bbox = "Mesa") %>%
add_osm_feature(key = "healthcare", value = c("clinic", "hospital", "pharmacy")) %>%
osmdata_sf()
## Police Stations
police_mesa <- opq(bbox = "Mesa") %>%
add_osm_feature(key = "amenity", value = c("police")) %>%
osmdata_sf()
## Fire Stations
fire_station <- opq(bbox = "Mesa") %>%
add_osm_feature(key = "amenity", value = c("fire_station")) %>%
osmdata_sf()
## Extract points for each resource
medical_mesa_points <- medical_mesa$osm_polygons %>%
st_set_crs(4326) %>%
dplyr::select(c("osm_id")) %>%
st_centroid() %>%
st_transform("ESRI:102249")
police_mesa_points <- police_mesa$osm_polygons %>%
st_set_crs(4326) %>%
st_make_valid() %>%
dplyr::select(c("osm_id")) %>%
st_centroid() %>%
st_transform("ESRI:102249")
fire_mesa_points <- fire_station$osm_polygons %>%
st_set_crs(4326) %>%
st_make_valid() %>%
dplyr::select(c("osm_id")) %>%
st_centroid() %>%
st_transform("ESRI:102249")
# Crop the points based on the Mesa Boundary
police_mesa <- police_mesa_points %>%
filter(sapply(st_intersects(police_mesa_points, mesaBoundaryNoAirport), length) > 0)
medical_mesa <- medical_mesa_points %>%
filter(sapply(st_intersects(medical_mesa_points, mesaBoundaryNoAirport), length) > 0)
fire_mesa <- fire_mesa_points %>%
filter(sapply(st_intersects(fire_mesa_points, mesaBoundaryNoAirport), length) > 0)
# Buffer of 500 meters for each resource
medicalBuffer <- st_buffer(medical_mesa, 500)
policeBuffer <- st_buffer(police_mesa, 500)
fireBuffer <- st_buffer(fire_mesa, 500)
# Combine all the buffers
current_resource_buffer <- rbind(medicalBuffer, policeBuffer, fireBuffer)# Map out current resource buffersbuffers
buffer_current_resources <- overdose_risk_sf %>%
na.omit() %>%
gather(Variable, Value, -label, -Risk_Category, -geometry) %>%
ggplot() +
geom_sf(aes(fill = Risk_Category), colour = NA) +
scale_fill_viridis(discrete = TRUE) +
geom_sf(data = policeBuffer, fill = "#98cbe3", color = NA, alpha = 0.7) +
geom_sf(data = medicalBuffer, fill = "#cfcfcf", color = NA, alpha = 0.7) +
geom_sf(data = fireBuffer, fill = "#e5bf68", color = NA, alpha = 0.7) +
labs(title="500m buffer around current resources",
caption = "Figure 12") +
theme_minimal() +
theme(panel.grid = element_blank(),
axis.text = element_blank(),
axis.title = element_blank(),
axis.ticks = element_blank())#Export files for creation of dashbobard mockup on Javascript
overdose_risk_sf_geojson <- overdose_risk_sf %>%
st_transform(crs = 4326)
current_resource_buffer_geojson <- current_resource_buffer %>%
st_transform(crs = 4326)
policeBuffer_geojson <- policeBuffer %>%
st_transform(crs = 4326)
medicalBuffer_geojson <- medicalBuffer %>%
st_transform(crs = 4326)
fireBuffer_geojson <- fireBuffer %>%
st_transform(crs = 4326)
st_write(overdose_risk_sf_geojson, "MesaHeroinResources.geojson", append = FALSE, quiet = TRUE)
st_write(current_resource_buffer_geojson, "ResourceBuffer.geojson", append = FALSE, quiet = TRUE)
st_write(policeBuffer_geojson, "policeBuffer.geojson", append = FALSE, quiet = TRUE)
st_write(medicalBuffer_geojson, "medicalBuffer.geojson", append = FALSE, quiet = TRUE)
st_write(fireBuffer_geojson, "fireBuffer.geojson", append = FALSE, quiet = TRUE)Mapping the gaps for vending machine placement
Overlaying the buffers of present resources over the risk prediction map, one can identify grid cells that have relatively higher predicted risk of overdose but not within a 1/3 mile radius of a current resource. Working within the limit of seven vending machines to be rolled out, we identified the seven most pressing grid cells to deploy the vending machines in.
# Map out current resources on risk prediction map base layer to identify grid cells to place vending machines
tmap_mode("view")
overdose_risk_sf_geojson <- overdose_risk_sf_geojson %>%
mutate(Unique_ID = row_number())
tm_shape(overdose_risk_sf_geojson %>% na.omit()) +
tm_fill(
"Risk_Category",
palette = "viridis",
alpha = 0.5,
title = "Risk Category",
popup.vars = c("Unique ID" = "Unique_ID", "Risk Category" = "Risk_Category")
) +
tm_shape(current_resource_buffer_geojson) +
tm_fill(fill = "#efefef", alpha = 0.3, title = "Resource Buffer") +
tm_layout(
title = "500m buffer around current resources",
title.size = 1.5,
title.position = c("center", "top"),
frame = FALSE,
legend.outside = TRUE,
legend.title.size = 0.8,
legend.text.size = 0.7
)Placing a vending machine in the centroid of an identified grid cell would be ideal, as this maximizes coverage within that cell and ensures the vending machine is equally accessible from all directions within the cell. However, this may not always be practical in reality, thus suitable public locations as close to the centroid as possible such as a park, social service organization or retail center would be the proposed sites for the vending machines.
#Select grid cells for vending machine placement
selected_polygons <- overdose_risk_sf_geojson %>%
filter(Unique_ID %in% c(140, 172, 165, 173, 210, 235, 242))
# Calculate centroids
centroids <- st_centroid(selected_polygons)
# Extract coordinates of centroids
centroid_coords <- st_coordinates(centroids)
# Create dataframe of coordinates of selected polygon centroids and export as file
result <- data.frame(
Unique_ID = selected_polygons$Unique_ID,
Longitude = centroid_coords[, "X"],
Latitude = centroid_coords[, "Y"]
)
centroid_points_sf <- result %>%
st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326)
st_write(centroid_points_sf, "centroid_points.geojson", append = FALSE, quiet = TRUE)
# Plot out vending machine placement centroid points
tm_shape(overdose_risk_sf) +
tm_fill("Risk_Category", palette = "viridis", title = "Risk Category") +
tm_shape(current_resource_buffer) +
tm_fill(fill = "#efefef", alpha = 0.7, title = "Resource Buffer") +
tm_shape(centroid_points_sf) +
tm_dots(size = 0.2, col = "red", title = "Centroids") +
tm_layout(
title = "Centroids of Selected Polygons",
title.size = 1.5,
title.position = c("center", "top"),
frame = FALSE,
legend.outside = TRUE,
legend.title.size = 0.8,
legend.text.size = 0.7
)Dashboard Use
The NARVEND dashboard will use the risk prediction map as the basemap. City officials can also click on individual vending machines to check the live stock count of the Narcan nasal sprays in each vending machine, enabling the monitoring of usage and necessity for restocking. This risk prediction map will update every quarter based on the most up-to-date risk factor values, such as the most recent crime counts, a risk factor category that is particularly highly correlated with heroin overdose counts. By comparing the most up-to-date risk predictions with current vending machine usage levels, the dashboard enables city health officials to adjust vending machine placements as needed, ensuring they can be relocated to more suitable areas to meet evolving demand.
Additionally, the 2020 heroin overdose dataset was selected to train the model due to the surge in cases during the COVID-19 pandemic, which highlighted vulnerable individuals at increased risk of overdose. However, this data set will be used only for the first three years of predictions. Afterwards, reported overdose cases from the previous year will be used to generate predictions for the following year.
Conclusion
Overall, the geospatial risk model is applied to the NARVEND dashboard, along with mapping current resources and monitoring live vending machine usage, to support health officials in identifying the most optimal locations for the Narcan vending machines.
City health officials should keep in mind that while the dashboard can assist in optimizing the placement of vending machines, its overall effectiveness hinges on the chosen locations, as well as any future changes, being clearly communicated and well publicized to the surrounding neighborhoods and broader community. This can be achieved through leveraging multiple outreach channels, such as local organizations, social media platforms and healthcare providers, to raise awareness about the vending machines and ensure that individuals in need can easily locate and access them.
As highlighted in the model evaluation section, the use of biased observational data of reported overdose cases is a limitation of the model that warrants further consideration of reporting practices and potential systemic inequities in overdose monitoring. Future efforts could involve integrating additional data sources, such as death certificates, harm reduction outreach reports, and qualitative data from interviews with stakeholders on the ground. This multi-source approach could help mitigate reporting bias and enhance the evaluation of risk predictions to produce more accurate and equitable representations of latent heroin overdose risk.
Heroin overdose is a serious issue in Mesa and while broader structural change is indisputably needed to address the root cause of this, increasing the availability and accessibility of Narvan can help to prevent devastating consequences should an overdose occur. With limited resources, the NARVEND dashboard is a unique and valuable resource for health officials to utilize to optimize allocation of these life-saving nasal sprays to areas with high predicted overdose risk but lacking accessibility, maximizing the number of lives that can be saved.
References
American Addiction Centers. (2022, July 26). Highest Drug Use By City. American Addiction Centers. https://americanaddictioncenters.org/blog/substance-abuse-by-city.
City of Mesa Data Portal. (2023, June 14). Unsheltered Point-In-Time (PIT) Count 2022 Details Mesa Only. Mesaaz.gov. https://citydata.mesaaz.gov/Community-Services/Unsheltered-Point-In-Time-PIT-Count-2022-Details-M/efjd-c5mi/about_data.
Mesa Data Hub. (2024a, September 30). Parks Locations And Amenities. Data.mesaaz.gov. https://data.mesaaz.gov/Parks-Recreation-and-Community-Facilities/Parks-Locations-And-Amenities/djym-pkpp/about_data.
Mesa Data Hub. (2024b, November 21). Fire and Medical Opioid Overdose Incidents. Data.mesaaz.gov. https://data.mesaaz.gov/Fire-and-Medical/Fire-and-Medical-Opioid-Overdose-Incidents/qufy-tzv6/about_data.
Mesa Data Hub. (2024c, December 3). Police Incidents. Mesaaz.gov. https://data.mesaaz.gov/Police/Police-Incidents/39rt-2rfj/about_data.