library(sf)           
library(ggplot2)      
library(dplyr)        
library(tidyr)        
library(scales)       
library(RColorBrewer) 
library(units)
library(cowplot)
library(here)
library(classInt)
library(leaflet)
library(lubridate)
library(tidycensus)
library(tidyverse)
library(plotly)
library(ggiraph)
library(htmlwidgets)
library(mapview)
library(tigris)
library(spdep)
library(dbscan)
library(MASS)
library(fields)
library(sfhotspot)
library(tmap)
tmap_mode("view")

Abstract

Highly regarded fire departments like Seattle Fire (SFD) score generally high respective of national response standards for Emergency Medical Services (EMS) calls but experience measurable variability across call type categories and year to year. There have been many studies examining the effects of EMS response times on certain incident types and patient survivability, studies aiming to quantify the effects of urbanicity or ruralness of a given area on response times, as well as consider the geographical and transportation impacts on response times (Yoon et. Al. 2017, Lee et al 2021, Mell et al 2017, Swan and Baumstark 2021, Nguyen et al 2023). This report attempts to join this conversation and explore the impacts of Seattle’s urbanicity relative to the temporal and spatial variability of SFD’s response times from January 2021 to March 2024. This project suggests that the department is generally meeting their response time standards and there is slight clustering of like mean response times with neighboring census block groups. These analyses can help guide training, organization, and funding for departments like Seattle to ensure emergency responses are timely in keeping with industry standards as well as ensuring care is given to all communities in a department’s service area in an equitable manner.

Introduction

The city of Seattle, Washington is a modern densely populated city served by a large award-winning fire department.

*“The Seattle Fire Department became in 2022 the first fire department in the state’s history to receive the Protection Class 1 rating from the WSRB (Washington Rating and Surveying Bureau)... This places Seattle Fire in the top 1% of fire   departments in the nation. WSRB serves Washington residents as the nonprofit independent insurance rating agency for fire departments, fire districts and regional fire authorities within the state of Washington. With the new rating, the SFD will join approximately 460 other fire districts nationwide with a Class 1 rating.*” (Annual Report Seattle Fire Department, 2022).

Even a highly regarded fire department like SFD does not always meet their response time standards. Seattle Fire Department is held to a call processing standard that is guided by National Fire Protection Association Standard 1221(2019). According to their 2022 Annual Report, “SFD response time standard for the arrival of the first Basic Life Support unit (aid car, fire engine, ladder truck) is four minutes, 90% of the time. This is the time span between a unit being en route to on-scene”, in 2022 55% of the time EMS turnout time was within 60 seconds (Annual Report Seattle Fire Department, 2022). “SFD response time standard for the arrival of an Advanced Life Support unit (ALS) (medic unit with two firefighter/paramedics) is eight minutes, 90% of the time. This is the time span between a unit being en route to on-scene. In 2022, 82% of the time the first ALS unit arrival was 8 minutes (Annual Report Seattle Fire Department, 2022).

There have been many studies examining the effects of response times on certain incident types and patient survivability, studies aiming to quantify the effects of urbanicity or ruralness of a given area on response times, as well as consider the geographical and transportation impacts on response times (Yoon et. Al. 2017, Lee et al 2021, Mell et al 2017, Swan and Baumstark 2021, Nguyen et al 2023). My project aims to join this conversation and explore the impacts of Seattle’s urbanicity relative to the temporal and spatial variability of SFD’s response times from January 2021 to March 2024. This study seeks to discover what is the temporal and spatial variability of SFD’s response times in this dataset? Is there any clustering of call times in the study area?

Data

Data was requested from the SFD via a Public Data Request within the City of Seattle’s Open Data Portal. This data includes all 911 calls received by SFD from January 2021 through March 2024. The calls are described with the following attributes, standard incident ID number, date/time fields for when the call was received, assigned and the unit(s) arrived, incident type, address of the incident location, latitude and longitude for the incident location, and the unit(s) dispatched.

Census block group geometry was downloaded using the tigris package (Walker 2024). This geometry is provided by U.S Census Bureau TIGER line shapefile (U.S Census Bureau 2020).

Data Cleaning

For the 911 call dataset, limitations of the data include, of the total 434,832 calls included, 3,687 of them include null coordinate values with vague or “unknown” descriptions in the corresponding incident type fields. These are removed. Additionally, the latitude and longitude values were unusable in their original format, these were converted to valid geometries to plot. Next, all the date time fields were unusable as character strings, these are formatted into date time values using lubridate package (Grolemund 2011) to perform time calculations. A year field to easily compare call response values in the analysis. As mentioned, the incident descriptions are not standardized, and some are unclear and vague. To avoid response time outliers and to more accurately compare real response times with the standard, medical related calls were isolated into one dataset, and “other” call types were isolated into another. After these steps, the dataset containing 434,831 records was reduced to 219,862 for medical calls and 60,518 for other call types.

#format datetime fields using lubridate, create new year only field----
sfdcalls <- sfdcalls %>%
  mutate(timeassigned = mdy_hm(Time_First_Unit_Assigned),
         timearrive = mdy_hm(Time_First_Unit_Arrived),
         responsedate = mdy_hm(Response_Date),
         timepickup = mdy_hm(Time_PhonePickUp),
         responseyear = year(responsedate))

#format lat and long------
#found method https://gis.stackexchange.com/questions/429921/formatting-decimal-point-for-x-y-coordinates-using-r
#there are some random latitude values that have 9 digits. all other lat's have 8 
#and all longs have 9 so I will remove the 9th digit of the lats to standardize.

