Geospatial Risk Modeling - Predicting Heroin Overdose in Mesa, AZ

Introduction

Of the few public sector machine learning algorithms in use today, none are as prevalent as the collection of algorithms that we refer to as “Predictive Response.” Predictive Response came about from city authorities asking questions regarding planning and resource allocation, and as such, has the ability to forecast specific events, where the event will happen, the city’s response to said event, as well as outcomes of responses. The context of our study is Mesa, Arizona, where we will develop a predictive response algorithm to predict heroin overdose incidents before they happen in order to better allocate overdose prevention services and resources. In 2018, the Mesa Fire and Medical Department Emergency Medical Support program recorded 548 incidents of overdoses, which is equivalent to almost 2 incidents per day. With less than 500 sworn-in Mesa Fire and Medical Department personnel to patrol Mesa’s 133 square mile area, what is the best way for the MFMF to allocate these personnel across time and space? In this report, we will be creating our own geospatial risk model, a type of regression model, to predict where incidents of overdoses will likely occur in Mesa. Using overdoses as our dependent variable, our hypothesis is that overdose risk is a function of exposure to a number of geospatial risk factors, with the assumption that exposure to said risk factors increases the rates of heroin-related overdoses. The risk factors that we have selected to test and implement in our model are: presence of graffiti, abandoned buildings, street light service requests, homeless person sightings, mentally ill person sightings, liquor sites (retailers or bars), tobacco retailers, subsidized housing, sex-related crime, drug-related crime, theft_related crime, and suicides. After validating our model through cross-validation, we will assess our model’s accuracy and generalizability across different contexts (space/time and socio-economic) in order to compare how useful our model is when compared to Mesa’s current approach. Depending on our results, we will choose whether or not to develop our predictive model into an app, Mesa Guardian, that could be used by the Mesa Fire and Medical Department to better allocate overdose prevention resources and services.

Learn more about Mesa Guardian here

Read in data from Mesa

overdoses <- read.csv("Data/Fire_and_Medical_Opioid_Overdose_Incidents.csv") %>%
  filter(Year == 2018) %>%
  na.omit() %>%
  st_as_sf(coords = c("Longitude","Latitude"), crs = 4326) %>%
  st_transform('EPSG:3857')

mesaBoundary <- st_read("Data/Mesa Census Tracts To City Boundary.geojson") %>%
  st_transform('EPSG:3857')

councilDistricts <- st_read("Data/CouncilDistrict.geojson") %>%
  st_transform('EPSG:3857')

The code block above imports 2018 overdose data, as well as Mesa’s city and council district boundaries, from Mesa’s Open Data Portal.

Visualizing Point Data

# filter out points that don't lie within Mesa boundary 
overdoses <- overdoses %>% mutate(in_bounds = lengths(st_within(overdoses,mesaBoundary)))
overdoses.in.boundary <- overdoses%>%
  filter(in_bounds == 1)

# uses grid.arrange to organize independent plots
grid.arrange(ncol=2,
ggplot() + 
  geom_sf(data = mesaBoundary) +
  geom_sf(data = overdoses.in.boundary, colour="red", size=0.1, show.legend = "point") +
  labs(title= "Overdoses, Mesa - 2018") +
  mapTheme(title_size = 14),

ggplot() + 
  geom_sf(data = mesaBoundary, fill = "grey40") +
  stat_density2d(data = data.frame(st_coordinates(overdoses.in.boundary)), 
                 aes(X, Y, fill = ..level.., alpha = ..level..),
                 size = 0.01, bins = 40, geom = 'polygon') +
  scale_fill_viridis() +
  scale_alpha(range = c(0.00, 0.35), guide = FALSE) +
  labs(title = "Density of Overdoses") +
  mapTheme(title_size = 14) + theme(legend.position = "none"))

Figure 1. The map on the left illustrates my outcome of interest, overdoses, in point form within Mesa’s city boundaries. The map on the right shows the density of overdoses in Mesa, with darker regions meaning less relative density compared to lighter regions. We can assume that this dataset is missing some incidents of overdoses because, for an overdose to appear in the data, it must be observed by a fire or medical team. Since not all overdoses can be observed and since fire and medical teams may be allocated more fervently in some communities than others (West Mesa), our model’s overdose outcomes may be influenced by selection bias and should be observed more closely.

Creating a fishnet grid

## using {sf} to create the grid
fishnet <- 
  st_make_grid(mesaBoundary,
               cellsize = 750, 
               square = TRUE) %>%
  .[mesaBoundary] %>%            # <- MDH Added
  st_sf() %>%
  mutate(uniqueID = rownames(.))

The code block above creates a grid cell lattice, called a fishnet, with 750 ft by 750 ft grid cells laid over the city of Mesa. A good way of understanding this process is relating it to elevation; overdoses cluster in space and overdose risk gets smaller as you move outward from cluster hotspots, just as elevation dips from the tops of mountains to the bottom of valleys.

Aggregate points to the fishnet

## add a value of 1 to each overdose, sum them with aggregate
overdoses_net <- 
  dplyr::select(overdoses.in.boundary) %>% 
  mutate(countOverdoses = 1) %>% 
  aggregate(., fishnet, sum) %>%
  mutate(countOverdoses = replace_na(countOverdoses, 0),
         uniqueID = rownames(.),
         cvID = sample(round(nrow(fishnet) / 8), 
                       size=nrow(fishnet), replace = TRUE))

ggplot() +
  geom_sf(data = overdoses_net, aes(fill = countOverdoses), color = NA) +
  scale_fill_viridis() +
  labs(title = "Count of Overdoses for the fishnet") +
  mapTheme()

xx <- mapview::mapview(overdoses_net, zcol = "countOverdoses")
yy <- mapview::mapview(mutate(overdoses.in.boundary, ID = seq(1:n())))
xx + yy

Figure 2. This map joins overdose counts to the fishnet to return the count of overdoses per 750 ft by 750 ft grid cell. This illustrates the clustered spatial process of overdoses in Mesa. We also split the number of grids, a total of 3639, by 26 and randomly assign each group a dummy id that will be used later in cross-validation. Because we are dealing with relatively rare events (overdoses), we make the fishnet grid size 750 ft by 750 ft (relatively small) in order to capture a realistic count of overdoses in each grid.

Modeling Spatial Features