cleancalls <- sfdcalls %>% 
  filter(Latitude != 0 & Longitude != 0) %>%
  mutate(Latitude = ifelse(nchar(Latitude) == 9, floor(Latitude / 10), Latitude),
         lat = as.numeric(Latitude) * 0.000001, 
         long = as.numeric(Longitude) * -.000001,
         longitude = as.numeric(Longitude) * -.000001,
         latitude = as.numeric(Latitude) * 0.000001) %>%
  st_as_sf(coords = c("long", "lat"), crs = 4326)

Data standardization and Calculation Preparation

The incident type descriptions are not standardized, and some are unclear and vague. To avoid response time outliers and to more accurately compare response times with the industry standard, medical calls were isolated and stored in one dataset, and all other call types were isolated and stored in a separate dataset. The medical call incident values kept for analysis had the key words: “Encampment Aid”, “Hang-up”, “Aid Response”, “BC Aid Response”, “Yellow”, “Medical”, “Medic”, “Aid Response Freeway”. The “other” call type incident values kept for analysis had the key words: “Violence”, “Acuity”, “Water”, “Marine”, “Tunnel”, “Fire”, “Mutual”, “Rescue”, “Spill”. The following key words were not included in either dataset for analysis as these call types required transport, did not include providing aid, or required longer travel time to assist another department: “mutual aid”, “nurseline/amr”, and “trans to amr”. Here is a breakdown of the distinct incident type values for the medical call dataset and for other call types dataset.

#incident type wrangling

#store medical calls
medtypes <- c("Encampment Aid", "Hang-up","Aid Response", "BC Aid Response", "Yellow", "Medical", "Medic", "Aid Response Freeway")

# store other calls, except transports and other outliers
othertypes <- c("Violence", "Acuity", "Water", "Marine", "Tunnel", "Fire", "Mutual", "Rescue", "Spill")


#filter for all those stored values I want. https://www.statology.org/r-grepl-multiple-patterns/

medcleancalls <- cleancalls %>%
  filter(grepl(paste(medtypes, collapse='|'), Incident_Type)) %>%
  dplyr::select(-Latitude, -Longitude, -Time_First_Unit_Assigned, -Response_Date, -Time_First_Unit_Arrived, -Time_PhonePickUp)

othercleancalls <- cleancalls %>%
  filter(grepl(paste(othertypes, collapse='|'), Incident_Type)) %>%
  dplyr::select(-Latitude, -Longitude, -Time_First_Unit_Assigned, -Response_Date, -Time_First_Unit_Arrived, -Time_PhonePickUp)

#see all distinct incident types
distinctmedtype <- unique(medcleancalls$Incident_Type)

#see distinct indident types for other
distinctothertype <- unique(othercleancalls$Incident_Type)

distinctmedtype
##  [1] "Aid Response"                  "Medic Response"               
##  [3] "Aid Response Yellow"           "Medic Response, Overdose"     
##  [5] "Automatic Medical Alarm"       "Ladder Code Yellow"           
##  [7] "Medic Response, 7 per Rule"    "Mutual Aid, Medic"            
##  [9] "Medic Response, 6 per Rule"    "BC Medic Response, 6 per rule"
## [11] "Encampment Aid"                "Engine Code Yellow"           
## [13] "Aid Response Freeway"          "BC Aid Response"              
## [15] "Single Medic Unit"             "Medic Response, Overdose 14"  
## [17] "BC Medic Response"             "MVI Medic"                    
## [19] "MVI Freeway Medic"             "Medic Response Freeway"       
## [21] "Multiple Medic Resp 14 Per"    "Encampment Medic"
distinctothertype
##  [1] "Scenes Of Violence 7"        "Automatic Fire Alarm Resd"  
##  [3] "Fuel Spill"                  "Rubbish Fire"               
##  [5] "Auto Fire Alarm"             "Low Acuity Referral"        
##  [7] "Rescue Elevator"             "Low Acuity Response"        
##  [9] "Dumpster Fire"               "Brush Fire"                 
## [11] "Automatic Fire Alarm False"  "Mutual Aid, Medic"          
## [13] "Spill, Non-Hazmat"           "Rescue Extrication"         
## [15] "Encampment Fire"             "Car Fire"                   
## [17] "Scenes Of Violence Aid"      "Water Job Minor"            
## [19] "Water Rescue Response"       "Marine Service Response"    
## [21] "Car Fire Freeway"            "Fire in Building"           
## [23] "Bark Fire"                   "Rescue Lock In/Out"         
## [25] "Tunnel Standby"              "Tunnel Aid"                 
## [27] "Scenes Of Violence Major"    "Water Job Major"            
## [29] "Marine Fire On Shore"        "Tranformer Fire"            
## [31] "Chimney Fire"                "Fire Response Freeway"      
## [33] "Scenes Of Violence 14"       "Shed Fire"                  
## [35] "Fire In A Highrise"          "Mutual Aid, Ladder"         
## [37] "Hazardous Mat, Spill-Leak"   "Dumpster Fire W/Exp."       
## [39] "Rescue Standby"              "Vault Fire (Electrical)"    
## [41] "Mutual Aid, Engine"          "Rescue Rope"                
## [43] "Car Fire W/Exp."             "Hang-Up, Fire"              
## [45] "Garage Fire"                 "Brush Fire Freeway"         
## [47] "Help the Fire Fighter"       "Mutual Aid, Task Force"     
## [49] "Mutual Aid, Marine"          "Rescue Trench"              
## [51] "Mutual Aid, Aid"             "Water Rescue Standby"       
## [53] "Mutual Aid, Hazmat"          "Rescue Heavy Major"         
## [55] "Mutual Aid, Adv. Life"       "Vessel Sinking On Water"    
## [57] "Brush Fire W/Exp."           "Mutual Aid, Strike Eng."    
## [59] "Marine Fire On Water"        "Brush Fire Major"           
## [61] "Tunnel MVI"                  "Rescue Confined Space"      
## [63] "Tunnel Fire"                 "Water Rescue Response Major"
## [65] "Hazardous Material w/Fire"   "Scenes Of Violence MCI"     
## [67] "Tunnel North Ops Bldg"       "Tunnel Rescue Standby"      
## [69] "Train Derailment wFireHzmt"
# only medical calls
medcleancalls <- medcleancalls %>% 
  filter(timearrive != 0) %>%
  mutate(responsetime= (timearrive - timeassigned)/60, #date type response time for mapping
         num_responsetime = as.numeric(responsetime)) %>% #need this for statistics and for moran's
  filter(responsetime != 0)