# load and wrangle risk factor datasets
drug_crime <- 
  read.socrata("https://data.mesaaz.gov/Police/Police-Incidents/39rt-2rfj") %>%
    filter(city == "MESA") %>%
    filter(report_year == 2018) %>%
    filter(crime_type == "DUI - LIQUOR - DRUGS" | crime_type == "DRUG PARAPHERNALIA-POSSESS-USE" | crime_type == "DANGEROUS DRUG-POSS-USE") %>%
    dplyr::select(Y = latitude, X = longitude) %>%
    na.omit() %>%
    st_as_sf(coords = c("X", "Y"), crs = 4326, agr = "constant") %>%
    st_transform(st_crs(fishnet)) %>%
    mutate(Legend = "Crime_Drugs")


suicide_crime <- 
  read.socrata("https://data.mesaaz.gov/Police/Police-Incidents/39rt-2rfj") %>%
    filter(city == "MESA") %>%
    filter(report_year == 2018) %>%
    filter(crime_type == "DEATH-SUICIDE" | crime_type == "SUICIDE ATTEMPT") %>%
    dplyr::select(Y = latitude, X = longitude) %>%
    na.omit() %>%
    st_as_sf(coords = c("X", "Y"), crs = 4326, agr = "constant") %>%
    st_transform(st_crs(fishnet)) %>%
    mutate(Legend = "Crime_Suicides")


theft_crime <- 
  read.socrata("https://data.mesaaz.gov/Police/Police-Incidents/39rt-2rfj") %>%
    filter(city == "MESA") %>%
    filter(report_year == 2018) %>%
    filter(crime_type == "THEFT-CONTROL LOST PROPERTY" | crime_type == "SHOPLIFTING-REMOVAL OF GOODS" | crime_type == "THEFT-CONTROL PROPERTY" | crime_type == "THEFT-MEANS OF TRANSPORTATION" | crime_type == "SHOPLIFTING-CONCEALMENT" | crime_type == "BURGLARY, 3RD DEGREE" | crime_type == "BURGLARY 1ST DEGREE" | crime_type == "BURGLARY 2ND DEGREE") %>%
    dplyr::select(Y = latitude, X = longitude) %>%
    na.omit() %>%
    st_as_sf(coords = c("X", "Y"), crs = 4326, agr = "constant") %>%
    st_transform(st_crs(fishnet)) %>%
    mutate(Legend = "Crime_Thefts")


sex_crime <- 
  read.socrata("https://data.mesaaz.gov/Police/Police-Incidents/39rt-2rfj") %>%
    filter(city == "MESA") %>%
    filter(report_year == 2018) %>%
    filter(crime_type == "PUBLIC SEXUAL INDEC-CONTACT" | crime_type == "SEXUAL ASSAULT - DV" | crime_type == "SEXUAL ASSAULT") %>%
    dplyr::select(Y = latitude, X = longitude) %>%
    na.omit() %>%
    st_as_sf(coords = c("X", "Y"), crs = 4326, agr = "constant") %>%
    st_transform(st_crs(fishnet)) %>%
    mutate(Legend = "Crime_Sex")


mental_ill_assistance <- 
  read.socrata("https://data.mesaaz.gov/Police/Police-Incidents/39rt-2rfj") %>%
    filter(city == "MESA") %>%
    filter(report_year == 2018) %>%
    filter(crime_type == "ASSIST-MENTALLY ILL PERSON") %>%
    dplyr::select(Y = latitude, X = longitude) %>%
    na.omit() %>%
    st_as_sf(coords = c("X", "Y"), crs = 4326, agr = "constant") %>%
    st_transform(st_crs(fishnet)) %>%
    mutate(Legend = "Mentally_Ill")


graffiti <- 
  read.socrata("https://data.mesaaz.gov/Transportation/Transportation-Graffiti/9spb-749m") %>%
    filter(year_reported == 2018) %>%
    dplyr::select(Y = lat, X = lon) %>%
    na.omit() %>%
    st_as_sf(coords = c("X", "Y"), crs = 4326, agr = "constant") %>%
    st_transform(st_crs(fishnet)) %>%
    mutate(Legend = "Graffiti")


streetLightsOut <- 
  read.socrata("https://data.mesaaz.gov/Transportation/Transportation-Work-Orders-Public-View/ay68-xuca") %>%
    mutate(year = substr(initiated_date,1,4)) %>% filter(year == "2018") %>%
    filter(description == "Contract Streetlight Repair" | description == "Repair Streetlight - Residential" | description == "Repair Streetlight") %>%
    dplyr::select(Y = latitude, X = longitude) %>%
    na.omit() %>%
    st_as_sf(coords = c("X", "Y"), crs = 4326, agr = "constant") %>%
    st_transform(st_crs(fishnet)) %>%
    mutate(Legend = "Street_Lights_Out")


homeless <-
  read.socrata("https://data.mesaaz.gov/Community-Services/Unsheltered-Point-in-Time-PIT-Count/jagk-fkkw") %>%
    filter(reporting_year == 2018) %>%
    filter(city == "Mesa") %>%
    dplyr::select(Y = latitude, X = longitude) %>%
    na.omit() %>%
    st_as_sf(coords = c("X", "Y"), crs = 4326, agr = "constant") %>%
    st_transform(st_crs(fishnet)) %>%
    mutate(Legend = "Homeless")


abandonedBuildings <- 
  read.csv("Data/City_Owned_Property.csv") %>% 
    filter(City.of.Mesa.Description == "Vacant (ADOT remnant)" | City.of.Mesa.Description == "Vacant") %>%
    dplyr::select(Y = Latitude, X = Longitude) %>%
    na.omit() %>%
    st_as_sf(coords = c("X", "Y"), crs = 4326, agr = "constant") %>%
    st_transform(st_crs(fishnet)) %>%
    mutate(Legend = "Abandoned_Buildings")


smokeshops <- 
  read.csv("Data/AZ.csv") %>% 
    filter(sg_c__sub_category == "Tobacco Stores") %>%
    dplyr::select(Y = sg_c__latitude, X = sg_c__longitude) %>%
    na.omit() %>%
    st_as_sf(coords = c("X", "Y"), crs = 4326, agr = "constant") %>%
    st_transform(st_crs(fishnet)) %>%
    mutate(Legend = "Tobacco_Retail")


liquor <- 
  read.csv("Data/AZ.csv") %>% 
    filter(sg_c__sub_category == "Beer, Wine, and Liquor Stores" | sg_c__sub_category == "Drinking Places (Alcoholic Beverages)") %>%
    dplyr::select(Y = sg_c__latitude, X = sg_c__longitude) %>%
    na.omit() %>%
    st_as_sf(coords = c("X", "Y"), crs = 4326, agr = "constant") %>%
    st_transform(st_crs(fishnet)) %>%
    mutate(Legend = "Alcohol_Site")


subsidized <- 
  st_read("Data/Mesa CARES - Rent and Mortgage Assistance.geojson") %>%
  dplyr::select(geometry) %>%
  st_transform('EPSG:3857') %>%
  na.omit() %>%
  mutate(Legend = "Subsidized_Housing") 


subsidized = subsidized[!st_is_empty(subsidized),drop=FALSE]


downtown <- 
  st_read("Data/DowntownMesa.geojson") %>%
  st_transform(st_crs(fishnet)) 

# load in and wrangle census data 
census_api_key("ea7bd3babfc0da6a4caa6d2536b66403a88929f2", overwrite = TRUE)

neighborhoods <- 
  get_acs(geography = "tract", variables = c("B01001_001E","B01001A_001E"), 
          year = 2018, state=04, county=013, geometry=T, output='wide') %>%
  st_transform('EPSG:3857') %>% 
  dplyr::rename(TotalPop = B01001_001E,
                NumberWhites = B01001A_001E) %>%
  dplyr::select(TotalPop, NumberWhites, GEOID) %>%
  mutate(percentWhite = NumberWhites / TotalPop,
         raceContext = ifelse(percentWhite > .5, "Majority_White", "Majority_Non_White")) %>%
  na.omit()


# filter out neighborhoods that don't lie within Mesa boundary 
neighborhoods <- st_filter(neighborhoods,mesaBoundary)

The code block above imports all of the risk factors that we will use to inform our predictive model: abandoned buildings, street light service requests, tobacco retailers, liquor sites, sex crime, drug crime, theft crime, subsidized housing, homeless person, mentally-ill person, grafitti, and suicides. We also make a variable called “downtown” to delineate Mesa between Downtown and the rest of Mesa, and import Census tract data to delineate Mesa between different neighborhoods and socio-economic contexts.

Feature Engineering

Count of Risk Factors by Grid Cell

# get count of each risk factor by grid cell in fishnet 
vars_net <- 
  rbind(abandonedBuildings,streetLightsOut,homeless,
        mental_ill_assistance, graffiti, sex_crime,
        theft_crime, suicide_crime, drug_crime, smokeshops,
        liquor, subsidized) %>%
  st_join(., fishnet, join=st_within) %>%
  st_drop_geometry() %>%
  group_by(uniqueID, Legend) %>%
  summarize(count = n()) %>%
    full_join(fishnet) %>%
    spread(Legend, count, fill=0) %>%
    st_sf() %>%
    dplyr::select(-`<NA>`) %>%
    na.omit() %>%
    ungroup()

# turn vars_net into long form and filter by risk factor to append to map  
vars_net.long <- 
  gather(vars_net, Variable, value, -geometry, -uniqueID)

vars <- unique(vars_net.long$Variable)
mapList <- list()

# map of risk factors
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="") +
      labs(title=i) +
      mapTheme()}

do.call(grid.arrange,c(mapList, ncol=4, top="Risk Factors by Fishnet"))

Figure 3. This is a multiple map of counts of each risk factor per each 750 ft by 750 ft grid cell. Here we see that some risk factors such as presence of graffiti and street light service requests have relatively large counts of incidents when compared to others like abandoned buildings and tobacco retailers. We also see that some risk factors are somewhat evenly distributed across the city - such as sex crime and street light service requests - while others are more concentrated in certain areas, like suicides and drug crime. Risk factors that are concentrated and not as spread out may introduce more selection bias into my model, since the increased counts in these areas may be derived from the factors being identified more fervently in some communities than others.

Nearest Neighbor Features

The code blocks below create nearest neighbor features for each of our model’s different risk factors and maps them according to the fishnet. Average nearest neighbor features are created by converting vars_net grid cells to centroid points then measuring to 3 k risk factor points.

# create nearest neighbor features for risk factors by fishnet 
st_c <- st_coordinates
st_coid <- st_centroid

vars_net <-
  vars_net %>%
    mutate(
      Abandoned_Buildings.nn =
        nn_function(st_c(st_coid(vars_net)), st_c(abandonedBuildings),3),
      Street_Lights_Out.nn =
        nn_function(st_c(st_coid(vars_net)), st_c(streetLightsOut),3),
      Graffiti.nn =
        nn_function(st_c(st_coid(vars_net)), st_c(graffiti),3),
      Tobacco_Retail.nn =
        nn_function(st_c(st_coid(vars_net)), st_c(smokeshops),3),
      Liquor.nn =
        nn_function(st_c(st_coid(vars_net)), st_c(liquor),3),
      Homeless.nn =
        nn_function(st_c(st_coid(vars_net)), st_c(homeless),3),
      Mentally_Ill.nn =
        nn_function(st_c(st_coid(vars_net)), st_c(mental_ill_assistance),3),
      Subsidized_Housing.nn =
        nn_function(st_c(st_coid(vars_net)), st_c(subsidized),3),
      Theft_Crime.nn =
        nn_function(st_c(st_coid(vars_net)), st_c(theft_crime),3),
      Sex_Crime.nn =
        nn_function(st_c(st_coid(vars_net)), st_c(sex_crime),3),
      Drug_Crime.nn =
        nn_function(st_c(st_coid(vars_net)), st_c(drug_crime),3),
      Suicide.nn =
        nn_function(st_c(st_coid(vars_net)), st_c(suicide_crime),3))

# turn vars_net into long form and filter by risk factor to append to map 
vars_net.long.nn <- 
  dplyr::select(vars_net, ends_with(".nn")) %>%
    gather(Variable, value, -geometry)

vars <- unique(vars_net.long.nn$Variable)
mapList <- list()

# map nearest neighbor features for risk factors by fishnet 
for(i in vars){
  mapList[[i]] <- 
    ggplot() +
      geom_sf(data = filter(vars_net.long.nn, Variable == i), aes(fill=value), colour=NA) +
      scale_fill_viridis(name="") +
      labs(title=i) +
      mapTheme()}