medcleancalls <- medcleancalls %>% 
  filter(timearrive != 0) %>%
  mutate(responsetime= timearrive - timeassigned, #date type response time for mapping
         num_responsetime = as.numeric(responsetime)) %>% #need this for statistics and for moran's
  filter(responsetime != 0)

# all other calls, except transports and other outliers
othercleancalls <- othercleancalls %>% 
  filter(timearrive != 0) %>%
  mutate(responsetime= (timearrive - timeassigned)/60, #date type response time for mapping
         num_responsetime = as.numeric(responsetime)) %>% #need this for statistics and for moran's
  filter(responsetime != 0)

othercleancalls <- othercleancalls %>% 
  filter(timearrive != 0) %>%
  mutate(responsetime= timearrive - timeassigned, #date type response time for mapping
         num_responsetime = as.numeric(responsetime)) %>% #need this for statistics and for moran's
  filter(responsetime != 0)

The 911 calls as a point value dataset was used for some calculations, but in order to compute additional temporal comparisons, global and local spatial autocorrelation values using neighbors, the point data needed to be aggregated up to area geometries. All the calls that are located within each block group were aggregated into/grouped by that block group. The block group geometries for King County Washington are used for this purpose. This polygon dataset was projected into a suitable local projected coordinate system (EPSG 32610), some geometries were repaired and some were removed because they couldn’t be repaired (i.e. linestrings). Additionally, six of the very large easternmost block groups that did not participate in SFD calls for the duration of the dataset were removed.

#------------------------
#AGGREGATE RAW POINTS TO POLYGON AREAS

# used chapter 5 from lesson 4 : https://bookdown.org/mcwimberly/gdswr-book/vector-geospatial-data.html#practice-4
# do join to get geoid, then calc number of calls per block group (call density to scale means) and mean call time per block group
# join back to geometry to map


# error Error in wk_handle.wk_wkb(wkb, s2_geography_writer(oriented = oriented,  : 
#  Loop 0 is not valid: Edge 8 is degenerate (duplicate vertex)
# ran sf_use_s2(FALSE) to fix 

#spatial join, to inherit block group attributes based on spatial intersection/join

med_call_all_intersect <- st_join(medcleancalls_all, king_block_groups, left = FALSE)

other_call_intersect <- st_join(othercleancalls_filter, king_block_groups, left = FALSE)

#calculate mean response time per block group, all medical calls
med_call_all_drop <- st_drop_geometry(med_call_all_intersect)

med_call_all_time_mean <- med_call_all_drop %>%
  group_by(GEOID) %>%
  summarize(mean_responsetime = mean(num_responsetime), 
            n = n(),
            sd_responsetime = sd(num_responsetime)) %>%
  mutate(cv_responsetime = sd_responsetime / mean_responsetime)

#calculate mean response time per block group, all other calls
other_call_drop <- st_drop_geometry(other_call_intersect)

other_call_time_mean <- other_call_drop %>%
  group_by(GEOID) %>%
  summarize(mean_responsetime = mean(num_responsetime), 
            n = n(),
            sd_responsetime = sd(num_responsetime)) %>%
  mutate(cv_responsetime = sd_responsetime / mean_responsetime)


# join back to geometry, calc area for density calc, all medical calls
med_call_all_mean_map <- king_block_groups %>%
  left_join(med_call_all_time_mean, by = "GEOID") %>%
  replace(is.na(.), 0) %>%
  mutate(area = sf::st_area(king_block_groups),
         calldens = 10^6 * n / area) #density as measured 1 call/ sq km
#%>%  drop_units()

# join back to geometry, calc area for density calc, all other calls
other_call_mean_map <- king_block_groups %>%
  left_join(other_call_time_mean, by = "GEOID") %>%
  replace(is.na(.), 0) %>%
  mutate(area = sf::st_area(king_block_groups),
         calldens = 10^6 * n / area) #density as measured 1 call/ sq km
#%>% drop_units()

Methodology

This analysis explores the temporal distribution of call times for the medical and other call types, using basic statistical calculations of mean values, standard deviation, and coefficient of variation. In the second section, this report explores spatial distribution of call response times, first using Moran’s I for global spatial autocorrelation, and lastly, using Local Moran’s to find local clusters of call response times within the study area.

Temporal Distribution Analysis

Answering the question of what is the response time variability for the medical calls and other call types within the dataset, the raw point call data was plotted to see the distribution of the response times through the timeline of the dataset. These raw point values are compared against the mean response time of each block group and the entire dataset, the standard deviation of the mean response time of each block group and the entire dataset, the coefficient of variation of the mean response time of each block group and the entire dataset. Mean was chosen as there are enough calls in the dataset to provide an accurate value for the most common response times within each block group, standard deviation of response times was chosen to quantify the dispersion of response times in each block group. Coefficient of variation was calculated to quantify the ratio of standard deviation from the mean response time. As another measure of quantifying dispersion of calls across each block group, call density was calculated which is the number of calls taken (across the entire dataset) divided by the block group area.