do.call(grid.arrange,c(mapList, ncol = 4, top = "Nearest Neighbor risk Factors by Fishnet"))

Figure 4. This is a multiple map of average nearest neighbor distance of each risk factor. Where our first feature engineering approach of counts of risk factors per grid cell produced a rigid spatial scale of exposure, average nearest neighbor distance hypothesizes a smoother exposure relationship. First, we converted the risk factor grid cells to centroid points, and then we measured the nearest neighbor to 3 k risk factor points. Now we have another way of seeing how the risk factors are distributed across space. Two commonalities across these maps are that the downtown area in the bottom right hand corner has high nearest neighbor distances and that the high-density overdose area in Figure 1 has low nearest neighbor distances, suggesting that most of our risk factors tend to occur where overdoses occur.

Measure Distance to One Point

In the code block below, we measure the distances of the risk factors to the centroid of Mesa’s downtown (high business activity) area and map them.

# measure distance of risk factors to the centroid of Mesa's Downtown 
downtownPoint <-
  downtown %>%
  st_centroid()

vars_net$downtownDistance =
  st_distance(st_centroid(vars_net),downtownPoint) %>%
  as.numeric() 

# map distances 
ggplot() + 
  geom_sf(data = vars_net, aes(fill=downtownDistance)) +
  scale_fill_viridis(name="downtownDistance") +
  labs(title="Euclidean distance to Downtown Mesa") +
  mapTheme()

Figure 5. This map shows distance of risk factors to the centroid of Mesa’s downtown area. What we see here is a confirmation of our findings in Figure 4, in that our risk factors generally occur outside of the downtown area and tend to be organized in space in a similar fashion to overdoses (in West Mesa).

Delineate Between Council Districts

In the code block below, we join our overdose data with our risk factor data and separate Mesa into council districts and neighborhoods.

# join crime_net and vars_net to make regression-ready dataframe 
final_net <-
  left_join(overdoses_net, st_drop_geometry(vars_net), by="uniqueID") 

# delineate between council districts to see how risk factors are spatially oriented 
final_net <-
  st_centroid(final_net) %>%
    st_join(dplyr::select(neighborhoods, GEOID)) %>%
    st_join(dplyr::select(councilDistricts, DISTRICT)) %>%
      st_drop_geometry() %>%
      left_join(dplyr::select(final_net, geometry, uniqueID)) %>%
      st_sf() %>%
  na.omit()

final_net$GEOID <- as.factor(final_net$GEOID)

grid.arrange(ncol=2,
ggplot() + 
  geom_sf(data = final_net, aes(fill=DISTRICT)) + 
  scale_fill_viridis() +
  theme(legend.position = "none") + 
  labs(title = "Council Districts, Mesa") +
  mapTheme(title_size = 20),
ggplot() + 
  geom_sf(data = final_net, aes(fill=GEOID)) +
  theme(legend.position = "none") + 
  labs(title = "Neighborhoods, Mesa") + 
  mapTheme(title_size = 20))

Figures 6. These maps shows how Mesa is broken down by council district and neighborhood. There are six council districts in all, each designated by 2020 Census data, and there are 139 different neighborhoods.

Moran’s I

# create neighbor list and spatial weights matrix 
final_net.nb <- poly2nb(as_Spatial(final_net), queen=TRUE)
final_net.weights <- nb2listw(final_net.nb, style="W", zero.policy=TRUE)

The code block above creates a neighbor list, final_net.nb, and a spatial weights matrix, final_net.weights using the “queen” contiguity, meaning every grid cell is related to the 8 grid cells surrounding it.

Local Moran’s I Statistics

# local Moran's I test for count of overdoses and spatial weights
local_morans <- localmoran(final_net$countOverdoses, final_net.weights, zero.policy=TRUE) %>% 
  as.data.frame()

final_net.localMorans <- 
  cbind(local_morans, as.data.frame(final_net)) %>% 
  st_sf() %>%
  dplyr::select(Overdose_Count = countOverdoses, 
                Local_Morans_I = Ii, 
                P_Value = `Pr(z != E(Ii))`) %>%
  mutate(Significant_Hotspots = ifelse(P_Value <= 0.001, 1, 0)) %>%
  gather(Variable, Value, -geometry)

  
vars <- unique(final_net.localMorans$Variable)
varList <- list()

# map it 
for(i in vars){
  varList[[i]] <- 
    ggplot() +
      geom_sf(data = filter(final_net.localMorans, Variable == i), 
              aes(fill = Value), colour=NA) +
      scale_fill_viridis(name="") +
      labs(title=i) +
      mapTheme() + theme(legend.position="bottom")}

do.call(grid.arrange,c(varList, ncol = 4, top = "Local Morans I statistics, Overdoses"))

Figure 7. This multiple map shows the Local Moran’s I statistics for overdoses in Mesa. The null hypothesis of a Local Moran’s I test is that the count of overdoses at a given grid cell is randomly distributed relative to its 8 surrounding grid cells. Here, “Significant Hotspots” are considered those grid cells with higher local counts than what otherwise would be expected under random conditions (p-values are less than or equal to 0.001). This local insight is instrumental for the generalizability of our model, as they can be altered into other spatial features to control for more local spatial processes.

Adjusting P-Value To Observe Appropriate Scale for Local Spatial Processes

# look deeper into local hotspots by varying p-value of Moran's I test 
final_net <-
  final_net %>% 
  mutate(overdose.isSig = ifelse(localmoran(final_net$countOverdoses, 
                                 final_net.weights, zero.policy=TRUE)[,5] <= 0.001, 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))

ggplot() + 
  geom_sf(data = final_net, aes(fill=overdose.isSig.dist)) +
  scale_fill_viridis(name="overdose.isSig.dist") +
  labs(title="Distance to highly significant overdose hotspots") +
  mapTheme()

Figure 8. This map explores local hotspots by changing the p-value included in the local Moran’s I test. The smaller the p-value, the more significant the clusters are, with a p-value of 0.001 conforming to significant overdose hotspots. What we see here is very similar to what we saw in Figure 5, with Downtown Mesa having the highest values for distance to significant hotspots and West Mesa having the lowest values.

Correlation Tests

The code block below executes correlation tests between each risk factor and our dependent variable of overdoses and plots them in scatterplot format.

# create multiple scatterplots of countOverdoses as a function of risk factors 
correlation.long <-
  st_drop_geometry(final_net) %>%
    dplyr::select(-uniqueID, -cvID, -DISTRICT, -GEOID) %>%
    gather(Variable, Value, -countOverdoses)

correlation.cor <-
  correlation.long %>%
    group_by(Variable) %>%
    summarize(correlation = cor(Value, countOverdoses, use = "complete.obs"))


ggplot(correlation.long, aes(x=Value, y=countOverdoses)) +
  geom_point(size = 0.1) +
  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 = 7, scales = "free_x", labeller=(Variable=label_wrap_gen(width=100))) +
  theme_classic() +
  theme(axis.text=element_text(size=14),
        axis.text.x=element_text(size=8),
        axis.title=element_text(size=16),
        title=element_text(size=8),
        legend.position="none",
        strip.text=element_text(size=10))+
  labs(title = "Overdose count as a function of risk factors") +
  plotTheme()

Figure 9. This is a collection of scatterplots with count of overdoses as a function of the different risk factors. Included in each of the scatterplots is the correlation coefficient, r. R can range from -1 to +1, with negative values representing negative linear relationships between overdose count and risk factor and vice versa.

Correlation Matrix

The code block below presents the correlation coefficients in an easier-to-read, matrix form.

correlation.sf <- final_net%>%
  dplyr::select(countOverdoses,Abandoned_Buildings,Alcohol_Site,Crime_Drugs,Crime_Sex,Crime_Suicides,Crime_Thefts,
                Graffiti,Homeless,Mentally_Ill,Street_Lights_Out,Subsidized_Housing,Tobacco_Retail,
                Abandoned_Buildings.nn, Street_Lights_Out.nn, Graffiti.nn, Tobacco_Retail.nn, Liquor.nn,
                Homeless.nn,Mentally_Ill.nn,Subsidized_Housing.nn,Theft_Crime.nn,Sex_Crime.nn,Drug_Crime.nn,
                Suicide.nn,downtownDistance,overdose.isSig, overdose.isSig.dist, geometry)

numericVars <- 
  select_if(st_drop_geometry(correlation.sf), is.numeric) %>% na.omit()

ggcorrplot(
  round(cor(numericVars), 1), 
  p.mat = cor_pmat(numericVars),
  type="lower",
  insig = "blank") +
  scale_fill_viridis() 

  labs(title = "Correlation across numeric variables") 

Figure 10. Using the ggmap package in R, we plotted relevant variables on the correlation matrix above. This matrix informed our variable selection for the final model. For example, Street_Lights_Out and Abandoned_Buildings appear to have little to no correlation with overdose count, so they were left out. On the other hand, overdose.isSig.dist and Suicide.nn have strong negative correlations and Graffiti and Crime_Drugs have strong positive correlations, so they were included in the model.

Histogram of Dependent Variable, Overdoses

Figure 11. This is a histogram of the dependent variable, overdoses, illustrating frequency distributions of count of overdoses per 400 ft by 400 ft grid cell. We can clearly see here that an overwhelming majority of our grid cells do not have any counts of overdoses, which makes sense considering how rare of an event overdoses generally are (2 per day in Mesa).

Modeling and Cross-Validation

# Make variable names for k-fold and LOGO-CV
reg.vars <- c("Drug_Crime.nn", "Graffiti.nn", "Liquor.nn", "Theft_Crime.nn", "Tobacco_Retail.nn",
              "Suicide.nn", "Sex_Crime.nn", "Homeless.nn", "Mentally_Ill.nn", "downtownDistance")

reg.ss.vars <- c("Drug_Crime.nn", "Graffiti.nn", "Tobacco_Retail.nn", "Liquor.nn", "Theft_Crime.nn",
                 "Suicide.nn", "Sex_Crime.nn", "Homeless.nn", "Mentally_Ill.nn", "downtownDistance",
                 "overdose.isSig", "overdose.isSig.dist")

This code block above selects the features that will be used in our four regressions. Figures 9 and 10 informed us about which risk factors we should use in our regression, as we took out risk factor features that had correlation coefficients that were relatively low when compared to other coefficients (any coefficient between -0.2 and +0.2 we tossed out). The reg.vars features are just risk factors, while the reg.ss.vars includes said risk factors and the spatial process features created from the local Moran’s I test.

# Make new cross-validate function 
crossValidateReal <- function(dataset, id, dependentVariable, indVariables) {
  
  allPredictions <- data.frame()
  cvID_list <- unique(dataset[[id]])
  
  for (i in cvID_list) {
    
    thisFold <- i
    cat("This hold out fold is", thisFold, "\n")
    
    fold.train <- filter(dataset, dataset[[id]] != thisFold) %>% as.data.frame() %>% 
      dplyr::select(id, geometry, indVariables, dependentVariable)
    fold.test  <- filter(dataset, dataset[[id]] == thisFold) %>% as.data.frame() %>% 
      dplyr::select(id, geometry, indVariables, dependentVariable)
    
    regression <-
      glm(paste0(dependentVariable,"~."), family = "poisson", 
          data = fold.train %>% 
            dplyr::select(-geometry, -id))
    
    thisPrediction <- 
      mutate(fold.test, Prediction = predict(regression, fold.test, type = "response"))
    
    allPredictions <-
      rbind(allPredictions, thisPrediction)
    
  }
  return(st_sf(allPredictions))
}

The code block above creates a function called crossValidate to estimate the different regressions.

# first two are for k-fold validation using risk factors and nearest neighbor features
reg.cv <- crossValidateReal(
  dataset = final_net,
  id = "cvID",
  dependentVariable = "countOverdoses",
  indVariables = reg.vars) %>%
    dplyr::select(cvID = cvID, countOverdoses, Prediction, geometry)

reg.ss.cv <- crossValidateReal(
  dataset = final_net,
  id = "cvID",
  dependentVariable = "countOverdoses",
  indVariables = reg.ss.vars) %>%
    dplyr::select(cvID = cvID, countOverdoses, Prediction, geometry)
  
# second two are for LOGO-CV, doing spatial cross-validation with previous two sets of features 
reg.spatialCV <- crossValidateReal(
  dataset = final_net,
  id = "GEOID",
  dependentVariable = "countOverdoses",
  indVariables = reg.vars) %>%
    dplyr::select(cvID = GEOID, countOverdoses, Prediction, geometry)