# RAW POINTS RESPONSE TIME CALCS (USING MEDCLEANCALLS AND OTHERCLEANCALLS for points)

#dataset stats - medical
med_totalstats_raw <- medcleancalls %>%
  #filter(num_responsetime > 0 & num_responsetime < 140)
  summarize(mean_responsetime = mean(num_responsetime, na.rm = TRUE),
            sd_responsetime = sd(num_responsetime, na.rm=TRUE),
            cv_responsetime = sd_responsetime / mean_responsetime)

med_yearstats_raw <- medcleancalls %>%
  group_by(responseyear) %>%
  summarize(mean_responsetime = mean(num_responsetime, na.rm = TRUE),
            sd_responsetime = sd(num_responsetime, na.rm=TRUE),
            cv_responsetime = sd_responsetime / mean_responsetime)

#dataset stats - other calls
oth_totalstats_raw <- othercleancalls %>%
  summarize(mean_responsetime = mean(num_responsetime, na.rm = TRUE),
            sd_responsetime = sd(num_responsetime, na.rm=TRUE),
            cv_responsetime = sd_responsetime / mean_responsetime)

oth_yearstats_raw <- othercleancalls %>%
  group_by(responseyear) %>%
  summarize(mean_responsetime = mean(num_responsetime, na.rm = TRUE),
            sd_responsetime = sd(num_responsetime, na.rm=TRUE),
            cv_responsetime = sd_responsetime / mean_responsetime)

Spatial Distribution

Answering the question, where do these calls occur? Where do the call times differ? Where do longer response times occur? Where do shorter response times occur?

Global Spatial Autocorrelation using Moran’s I

Moran’s I is calculated for each dataset, to determine the presence of clustering mean response times across the entire study area as a measure of global spatial autocorrelation. Moran’s I considers the degree to which the mean response time of one block group is similar or dissimilar to its neighboring block group (Farbar 2013). This statistic is computed using the spdep package (Bivand 2013) and requires the calculation of a list of neighboring block groups, a weighted neighbors matrix based on relative proximity, then quantifying how similar each block group’s mean response time is in relation to its neighbors and averaging these assessments (Moraga 2024). The statistic ranges from –1 to 1, values closer to 1 indicate positive spatial autocorrelation or clustering, indicating that like mean values tend to be closer together. Moran’s I values closer to –1 indicate negative spatial autocorrelation or alternating values, suggesting higher degree of heterogeneity (highs next to lows and vise versa). Values closer to 0 indicate weak spatial autocorrelation or randomness, that is, absence of spatial pattern (Farbar 2013).

#------------------------------
# ALL MEDICAL CALLS
#calc neighborhoods
# https://www.paulamoraga.com/book-spatial/spatial-neighborhood-matrices.html
# https://mgimond.github.io/simple_moransI_example/#step_2:_assign_weights_to_the_neighbors
# https://mgimond.github.io/Spatial/spatial-autocorrelation-in-r.html

#prep for calc
# Filter the data to keep only POLYGON and MULTIPOLYGON geometries and convert GEOMETRYCOLLECTION which are invalid geoms
med_call_all_mean_map1 <- med_call_all_mean_map %>% 
  filter(!st_is(.,"MULTILINESTRING")) %>%
  filter(!st_is(.,"GEOMETRYCOLLECTION")) %>%
  #st_cast("POLYGON") %>%
  st_make_valid()

#remove block group geometries that have no calls to not affect statistic 
med_call_all_mean_map1 <- med_call_all_mean_map1 %>% dplyr::filter(n != 0)

#compute neighbors list
nball <- spdep::poly2nb(med_call_all_mean_map1, queen = TRUE)

#---
# calc spatial weights matrix, assign to neighbors
# each neighboring polygon will be assigned equal weight when computing the neighboring mean response time
# W is row standardised (sums over all links to n), what we need for moran's

nbweightsall <- nb2listw(nball, style="W", zero.policy = TRUE)

# calc average neighbor mean resopnsetime for each block group (spatially lagged values)
#allmeanresponsetime_lag <- lag.listw(nbweightsall, med_call_all_mean_map1$mean_responsetime)

# moran's i result is .815

medM <- spdep::moran(med_call_all_mean_map1$mean_responsetime, nbweightsall, length(nball),  zero.policy= TRUE, Szero(nbweightsall))
#spdep::moran.test(med_call_all_mean_map$mean_responsetime, nbweightsall, zero.policy= TRUE)
#------------------------------
# ALL "OTHER" NON MEDICAL CALLS
#calc neighborhoods
# https://www.paulamoraga.com/book-spatial/spatial-neighborhood-matrices.html
# https://mgimond.github.io/simple_moransI_example/#step_2:_assign_weights_to_the_neighbors

# Filter the data to keep only POLYGON and MULTIPOLYGON geometries and convert GEOMETRYCOLLECTION which are invalid geometries
other_call_mean_map1 <- other_call_mean_map %>% 
  filter(!st_is(.,"MULTILINESTRING")) %>%
  filter(!st_is(.,"GEOMETRYCOLLECTION")) %>%
  st_make_valid()

#remove block group geometries that have no calls to not affect statistic 
other_call_mean_map1 <- other_call_mean_map1 %>% dplyr::filter(n != 0)