reg.ss.spatialCV <- crossValidateReal(
  dataset = final_net,
  id = "GEOID",
  dependentVariable = "countOverdoses",
  indVariables = reg.ss.vars) %>%
    dplyr::select(cvID = GEOID, countOverdoses, Prediction, geometry)

This code block crossvalidates our four regressions, and converts them to sf layers with both predicted and observed overdose counts. We do 139 (2374 grids/17) different k-folds because we have 139 different neighborhoods.

# make long form of observed and predicted counts and errors 
reg.summary <- 
  rbind(
    mutate(reg.cv,           Error = Prediction - countOverdoses,
                             Regression = "Random k-fold CV: Just Risk Factors"),
                             
    mutate(reg.ss.cv,        Error = Prediction - countOverdoses,
                             Regression = "Random k-fold CV: Spatial Process"),
    
    mutate(reg.spatialCV,    Error = Prediction - countOverdoses,
                             Regression = "Spatial LOGO-CV: Just Risk Factors"),
                             
    mutate(reg.ss.spatialCV, Error = Prediction - countOverdoses,
                             Regression = "Spatial LOGO-CV: Spatial Process")) %>%
    st_sf()

This code block above delineates the four different regressions between k-fold crossvalidations and spatial LOGO crossvalidations and puts them all into one dataframe, reg.summary.

Accuracy and Generalizability

# calculate and visualize MAE for each fold and across each regression 
data_group <- reg.summary %>%                                
  group_by(cvID, Regression) %>%
  dplyr::summarize(Mean_Error = mean(Prediction - countOverdoses, na.rm = T),
                   MAE = mean(abs(Mean_Error), na.rm = T), 
                   SD_MAE = mean(abs(Mean_Error), na.rm = T)) %>% 
  as.data.frame()%>%
  st_sf()%>%
  na.omit()

data_group %>%
  ggplot(aes(MAE)) + 
    geom_histogram(colour="black", fill = "#FDE725FF") +
    facet_wrap(~Regression) +  
    geom_vline(xintercept = 0) + 
    labs(title="Distribution of MAE", subtitle = "k-fold cross validation vs. LOGO-CV",
         x="Mean Absolute Error", y="Count") +
    plotTheme()

Figure 12. This is a multiple histogram for the distributions of Mean Absolute Error (MAE) for each type of regression. LOGO-CV regressions assume that the local spatial process from all neighborhoods generalizes to the hold-out. When the local spatial process is not taken into consideration - as in the case of “Just Risk Factors” regressions - some hold-outs have large MAEs, however these errors are remediated when “Spatial Process” features from the local Moran’s I test are included.

# summary table of mean and standard deviation for MAE by Regression 
table <-  st_drop_geometry(data_group) %>%
  group_by(Regression) %>%
  dplyr::summarize(Mean_MAE = round(mean(MAE), 2),
                   SD_MAE = round(sd(MAE), 2)) %>%
  kable() %>%
    kable_styling("striped", full_width = F) %>% 
    row_spec(2, color = "black", background = "#FDE725FF") %>%
    row_spec(4, color = "black", background = "#FDE725FF")
table 

Regression	Mean_MAE	SD_MAE
Random k-fold CV: Just Risk Factors	0.32	0.37
Random k-fold CV: Spatial Process	0.33	0.39
Spatial LOGO-CV: Just Risk Factors	0.44	0.48
Spatial LOGO-CV: Spatial Process	0.44	0.48

Figure 13. This is a summary table of MAE and standard deviation for MAE by each regression. Here, we see that what we said in Figure 13 was correct; that regressions with spatial features included have lower errors than those with just risk factors.

# visualize LOGO-CV errors spatially 
data_group %>%
  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") +
    mapTheme() + theme(legend.position="bottom")

Figure 14. These maps visualize the LOGO-CV errors spatially, showing where high errors occurred and vice versa. When comparing these with the hotspots designated in Figure 7, it is not surprising that the largest errors are in the hotspot locations. Again, we see in the map on the right that when spatial features are included in the regression along with the risk factors, MAE is lower.

Generalizability by Neighborhood Context

# load in and wrangle census data 
census_api_key("ea7bd3babfc0da6a4caa6d2536b66403a88929f2", overwrite = TRUE)
 
tracts18 <- 
  get_acs(geography = "tract", variables = c("B01003_001E","B02001_002E","B02001_003E",
                                             "B03002_003E","B02001_005E","B19013_001E",
                                             "B15002_015E","B15002_032E","B06012_002E"), 
          year=2018, state=04, county=013, geometry=T, output='wide') %>%
  st_transform('EPSG:3857') %>%
  rename(TotalPop = B01003_001E, 
         Whites = B02001_002E,
         FemaleBachelors = B15002_032E, 
         MaleBachelors = B15002_015E,
         MedHHInc = B19013_001E, 
         AfricanAmericans = B02001_003E,
         Latinos = B03002_003E,
         Asians = B02001_005E,
         TotalPoverty = B06012_002E) %>%
  dplyr::select(-NAME, -starts_with("B")) %>%
  mutate(pctWhite = ifelse(TotalPop > 0, Whites / TotalPop,0),
         pctBlack = ifelse(TotalPop > 0, AfricanAmericans / TotalPop,0),
         pctLatino = ifelse(TotalPop > 0, Latinos / TotalPop,0),
         pctAsian = ifelse(TotalPop > 0, Asians / TotalPop,0),
         pctBachelors = ifelse(TotalPop > 0, ((FemaleBachelors + MaleBachelors) / TotalPop),0),
         pctPoverty = ifelse(TotalPop > 0, TotalPoverty / TotalPop, 0),
         pctMinority = ifelse(TotalPop > 0, (AfricanAmericans + Latinos + Asians) / TotalPop,0),
         raceContext = ifelse(pctMinority > 0.5, "Majority Minorities", "Majority White"),
         incomeContext = ifelse(MedHHInc >= 50000, "Relatively-High Income", "Relatively-Low Income"))

 tracts18_area <- tracts18 %>%
  mutate(area = as.numeric(st_area(tracts18)),
         areaContext = ifelse(area >= 35054741, "Relatively-High Area", "Relatively-Low Area"),
         year = "2018") %>%
  dplyr::select(-Whites, -FemaleBachelors, -MaleBachelors, -TotalPoverty, -AfricanAmericans, -Asians, -Latinos) %>%
  na.omit() %>%
  .[data_group,]

finished18 <- subset(tracts18_area, GEOID != "04013941300") %>%
  na.omit()

real18 <- subset(finished18, GEOID != "04013010102") %>%
  na.omit() %>%
  .[data_group,]
  
# map race and income context 

grid.arrange(ncol=2,
ggplot() + 
  geom_sf(data = st_union(data_group)) +
  geom_sf(data = real18, aes(fill=raceContext)) + 
  labs(title = "Race Context",
       subtitle = "Mesa, AZ") +
  mapTheme(title_size = 20),

ggplot() +
  geom_sf(data = st_union(data_group)) +
  geom_sf(data = real18, aes(fill=incomeContext))+
  labs(title="Income Context", 
       subtitle="Mesa, AZ") +
  mapTheme(title_size = 20))

Figure 15. These maps visualize 2018 Census data by tract to delineate Mesa between majority white and majority non-white populations and relatively-high and relatively-low median household incomes. Here, majority non-white is any neighborhood that is less than 50% white and relatively-low income is any neighborhood whose median household income is less than $50,000. Comparing this map with the hotspot map in Figure 7, we see that West Mesa (the overdose hotspot) has a generally more white population and has a lower median household income than the rest of the city.

# summary table of error by racial context 
reg.summary %>% 
  filter(str_detect(Regression, "LOGO")) %>%
    st_centroid() %>%
    st_join(real18) %>%
    na.omit() %>%
      st_drop_geometry() %>%
      group_by(Regression, raceContext) %>%
      dplyr::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) %>%
        row_spec(2, color = "black", background = "#FDE725FF") 

Mean Error by neighborhood racial context

Regression	Majority Minorities	Majority White
Spatial LOGO-CV: Just Risk Factors	0.0290985	-0.2992854
Spatial LOGO-CV: Spatial Process	0.0188130	-0.2096411

Figure 16. This is a summary table of MAE with respect to neighborhood racial context. This is particularly useful when seeing how well our model generalizes across neighborhoods. What we can gain from this table is that our model generally under-predicts MAE in majority non-white neighborhoods and over-predicts in majority white neighborhoods; this goes hand-in-hand with our aforementioned comments on selection bias within the dependent variable, overdoses. The Spatial Process model reports lower errors for both racial contexts, which leads us to assume that our model generalizes well with respect to racial context.

# summary table of error by income context 
reg.summary %>% 
  filter(str_detect(Regression, "LOGO")) %>%
    st_centroid() %>%
    st_join(finished18) %>%
    na.omit() %>%
      st_drop_geometry() %>%
      group_by(Regression, incomeContext) %>%
      dplyr::summarize(mean.Error = mean(Error, na.rm = T)) %>%
      spread(incomeContext, mean.Error) %>%
      kable(caption = "Mean Error by neighborhood median household income context") %>%
        kable_styling("striped", full_width = F) %>%
        row_spec(2, color = "black", background = "#FDE725FF") 

Mean Error by neighborhood median household income context

Regression	Relatively-High Income	Relatively-Low Income
Spatial LOGO-CV: Just Risk Factors	0.0181854	-0.0434386
Spatial LOGO-CV: Spatial Process	-0.0048856	0.0025272

Figure 17. This is a summary table of MAE with respect to neighborhood median household income context. This is yet another socio-economic context to assess our model’s generalizability across neighborhoods. What we can gain from this table is that our model generally over-predicts MAE in relatively-low income neighborhoods and under-predicts in relatively-high income neighborhoods; this also goes directly with our aforementioned comments on selection bias within the dependent variable, overdoses, as we excepted higher counts of overdoses in the low income areas. Not only does the Spatial Process model report lower errors overall and a smaller difference between errors, but the errors are also very low, which suggests that our model may also generalize well with respect to median household income.

Traditional Kernel Density vs Risk Predictions

# demo of kernel width
overdose_ppp <- as.ppp(st_coordinates(overdoses.in.boundary), W = st_bbox(final_net))
overdose_KD.1000 <- spatstat.core::density.ppp(overdose_ppp, 1000)
overdose_KD.1500 <- spatstat.core::density.ppp(overdose_ppp, 1500)
overdose_KD.2000 <- spatstat.core::density.ppp(overdose_ppp, 2000)
overdose_KD.df <- rbind(
  mutate(data.frame(rasterToPoints(mask(raster(overdose_KD.1000), as(data_group, 'Spatial')))), Legend = "1000 Ft."),
  mutate(data.frame(rasterToPoints(mask(raster(overdose_KD.1500), as(data_group, 'Spatial')))), Legend = "1500 Ft."),
  mutate(data.frame(rasterToPoints(mask(raster(overdose_KD.2000), as(data_group, 'Spatial')))), Legend = "2000 Ft.")) 

overdose_KD.df$Legend <- factor(overdose_KD.df$Legend, levels = c("1000 Ft.", "1500 Ft.", "2000 Ft."))

ggplot(data=overdose_KD.df, aes(x=x, y=y)) +
  geom_raster(aes(fill=layer)) + 
  facet_wrap(~Legend) +
  coord_sf(crs=st_crs(final_net)) + 
  scale_fill_viridis(name="Density") +
  labs(title = "Kernel density with 3 different search radii") +
  mapTheme(title_size = 14)

Figure 18. These maps illustrate Kernel density at three different scales. Kernel density works by centering a curve over each theft point so that the curve is at its highest when it is directly above the points and the lowest at the circular search radius.

# convert ppp to dataframe and then sf layer and map kernel density with overlayed sample points 
as.data.frame(overdose_KD.1000) %>%
  st_as_sf(coords = c("x", "y"), crs = st_crs(final_net)) %>%
  aggregate(., final_net, mean) %>%
   ggplot() +
     geom_sf(aes(fill=value)) +
     geom_sf(data = sample_n(overdoses.in.boundary, 531), size = .5) +
     scale_fill_viridis(name = "Density") +
     labs(title = "Kernel density of 2018 overdoses") +
     mapTheme(title_size = 14)

Figure 19. This is a Kernel density map with a search radius of 1000 feet. Mean density is calculated and visualized with a sample of over 500 random points overlayed above it.

Comparing Kernel Density and Risk Factor Predictions for 2019 Thefts