nboth <- poly2nb(other_call_mean_map1, queen = TRUE)

#---
# calc spatial weights matrix, assign to neighbors
# each neighboring polygon will be assigned equal weight when computing the neighboring mean response time
# W is row standardised (sums over all links to n), what we need for moran's
nbweightsoth <- nb2listw(nboth, style="W", zero.policy = TRUE)

# calc average neighbor mean resopnsetime for each block group (spatially lagged values)
#othmeanresponsetime_lag <- lag.listw(nbweightsoth, other_call_mean_map1$mean_responsetime)

# moran's i result is .815
othM <- spdep::moran(other_call_mean_map1$mean_responsetime, nbweightsoth, length(nboth), zero.policy= TRUE, Szero(nbweightsoth))[1]
#spdep::moran.test(other_call_mean_map$mean_responsetime, nbweightsoth, zero.policy= TRUE)

Local spatial autocorrelation using Local Moran’s

“Local Indicators of Spatial Association (LISA) (Anselin 1995) are designed to provide an indication of the extent of significant spatial clustering of similar values around each observation” (Moraga 2024). This report uses a popular LISA, the local version of Moran’s I, to calculate which regions in the study area are similar or dissimilar (in mean response times) and their relationship via proximity (as neighbors). This statistic is also computed using the spdep package (Bivand 2013) where the neighbors weight matrix (calculated during global spatial autocorrelation), mean response time, and a two sided alternative hypothesis are used with the localmoran() function to create a local moran’s statistic for each block group, the expectation of the local moran’s statistic, the variance of the local moran’s statistic, a z-score and a p-value for an alternative hypothesis thats two sided (Moraga 2024). A two-sided alternative hypothesis tests no spatial autocorrelation vs positive or negative spatial autocorrelation (Moraga 2024). The classification of clusters is combined with (or compared with rather) an indication of significance, a pseudo p-value for each block group. So each block group’s scaled value for mean response time is compared with its spatially lagged value and its pseudo p-value is compared with a relative significance cut-off of 0.05, as recommend by Moraga 2024. Local Moran’s I identifies clusters of the following types: High-High, which are areas of high mean response time values with neighbors of high values; High-Low, areas of high mean response time values with neighbors of low values; Low-High, areas of low mean response time values with neighbors of high values; Low-Low, areas of low mean response time values with neighbors of low values (Moraga 2024).

Results

Temporal Distribution Analysis Results

Examining figures 1 and 2, most of calls each year fall well below 100 minutes, but there are few outlier calls over 100 minutes. There are also a few negative response times, this is likely due to SFD already being on scene and receiving an official aid call and dispatch after. The mean response time for medical calls across the entire dataset is 4.6 minutes, the standard deviation is 2.7 minutes and the coefficient of variation of the mean is .58, or the standard deviation is 58% of the mean. The mean response time for other call types across the entire dataset is 5.9 minutes, the standard deviation is 5.4 minutes, and the coefficient of variation is .92, or the standard deviation is 92% of the mean. Both datasets had a large amount of variability in response times. In the medical dataset, each year’s mean response time was around 4.5 minutes, compared to the other call type dataset where each year tended to have a mean around 6 minutes. Of note, in the other call dataset in 2023, the standard deviation and the mean response time were same, denoting that there is an extreme amount of variability in this year. To further examine the temporal variability, figures 3 and 4 are choropleth maps depicting the mean response times per block group for both datasets. The block groups with the largest mean response times tend to have the smallest call density, pointing to the large variability of response times within each block group.

# MAP MAP
max_width <- 50

scattermed <- medcleancalls %>% filter(num_responsetime != 0) %>%
  ggplot(., aes(x = as.factor(responseyear), y = num_responsetime)) +
  geom_point() +
  labs(#title = "Response Time Distribution per Year, All Medical Calls",
    x = "Year", y = "Response Time (minutes)") +
  theme_gray()

barmed <- ggplot(med_yearstats_raw, aes(x = responseyear)) + 
  geom_col(aes(y = mean_responsetime), fill = "yellow", position= "dodge") +
  geom_col(aes(y = sd_responsetime), fill = "green", position = "dodge") +
  geom_col(aes(y = cv_responsetime), fill = "lightblue", position = "dodge") + 
  labs(#title = "Comparison of Medical Response Time Variability by Year and Block Group",
       #subtitle = "CV of Mean in blue, SD of Mean in green, Mean in Yellow",
       caption = str_wrap("Figure 1. Medical Calls, 2021-2024. Comparing raw point response times with statistics derived from mean response times aggregated by census block groups. Block group mean response times derived from raw 911 point data. Standard deviation and coefficiet of variation derived from block group mean response time. All time in minutes. Census block group geometry from US Census Bureau 2020", width = max_width),
       x = "Year",
       y = str_wrap("CV of Mean in blue, SD of Mean in green, \nMean in Yellow", max_width))


#composite med 
commed <- plot_grid(scattermed, barmed, 
          labels = c(
            "A) Response Times Distribution per Year as Raw Points, All Medical Calls",
            "B) Response Time Variability by Year and Block Group, All Medical Calls",
            label_size = 8),
          ncol = 1, 
          label_x = 0,
          align = "hv",
          hjust = -0.15,
          vjust = .85)
commed

#--- raw stats and plots### 
scatteroth <- othercleancalls %>% filter(num_responsetime != 0) %>%
  ggplot(., aes(x = as.factor(responseyear), y = num_responsetime)) +
  geom_point() +
  labs(#title = "Response Time Distribution per Year, Other Call Types",
    x = "Year", y = "Response Time (minutes)") +
  theme_gray()

othmed <- ggplot(oth_yearstats_raw, aes(x = responseyear)) + 
  geom_col(aes(y = mean_responsetime), fill = "yellow", position= "dodge") +
  geom_col(aes(y = sd_responsetime), fill = "green", position = "dodge") +
  geom_col(aes(y = cv_responsetime), fill = "lightblue", position = "dodge") + 
  labs(#title = "Comparison of Medical Response Time Variability by Year and Block Group",
    #subtitle = "CV of Mean in blue, SD of Mean in green, Mean in Yellow",
    caption = str_wrap("Figure 2. Other Call Types, 2021-2024. Comparing raw point response times with statistics derived from mean response times aggregated by census block groups. Block group mean response times derived from raw 911 point data. Standard deviation and coefficiet of variation derived from block group mean response time. All time in minutes. Census block group geometry from US Census Bureau 2020.", width = max_width),
    x = "Year",
    y = str_wrap("CV of Mean in blue, SD of Mean in green, \nMean in Yellow", max_width))

#composite other
comoth <- plot_grid(scatteroth, othmed, 
                    labels = c(
                      "A) Response Times Distribution per Year as Raw Points, Other Call Types",
                      "B) Response Time Variability by Year and Block Group, Other Call Types",
                               label_size =8),
                    ncol = 1, 
                    label_x = 0,
                    align = "hv",
                    hjust = -0.15,
                    vjust = 1)
comoth

# map response time based on polygon areas + call density (using med_call_all_mean_map and other_call_mean_map for mean, sd, and cv by block group!)

#MEDICAL leaflet MAP

med_map <- med_call_all_mean_map %>% sf::st_as_sf() %>% 
  filter(mean_responsetime != 0) %>% st_transform(crs= 4326) # get rid of zeros and use geographic coord sys for leaflet

palette <- brewer.pal(9, "Greens")

# Classify mean response times using Jenks natural breaks
jenks_breaks <- classIntervals(med_map$mean_responsetime, n = length(palette), style = "jenks")
jenks_classes <- cut(med_map$mean_responsetime, breaks = jenks_breaks$brks, include.lowest = TRUE)

# Create a color palette with the classification
pal <- colorFactor(palette = palette, domain = jenks_classes)

# Create the leaflet map
leaflet(data = med_map) %>%
  addTiles() %>%
  addPolygons(
    fillColor = ~pal(jenks_classes),  # Use discrete color scheme based on Jenks breaks
    color = "black",
    weight = 1,
    fillOpacity = 1,
    label = ~mean_responsetime,
    popup = paste0(
      "GEOID: ", med_map$GEOID, "<br>",
      "Standard Deviation of Mean Response Time: ", med_map$sd_responsetime, "<br>",
      "Coefficient of Variation of Mean Response Time: ", med_map$cv_responsetime, "<br>",
      "Block Group Call Density: ", med_map$calldens
    )
  ) %>% 
  leaflet::addLegend('bottomright',
                     pal = pal,
                     values = jenks_classes,
                     opacity = 0.9, 
                     title = str_wrap("Mean Response Time (minutes) for Medical Calls", max_width))

Figure 3. Mean Response Times in minutes, for Medical Calls, 2021-2024. Block group mean response times and (pop-up) statistics derived from raw 911 point data. All time in minutes. Choropleth classified via Jenks Natural breaks in mean response times for each block group. Block group label is mean response time. Census block group geometry from US Census Bureau 2020.

#OTHER LEAFLET MAP
oth_map <- other_call_mean_map %>% sf::st_as_sf() %>% 
  filter(mean_responsetime != 0) %>% st_transform(crs= 4326)  # get rid of zeros and use geographic coord sys for leaflet

jenks_breaks2 <- classIntervals(oth_map$mean_responsetime, n = length(palette), style = "jenks")
jenks_classes2 <- cut(oth_map$mean_responsetime, breaks = jenks_breaks$brks, include.lowest = TRUE)

# Create a color palette with the classification
pal2 <- colorFactor(palette = palette, domain = jenks_classes2)

leaflet(data = oth_map) %>%
  addTiles() %>% 
  #addMarkers(clusterOptions = markerClusterOptions()) %>%
  addPolygons(fillColor = ~pal(jenks_classes2),  # Use continuous color scale for fill color
              color = "black", 
              #opacity = .9,
              weight = 1,
              fillOpacity = 1, 
              label = ~mean_responsetime, 
              popup = paste0("GEOID: ", oth_map$GEOID,"<br>",
                             "Standard Deviation of Mean Response time: ", oth_map$sd_responsetime, "<br>",
                             "Coefficient of Variation of Mean Response time: ", oth_map$cv_responsetime, "<br>",
                             "Block Group Call Density: ", oth_map$calldens)
  ) %>% 
  leaflet::addLegend('bottomright',
                     pal = pal2,
                     values = jenks_classes2,
                     opacity = 0.9, 
                     title = str_wrap("Mean Response Time (minutes) for Other Call Types",max_width))

Figure 4. Mean Response Times in minutes, for Other Call Types, 2021-2024. Block group mean response times and (pop-up) statistics derived from raw 911 point data. All time in minutes. Choropleth classified via Jenks Natural breaks in mean response times for each block group. Block group label is mean response time. Census block group geometry from US Census Bureau 2020.

Spatial Distribution Analysis Results

Global