# load in 2019 theft data 
overdoses19 <- read.csv("Data/Fire_and_Medical_Opioid_Overdose_Incidents.csv") %>%
  filter(Year == 2019) %>%
  na.omit() %>%
  st_as_sf(coords = c("Longitude","Latitude"), crs = 4326) %>%
  st_transform('EPSG:3857') %>%
  distinct() %>%
  .[data_group,]

# filter out points that don't lie within Mesa boundary 
overdoses19 <- overdoses19 %>% mutate(in_bounds = lengths(st_within(overdoses19,mesaBoundary)))
overdoses19.in.boundary <- overdoses19 %>%
  filter(in_bounds == 1)

The code block above reads in Mesa overdose data from 2019 and filters out the points that did not occur inside Mesa’s city boundaries.

# do kernel density on 2018 overdoses, make sf, and add count of overdoses in 2019 
overdose_ppp <- as.ppp(st_coordinates(overdoses.in.boundary), W = st_bbox(final_net))
overdose_KD <- spatstat.core::density.ppp(overdose_ppp, 1000)

overdose_KDE_sf <- as.data.frame(overdose_KD) %>%
  st_as_sf(coords = c("x", "y"), crs = st_crs(final_net)) %>%
  aggregate(., final_net, mean) %>% 
  mutate(label = "Kernel Density",
         Risk_Percent = ntile(value, 100),
         Risk_Category = case_when(Risk_Percent >= 90 ~ "90% to 100%",
                                   Risk_Percent >= 70 & Risk_Percent <= 89 ~ "70% to 89%",
                                   Risk_Percent >= 50 & Risk_Percent <= 69 ~ "50% to 69%",
                                   Risk_Percent >= 30 & Risk_Percent <= 49 ~ "30% to 49%",
                                   Risk_Percent >= 1 & Risk_Percent <= 29 ~ "1% to 29%")) %>%
  cbind(
    aggregate(
      dplyr::select(overdoses19.in.boundary) %>% mutate(overdoseCount = 1), ., sum) %>%
    mutate(overdoseCount = replace_na(overdoseCount, 0))) %>%
  dplyr::select(label, Risk_Category, overdoseCount)

This code block above computes Kernel density for 2018 data, converts the densities into 5 different risk categories, and joins the count of overdoses in 2019.

# do same thing for risk predictions 
overdose_risk_sf <-
  filter(reg.summary, Regression == "Spatial LOGO-CV: Spatial Process") %>%
  mutate(label = "Risk Predictions",
         Risk_Percent = ntile(Prediction, 100),
         Risk_Category = case_when(Risk_Percent >= 90 ~ "90% to 100%",
                                   Risk_Percent >= 70 & Risk_Percent <= 89 ~ "70% to 89%",
                                   Risk_Percent >= 50 & Risk_Percent <= 69 ~ "50% to 69%",
                                   Risk_Percent >= 30 & Risk_Percent <= 49 ~ "30% to 49%",
                                   Risk_Percent >= 1 & Risk_Percent <= 29 ~ "1% to 29%")) %>%
  cbind(
    aggregate(
      dplyr::select(overdoses19.in.boundary) %>% mutate(overdoseCount = 1), ., sum) %>%
      mutate(overdoseCount = replace_na(overdoseCount, 0))) %>%
  dplyr::select(label,Risk_Category, overdoseCount)

The code block above does the same as the last one, except this time on risk predictions instead of Kernel density.

Figure 20. These maps show the risk categories for both model types with a sample of close to 1000 2019 overdose points on top. Here, we see that our highest risk category seems to be on-target to locations with a relatively high density of overdose points, suggesting that our model may fit.

Figure 21. This bar graph illustrates the rate of 2019 overdose points in relation to both risk category and model type. Well-fit models will show risk predictions as having a higher share of 2019 overdose points in the highest risk categories than Kernel density. In this case, the risk prediction model is more than the Kernel density in three of the five risk categories, suggesting that this model has more value than the business-as-usual hotspot approach.

Conclusion

While running through our model, we ran into some results that suggested that our algorithm not be put into production. For one, the data that we used for both our dependent variable, overdoses, and some risk factors may have been influenced by selection bias. Since not all overdoses can be observed and since Mesa Fire and Medical Department personnel may be more willing to allocate themselves and overdose prevention resources more fervently in some communities than others, our model’s overdose outcomes may be influenced by selection bias. Additionally, risk factors -like drug crime and graffiti sightings -that are concentrated and not as spread out across Mesa- may introduce more selection bias into our model, since the increased counts in these areas may be derived from the factors being identified more fervently in some communities than others (Figures 3 and 4). After observing the correlation coefficients in the scatterplots between count of overdoses and each risk factor (Figure 9) and the correlation matrix (Figure 10), it looks as though some of the risk factors that we included in our model are not strongly correlated (i.e abandoned buildings and service requests for out street lights); this leads us to think that we could have used different risk factors in our model to produce more accurate predictions.Lastly, Figures 13 and 21 suggest that our model barely outperforms the random k-fold model and the business-as-usual, hotspot approach in predicting future overdoses. If we had a well-fit model, there would be a huge discrepancy in values between the rate_of_test_set_overdoses for our risk prediction model and the Kernel density model in the highest risk categories of Figure 21, however, the Kernel density approach does only a slightly worse job at this than our risk predictions. This is confirmed in Figure 13, where our spatial model’s MAE is only slightly less than that of the just-risk-factor model.

Despite the aforementioned reasons, we think that some parts of our algorithm could be put into production. While some of the correlation coefficients in Figures 9 and 10 were in-significant, some others - like graffiti and overdose.isSig.dist - had relatively high correlations. If we could run our model again, we would spend more time to identify risk factors with little selection bias and that we thought would better predict future overdoses. Additionally, our summary tables in Figures 13, 16 and 17 show how the Spatial Features regression had less error and latent risk in general and for both the racial and income contexts than the Just Risk Factors regression; this means that our model is more generalizable across socio-economic contexts, and that our Moran’s I spatial features generally improve our model and subsequent overdose prevention resource/service allocation. Finally, while there is only a slight discrepancy between the share of 2019 overdoses that are captured by the Kernel density approach and our risk prediction model (Figures 13 and 21), our risk prediction model still yields lower mean absolute errors and outperforms, confirming that our model is more valuable than the business-as-usual method and could be implemented into our app, MesaGuardian, to help Mesa’s Fire and Medical Department better allocate heroin prevention resources and services across the city.