The Moran’s I statistic for the medical call dataset is .34, suggesting a slight positive global spatial autocorrelation exists in the dataset for mean response times. As visible in the scatterplot in figure 5, this slight clustering in the medical dataset still has many outlier mean response times (all zeros were removed). The Moran’s I statistic for the other call type dataset is .29, suggesting also a slight positive global spatial autocorrelation exists in the dataset for mean response times. As seen in figure 6, and similar to the medical dataset, there are many outlier response times, but as stated earlier in the methodology, these were included to capture the true values present in the data paired with the large number of data (calls) available to still compute accurate mean values for analysis.

# scatter plot
# https://r-spatial.github.io/spdep/reference/moran.plot.html

#print(paste("Moran's I for Medical Calls: ", medM))

allmedicalscatter <- moran.plot(med_call_all_mean_map1$mean_responsetime, nbweightsall, zero.policy = TRUE,
           ylab= "Spatial Lag of Mean Response Times",
           xlab= "Mean Response Times by Block Group",
           main= "Moran's I Scatter Plot, All Medical Calls",
           pch=19)

Figure 5. Moran’s I Value for mean response times of Medical Calls by Census Block group, 2021-2024. All time in minutes. Suggests a minor positive spatial autocorrelation within the data.

# scatter plot
#print(paste("Moran's I for Other Call Types: ", othM))


otherscatter <- moran.plot(other_call_mean_map1$mean_responsetime, nbweightsoth, zero.policy = TRUE,
     ylab = "Spatial Lag of Mean Response Times",
     xlab = "Mean Response Times by Block Group",
     main = "Moran's I Scatter Plot, All Non-Medical Calls",
pch=19)

Figure 6. Moran’s I Value for mean response times of Other Call Types by Census Block group, 2021-2024. All time in minutes. Suggests a minor positive spatial autocorrelation within the data.

Local

Given the slight clustering of mean response times present in both datasets (slightly positive Moran’s I value), the few significant clusters of mean response times visible in figures 7 and 8 are expected. Note, each block group can be selected to view statistics for each geometry. Both datasets had block groups with no neighbors, and these were removed from the local Moran’s I analysis. Shown in figure 7, the medical dataset’s most prominent clusters of high mean response times near high mean response times (High-High) are located near White Center, just north of Burien and near the city of Shoreline where there is an additional cluster of low mean response times near other low mean response times. Additional Low-Lows are located in Mercer Island, downtown Seattle near the neighborhood of Capital Hill, and near the Ravenna neighborhood near Greenlake. As shown in figure 8, the other call type dataset’s most prominent clusters of High-Highs are located North near Shoreline with a few Low-High and two Low-Low neighbors. Suggesting a similar distribution of medical and other call type response times. On Mercer Island there are block groups with the same Low-Low classification in both datasets, but differing neighbors. West and (South) West Seattle’s cluster classification differs greatly in both datasets, however the clusters of High-Highs in the White Center area of South Seattle are similar. Additionally, while different in their specific block group distribution, both datasets have contrasting clusters of High-Highs near Low-Lows in the Ravenna neighborhood just south of Greenlake.

medlocal_map <- tm_shape(med_call_all_mean_map2) + tm_fill(col = "quadrant", title = "Clusters",
                        breaks = c(1, 2, 3, 4, 5, 6),
                        palette =  c("red", "blue", "lightpink", "skyblue2", "white"),
                        labels = c("High-High", "Low-Low", "High-Low",
                                   "Low-High", "Non-significant"),
                        popup.vars = c("GEOID" = "GEOID", "Mean Call Response Time" = "mean_responsetime", "Standard Deviation of Mean Response time" = "sd_responsetime", "Coefficient of Variation of Mean Response time" = "cv_responsetime", "Block Group Call Density" = "calldens" )) +
  tm_legend(text.size = 1)  + tm_borders(alpha = 0.5) +
  tm_layout(frame = FALSE,  title = "Clustering of Medical Response Times, Using Local Moran's")  +
  tm_layout(legend.outside = TRUE) 
medlocal_map

Figure 7. Local Moran’s Cluster Analysis on Mean Response times of Medical Calls by Census Block group, 2021-2024.(Pop-up) Block group mean response times and statistics derived from raw 911 point data. All time in minutes. Block group label is state id. Census block group geometry from US Census Bureau 2020.

othlocal_map <- tm_shape(other_call_mean_map2) + tm_fill(col = "quadrant", title = "Clusters",
                                           breaks = c(1, 2, 3, 4, 5, 6),
                                           palette =  c("red", "blue", "lightpink", "skyblue2", "white"),
                                           labels = c("High-High", "Low-Low", "High-Low",
                                                      "Low-High", "Non-significant"),
                                           popup.vars = c("GEOID" = "GEOID", "Mean Call Response Time" = "mean_responsetime", "Standard Deviation of Mean Response time" = "sd_responsetime", "Coefficient of Variation of Mean Response time" = "cv_responsetime", "Block Group Call Density" = "calldens" )) +
  tm_legend(text.size = 1)  + tm_borders(alpha = 0.5) +
  tm_layout(frame = FALSE,  title = "Clustering of Other Call Types Response Times Using Local Moran's")  +
  tm_layout(legend.outside = FALSE, legend.position = c('right', 'bottom')) 

othlocal_map

Figure 8. Local Moran’s Cluster Analysis on Mean Response times of Other Call Types by Census Block group, 2021-2024.(Pop-up) Block group mean response times and statistics derived from raw 911 point data. All time in minutes. Block group label is state id. Census block group geometry from US Census Bureau 2020.

Conclusion

911 calls taken by and responded by SFD from January 2021 through March 2024 containing 219,862 medical calls and 60,518 other call types were analyzed in this report. Generally, there is a lot of variability of call response time within the medical dataset and within the other call type dataset, the standard deviation of the mean for the entire medical dataset was 2.7 minutes, and is 5.4 for other call type dataset. The other call types dataset has higher mean values for each block group and across the dataset, 4.6 minutes for medical VS 5.9 minutes for other call types. In contrast to the distribution of mean response times for both datasets generated in the temporal analysis section, where similar high or low response times appeared to be located in close proximity, both datsets only had a slightly positive Moran’s I value, suggesting the datasets contained slight clustering of mean response times. Localized clusters generated as a result of computing Local Moran’s I for block group, resulted in a few significant clusters of mean response times. Notably both datasets had similar distribution of clustering classes in some areas, suggesting some consistency in response times across call types relative to Seattle neighborhoods.

Analysis was not conducted to measure the percentage of response time achieved relative to call volume (i.e percentage of the time SFD responded within their standards), but the medical dataset’s mean response times (~ four minutes each year) are within SFD’s response time standards as outlined earlier in this report, which is four minutes (90% of the time) for basic life support and eight minutes (90% of the time). Note that there the medical dataset contains both basic and advanced life support calls.

Additional analysis should include sensitivity for neighbors weight matrix chosen for both global and local spatial autocorrelation processes as well as randomness tests for both statistics. Also, temporal and spatial autocorrelation can be conducted on the raw call point data and compared with the call data aggregated to block group geometries (as done in this report). Further specification and standardization can be computed to delineate further between call types in each dataset.

This report analyses the temporal and spatial distribution of SFD’s 911 response times and suggests that the department is generally meeting their response time standard.

This report attempts to join this conversation and explore the impacts of Seattle’s urbanicity relative to the temporal and spatial variability of SFD’s response times from January 2021 to March 2024. This project seeks to provide greater insight into SFD’s current operations and possibly direct future training and funding, and ultimately suggests that the department is generally meeting their response time standard. These analyses can help guide training, organization, and funding for departments like Seattle to ensure emergency responses are timely in keeping with industry standards as well as ensuring care is given to all communities in a department’s service area in an equitable manner.

References

Abdishakur. (2019, November 5). What is Exploratory Spatial Data Analysis (ESDA)?
Medium. https://towardsdatascience.com/what-is-exploratory-spatial-data- analysis-esda-335da79026ee

Annual Report Seattle Fire Department. (n.d.). Retrieved March 24, 2024, from https://www.seattle.gov/documents/Departments/Fire/About/2022_Annual_Report_Web.pdf

Bureau, U.S. Census. (2020, January 1). Tiger/line shapefiles. Census.gov. https://www.census.gov/geographies/mapping-files/2020/geo/tiger-line-file.html

Cho, J., You, M., & Yoon, Y. (2017). Characterizing the influence of transportation infrastructure on Emergency Medical Services (EMS) in urban area—A case study of Seoul, South Korea. PLOS ONE, 12(8), e0183241. https://doi.org/10.1371/journal.pone.0183241

Emergency Medical Services Division - 2023 Annual Report to the King County Council, Public Health, Seattle & King County, Sept. 2023, https://cdn.kingcounty.gov/-/media/king-county/depts/dph/documents/health-safety/health-programs-services/emergency-medical-services/reports/2023-ems-annual-report.pdf?rev=114e9f3e1ff843dfb7374ffdd895b5eb&hash=A852B5EC7A2CF5D7EEC587CC1F3C154B

Farbar, S. (2013, June 11). GEOG 3020 Geographical Analysis (3020). Www.youtube.com; University of Utah. https://www.youtube.com/watch?v=aPYAiohGo74&list=PLas_mS2V1XIRsJZ9yACVeg8cvZucNnhvI&index=152

Gimond, Manuel. Chapter 13 Spatial Autocorrelation | Intro to GIS and Spatial Analysis. Mgimond.github.io, 15 Dec. 2023, mgimond.github.io/Spatial/spatial-autocorrelation.html#global-morans-i. Accessed 16 Apr. 2024.

Gimond, Manuel. “A Basic Introduction to Moran’s I Analysis in R.” Mgimond.github.io, 2019, mgimond.github.io/simple_moransI_example/#moran%E2%80%99s_i_analysis. Accessed 21 Apr. 2024.

Grolemund G, Wickham H (2011). “Dates and Times Made Easy with lubridate.” Journal of Statistical Software, 40(3), 1–25. https://www.jstatsoft.org/v40/i03/.

Moraga, Paula. Chapter 7 Spatial Neighborhood Matrices | Spatial Statistics for Data Science: Theory and Practice with R. Www.paulamoraga.com, www.paulamoraga.com/book-spatial/spatial-neighborhood-matrices.html#neighbors-based-on-k-nearest-neighbors. Accessed 21 Apr. 2024.

Pebesma E, Bivand R (2023). Spatial Data Science With Applications in R. Chapman & Hall. https://r-spatial.org/book/.

Seattle geodata. Seattle GeoData. (n.d.). https://data-seattlecitygis.opendata.arcgis.com/

Walker, K. (2024, January 24). Load Census Tiger/Line shapefiles [R package Tigris version 2.1]. The Comprehensive R Archive Network. https://cran.r-project.org/web/packages/tigris/index.html

Wimberly, Michael C. Chapter 5 Vector Geospatial Data | Geographic Data Science with R: Visualizing and Analyzing Environmental Change. Bookdown.org, Chapman & Hall/CRC, 19 Feb. 2023, bookdown.org/mcwimberly/gdswr-book/vector-geospatial-data.html. Accessed 21 Apr. 2024.

911 Calls, Seattle Fire Department, January 2021 – March 2024.