Introduction

Determining geospatial accessibility to a healthcare provider is more complicated than just determining if an individual living in a specific area could reach a provider within so many minutes by driving or by using public transportation methods. That only measures the existence of a healthcare provider within a certain amount of distance from the individual, but doesn’t measure if the provider is accessible to that individual. One big reason that the provider may not be accessible to that individual is the local demand for that provider, which is especially important with specialties that may bring patients from all over the region to that provider. Luo and Qi developed the Enhanced Two-Step Floating Catchment Area (E2SFCA) method to provide a measurement that intuitively makes sense for determining geospatial access to heathcare providers.

Using E2SFCA to measure spatial access begins with a collection of points that represent where the population seeking that healthcare lives. These should be central points in smaller geographic units (e.g. the centroid of a census tract or of a ZIP code–see notes below on using ZIP codes), and the analysis should include a point for each of those geographic units in that region. Isochrones are then built from each of these points with multiple subzones to represent different ranges of time that an individual may be willing to commute to that provider. For example, the graphic below depicts a catchment area with two sub-zones built around a house. The orange subzone could represent a 0-to-15 minute drive from that house, and the green subzone could represent a 15-to-30 minute drive from that house:


isochrone example


The locations of providers are then plotted over those catchment areas and a matrix identifying which subzone each provider fits into for each starting location is generated. That matrix is then converted into weights, which can be calculated through a variety of methods, where the higher weights are assigned to shorter commuting combinations, and lower weights are assigned to longer commuting combinations. For example, a provider that is the smallest subzone for a location in the analysis could have a weight of “1”, and a provider that is outside of the entire catchment area for that location could have a weight of “0”.

This matrix of weights is then combined with supply measurements for the providers, which could be the number of providers in that location, the number of prescriptions they can fill for that treatment, the number of patients they see in a specified period of time, etc., and demand measurements for those geographic units–which could be total population, adult population, population with a specific health condition, etc. All of that information is then plugged into a formula that returns a value for each geographic unit (called Spatial Access Index (SPAI)), which is difficult to interpret on its own, but the ranking of that value matters for identifying which geographic units have better access to that type of healthcare provider or service than other areas in their region. Standardizing that value by dividing it by the regional average (called Spatial Access Ratio (SPAR)) makes it easier to identify which areas have higher access than the regional average and which areas have lower access than the regional average.

The rest of this document provides a step-by-step example for using E2SFCA to measure spatial access in every ZIP code in Allegheny County, Pennsylvania to a chain of hospitals (UPMC) with Pittsburgh addresses–Pittsburgh is the county seat of Allegheny County and is located near the geographic middle of the county. The results are plotted in an interactive map at the end of this tutorial.

Notes about using ZIP codes for spatial analysis: ZIP codes are convenient smaller geographic units to use because address information usually includes a ZIP code; however, ZIP codes are a collection of addresses that do not represent a physical area, which means that there is both no real central location for a ZIP code and that the collection of addresses included in a ZIP code do not necessarily have a lot in common. Here are a few articles about the issues with using ZIP codes for this type of analysis:

Install and load libraries

This step uses the pacman package to install and load all of the packages for this tutorial.

# Use pacman package to load packages used for exercise
library(pacman)

p_load(tidyverse,   # Variety of packages for data manipulation
       sf,          # Variety of functions for manipulating spatial data
       kableExtra,  # For creating tables
       tidycensus,  # Package for downloading census data
       RCurl,       # For composing general HTTP requests
       rjson,       # Converts R objects to JSON objects and vice-versa
       leaflet,     # For creating interactive maps
       hereR,       # Package for accessing HERE isoline API
       rgdal,       # Geospatial bindings (used for writing and reading .shp files)
       sp,          # Spatial data package
       profvis)     # Includes pause() function to slow down API calls

1. Create catchment areas

The first step for using E2SFCA to determine spatial accessibility is to create catchment areas for all of the geographic units. This example uses all of the ZIP codes in Allegheny County, Pennsylvania. ZIP codes do not have a specific shape, so the Census Bureau’s estimated equivalents, called ZIP Code Tabulation Areas (ZCTA), will be used to represent ZIP codes and to help identify a central location for each ZIP code.

Find ZIP code central points

Set Census API key

ZCTA shapefiles can be downloaded directly from the Census Bureau using the code below. Before using that code, however, a Census API key must be obtained and installed in the R startup file.

A Census API key can be requested here.

# Set Census API key
census_api_key("[Insert API Key Here]",
               install = TRUE)

# Re-reads startup file to load Census API key
readRenviron("~/.Renviron")

Download ZIP code data

The first portion of the code below identifies each ZIP code with a presence in Allegheny County, using the 2019 4th quarter ZIP code to county crosswalk provided by the U.S. Department of Housing and Urban Development (HUD).

Nationwide ZCTA data with 2019 population estimates (the most recent release) and shapefiles are then downloaded from the Census Bureau and filtered to only include the Allegheny County ZIP codes identified in the HUD data.

Lastly, the HUD data includes proportion of residents with a ZIP code that live in that county (for ZIP codes in only one county, that is 100%), which is multipled by the population data downloaded from the Census Bureau to estimate the Allegheny County population for that ZIP code.

# List of ZIP codes in Allegheny county
# Note: some of these ZIP codes exist in multiple counties
allegheny_zips_list <- read.csv("./data/ZIP_COUNTY_122019.csv") %>%
  filter(COUNTY == 42003) %>%
  select(ZIP, RES_RATIO)

# Download Allegheny ZIP Codes
allegheny_zips <- get_acs(geography = "zcta",
                          variables = "DP05_0001E",
                          state = "PA",
                          geometry = TRUE) %>%
  filter(GEOID %in% allegheny_zips_list$ZIP) %>%
  mutate(ZIP = as.numeric(GEOID)) %>%
  rename(Population = estimate) %>%
  select(ZIP, Population)

# Estimate population using HUD residential percent
populations <- left_join(allegheny_zips,
                         allegheny_zips_list,
                         by = "ZIP") %>%
  mutate(Allegheny_Population = round(Population * RES_RATIO, 0)) %>%
  select(ZIP, Allegheny_Population) %>%
  
  # Keep only ZIP codes with a population in Allegheny County
  filter(Allegheny_Population > 0)

# Create a ZIP code dataframe
zips <- as.data.frame(populations$ZIP)
colnames(zips) <- "ZIP"

Geocode with Bing Maps

A central point for each ZIP code is then found using the Bing Maps API. As explained in greater detail here, plugging a ZIP code into the Bing Maps geocoding API returns the latitude and longitude of a street address near the center of that ZIP code. This is more useful than the exact geographic center of the ZIP code, which may not be on or even near a street.

Before running this step, first create a Bing Maps Developer account and then generate an API Key.

# Bing Maps API Key
BingMapsKey <- "[insert your API key here]"

# Loop through each ZIP code and add geographic coordinates
for(i in 1:nrow(zips)){
  
  url <- URLencode(paste0("http://dev.virtualearth.net/REST/v1/Locations?q=",
                          zips$ZIP[i],"&maxResults=1&key=",BingMapsKey))
  json <- fromJSON(getURL(url), simplify = FALSE)
  
  if (json$resourceSets[[1]]$estimatedTotal > 0) {
    lat <- json$resourceSets[[1]]$resources[[1]]$point$coordinates[[1]]
    lon <- json$resourceSets[[1]]$resources[[1]]$point$coordinates[[2]]
  }  else {
    lat <- lon <- NA
  }
  
  zips$Latitude[i] <- lat
  zips$Longitude[i] <- lon
}

# Plot map with ZIP code central points
leaflet() %>%
  addProviderTiles("CartoDB.Positron", group="Greyscale") %>% 
  addCircleMarkers(lng = zips$Longitude,
                   lat = zips$Latitude,
                   radius = 4,
                   fillOpacity = .75,
                   popup = paste0("ZIP code ",zips$ZIP))

Create isochrones

The steps below show how to create isochrones for all of the ZIP code central points geocoded above. This step saves each isochrone–with catchment areas for 0-10 minutes and 10-20 minutes–as .shp shapefiles in a folder, so that they can be reloaded one-by-one while filling out the catchment matrix below.

Create sf object for ZIP codes

The API used below to create isochrones works best with simple features. The step below converts the ZIP code dataframe to a simple feature object.

# Convert ZIP codes to SF object

zips.sf <- st_as_sf(x = zips,
                    coords = c("Longitude", "Latitude"),
                    crs = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")

Create isochrones using HERE API

The step below builds the catchment area isochrones for all of the ZIP code central points geocoded above from the simple feature object. Each isochrone is saved as a.shp file in the folder created in this step.

The first portion of this step sets the user’s API key for HERE’s isoline tool and then proceeds to create and save each catchment area shapefile. For more background on this step, see this guide on building isochrones with HERE.

# Set HERE key
here_key <- "[insert your API key here]"
set_key(here_key)

# Create directory to store isochrone shapefiles
dir.create("./data/Isochrones/", showWarnings = FALSE)

# Loop to create isochrones for each ZIP code
for(i in 1:nrow(zips.sf)){
  
  isochrone <- isoline(
    poi = zips.sf[i,],
    range = seq(10, 20, 10) * 60,
    range_type = "time",
    datetime <- as.POSIXct(paste0(Sys.Date()," 10:00"))
    ) %>%
    mutate(name = paste0((range - 600) / 60," to ", range / 60, " mins"))
  
  # Create spatial polygons data frame
  isochrone.sp <- as(isochrone, "Spatial")
  
  # Save spatial polygons data frame as a shapefile in a folder called "Isochrones"
  dsn <- paste0("./data/Isochrones")
  layer <- zips.sf$ZIP[i]
  writeOGR(obj = isochrone.sp,
           dsn = dsn,
           layer = layer,
           driver = "ESRI Shapefile")
  
  # Add a 1 second pause to prevent overload
  pause(1)
  
}

2. Create catchment matrix

The catchment areas created above can be used to locate roughly how far away each provider is from each ZIP code central point. They can also be used to determine roughly how far away each ZIP code central point is from each provider. Using this information in both directions is important for E2SFCA, measuring both which providers are accessible within a certain amount of drive time from a ZIP code central point and how large the demand on that provider is–based on which ZIP codes could access that provider.

Below is the isocrhone example from above but now there are also 3 hospitals added to the image (labeled “A”, “B”, and “C”). Hospital A is located within the orange subzone, which means it is within a 0-15 minute drive from that house. Hospital B is located within the green subzone, so it will take 15-30 minutes to drive from that house to that hospital. Lastly, hospital C is outside of the catchment area entirely, so it takes more than 30 minutes to drive there from that house.


isochrone example


The hospitals above aren’t just accessible by that single house, which impacts the amount of potential patients that utilize their services. Here is the same example but with a second house and its catchment area with the same subzones added. The original house above is now labeled as house “1”, and a new house is labeled as house “2”. The subzones for house 2 are 0-15 minutes in purple and 15-30 minutes in yellow.


isochrone example


Hospitals A and B are outside of the catchment area for house 2, but hospital C is within the 15-30 minute subzone for house 2.

Combining all of this information, we can fill out the following matrix of information:

Example Catchment Matrix
House 1 House 2
Hospital A 0-15 minutes > 30 minutes
Hospital B 15-30 minutes > 30 minutes
Hospital C > 30 minutes 15-30 minutes

The values in the matrix can be interpreted in two ways: which hospitals are accessible (and in which subzone) for each house, and which houses (and in which subzone) are in driving range for each hospital. For example, hospital B is located within a 15-30 minute drive from house 1, and house 1 is located within a 15-30 minute drive from hospital B. Together, this data shows which hospitals are accessible to which houses, and how many houses are able to access those hospitals.

Geocode providers

The next step for developing a catchment matrix to store this information is to make sure that the addresses of all of the healthcare providers are geocoded and can be plotted over the catchment areas previously created.

# Load UPMC hospitals in Pittsburgh
upmc_hospitals <- read.csv("./data/UPMC Locations.csV") %>%
  
  # Filter only those with "Pittsburgh" in the address
  filter(grepl("Pittsburgh", Address))

# Loop through each hospital address and add geospatial coordinates
for(i in 1:nrow(upmc_hospitals)){
  
  url <- URLencode(paste0("http://dev.virtualearth.net/REST/v1/Locations?q=",
                          upmc_hospitals$Address[i],"&maxResults=1&key=",BingMapsKey))
  json <- fromJSON(getURL(url), simplify = FALSE)
  
  if (json$resourceSets[[1]]$estimatedTotal > 0) {
    lat <- json$resourceSets[[1]]$resources[[1]]$point$coordinates[[1]]
    lon <- json$resourceSets[[1]]$resources[[1]]$point$coordinates[[2]]
  }  else {
    lat <- lon <- NA
  }
  
  upmc_hospitals$Latitude[i] <- lat
  upmc_hospitals$Longitude[i] <- lon
}

# Convert ZIP codes to SF object
upmc_hospitals.sf <- st_as_sf(x = upmc_hospitals,
                              coords = c("Longitude", "Latitude"),
                              crs = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")

# Convert to spatial points
upmc_hospitals.sp <- as(upmc_hospitals.sf, "Spatial")

Populate catchment matrix

Geocoded provider locations can be placed over a catchment area and retrieve information from the catchment area in that exact location uver the over() function from the sp package. The steps below create a catchment matrix with one column for every ZIP code catchment area and one row for every hospital geocoded above. The for() loop loads a single catchment area shapefile, plots each hospital location over it, and then records the subzone name in the appropriate cell of the catchment matrix. In this case, a cell will read “0-10 minutes”, “10-20 minutes” or “NA” for hospitals that are not located in a subzone for that ZIP code catchment area.

# Create list of isochrones
isochrone_list <- "./data/Isochrones/"

shape.extension <- list.files(path = isochrone_list,
                              pattern = ".shp",
                              full.names = TRUE)

shape.names <- list.files(path = isochrone_list,
                          pattern = ".shp",
                          full.names = FALSE) %>%
  str_remove(".shp")

# Build and label matrix
catchment_matrix <- matrix(ncol = length(shape.names), 
                           nrow = nrow(upmc_hospitals.sp))

colnames(catchment_matrix) <- shape.names
rownames(catchment_matrix) <- upmc_hospitals.sp$Location

# Fill out matrix
for(i in 1:length(shape.names)){
  
  # Loads isochrone shapefile
  isochrone <- readOGR(shape.extension[i])
  
  # Matches projection with hospital points projection
  isochrone@proj4string@projargs <- upmc_hospitals.sp@proj4string@projargs
  
  # Plots hospital points over catchment areas and takes name of subzone
  catchment_matrix[,i] <- over(upmc_hospitals.sp,
                               isochrone)$name
}

3. Create weighted matrix

The catchment matrix identifies the distance between each healthcare provider and each ZIP code central location, but that information itself doesn’t provide the values needed for the calculations for this analysis. Converting that information into weights–with higher values for closer combinations–to represent how much spatial influence each combination has on each other provides the values used for the rest of this analysis.

Create distance decay function

One way to create these weights is to use a distance decay function with predefined parameters to calculate weights for each subzone. These weights represent how much influence each location has on each other, and emphasize that closer combinations have a greater influence on each other than combinations that are further away from each other. The example below creates a Gaussian function (called an impedance function) that assigns a weight of 1 to combinations that are within a 10 minute drive of each other, a value of 0.13 to combinations between a 10 and 20 minute drive of each other, and a weight of 0 for all other combinations. These weights are calculated using the Gaussian impedance function.

# Gaussian impedance function
Gweight <- function(mintime, subzonetime, maxtime){
  beta <- -((maxtime-mintime)^2)/log(0.01)
  weight <- exp(-(subzonetime-mintime)^2/beta)
  weight <- round(weight, 2)
  return(weight)
}

# Notes for using impedance function: 
# 1. subzonetime = half of time between subzones (i.e. for 10/20 they are 5, 15)
# 2. mintime = subzonetime for smallest subzone (i.e. for 10/20 it is 5)
# 3. maxtime = maximum value of entire catchment area (i.e. for 10/20 it is 20)

# Create weights for catchment matrix
weight_NA <- 0
weight_10 <- Gweight(mintime = 5, subzonetime =  5, maxtime = 20)
weight_20 <- Gweight(mintime = 5, subzonetime = 15, maxtime = 20)

Fill out weight matrix

This step creates a new matrix, a weight matrix, which replaces the catchment matrix values with the weights developed above.

# Define names of subzones in catchment matrix 
sub1 <- "0 to 10 mins"
sub2 <- "10 to 20 mins"

# Fill out weight matrix values
weight_matrix <- data.matrix(catchment_matrix, rownames.force = NA)
weight_matrix[is.na(weight_matrix)]   <- weight_NA
weight_matrix[weight_matrix == sub1]  <- weight_10
weight_matrix[weight_matrix == sub2]  <- weight_20

# Make sure all values in matrix are numeric
class(weight_matrix) <- "numeric"

4. Demand and supply measurements

Demand vector (population)

To prepare for the final calculations in this analysis, a vector with populations for each of the ZIP codes needs to be prepared to be in the same order as the rows of the catchment and weight matrix (from left to right). It also has to be oriented so that it can multiplied with the weight matrix. If the weight matrix is x rows by y columns, where x is the number of healthcare providers, and y is the number of geospatial units, then the population vector must contain the population for all of those y geospatial units. The matrix calculations will be further explained below.

# The order of zip codes must match the order in your matrix
population_vector <- t(populations$Allegheny_Population)
colnames(population_vector) <- zips$ZIP
rownames(population_vector) <- "Population"

Supply-to-demand ratio

Multiplying the population vector on the left by the weight matrix on the right will provide a vector with information for each of the x providers. This will then be modified to capture the supply of services that each healthcare provider can provide.

In this example, the population_vector is 1 row by y columns, and the weight_matrix is x rows by y columns, so the setup for the multiplication is the population vector multiplied by the product of the weight matrix.

The results up to this point are the Allegheny County population within the catchment area of each hospital, weighted by distance from that hospital. The result is a 1 row by x columns vector.

# Create provider vector (i.e. weighted population for each provider)
provider_vector <- population_vector %*% t(weight_matrix)

# Display provider vector values
provider_vector %>%
  kbl(caption = "Provider Vector", 
      align = "r",
      format.args = list(big.mark = ",",
                         scientific = FALSE)) %>%
  kable_minimal("hover") %>%
  row_spec(0, font_size = 11)
Provider Vector
UPMC Children’s Hospital of Pittsburgh UPMC Magee-Womens Hospital UPMC Mercy UPMC Montefiore UPMC Passavant - McCandless UPMC Presbyterian UPMC Shadyside UPMC St. Margaret UPMC Western Psychiatric Hospital
Population 150,123.3 183,135.2 297,921.2 195,065.5 15,330.77 193,332.4 141,592.3 78,668.8 168,567

The product of that multiplication can then be modified to incorporate a value for supply. In this example, we’re just looking for accessibility to a hospital, so each hospital will receive a supply value of 1. It is the 1 hospital at that location. However, if the analysis was designed to use a different metric of supply, like number of buprenorphine prescriptions at that healthcare provider or the number of patients that provider can see in a given time period, then this is the step to add that information.

The result of this step is a 1 row by x columns vector with a value that represents numbers of hospitals per 100,000 residents for each of the nine hospitals in the data.

# Create supply for provider (use 1 for example)
supply <- 1

# Create provider to population ratio (per 100,000 residents)
for(k in 1:ncol(provider_vector)){
  if(provider_vector[1,k] != 0){
    provider_vector[1,k] <- ((100000 * supply) / provider_vector[1,k])
  }
}

5. Spatial Access Index

Calculate Spatial Access Index

The value that comes out of this analysis is a score called Spatial Access Index. This value is generated for each of the geographic units (e.g. each ZIP code in Allegheny county) and can be interpreted as the higher the score, the more access that geographic unit has to a healthcare provider.

SPAI is calculated by multiplying the 1 row by x columns provider vector from above with the x rows by y columns weight matrix.

# Spatial Access Index (SPAI)
SPAI_vector <- provider_vector %*% weight_matrix

Calculate Spatial Access Ratio

However, this value is difficult to interpret, and is somewhat meaningless without knowing the full distribution of values in the region. To resolve both of those issues, the next step of this type of analysis is to usually calculate a Spatial Access Ratio (SPAR), which normalizes SPAI values by dividing them by the average SPAI value for the entire region. A SPAR value of exactly 1 means that that geographic unit has exactly average access to a healthcare provider for that region. Geographic units with a SPAR above 1 have greater access than the regional average and geographic units with SPAR values below 1 have less access than the regional average. SPAR values also indicate rank, so if one geographic unit has a higher SPAR than another geographic unit, then the geographic unit with the higher SPAR value has greater access to a healthcare provider.

The step below calculates SPAR for each ZIP code.

# Format SPAI values and calculate Spatial Access Ratio (SPAR)
SPAI <- as.data.frame(t(SPAI_vector)) %>%
  rename(SPAI = Population) %>%
  rownames_to_column(var = "ZIP") %>%
  mutate(SPAR = SPAI / mean(SPAI)) %>%
  select(ZIP, SPAI, SPAR)

# Display first 10 values
SPAI[1:10,] %>%
  kbl(caption = "Check SPAI and SPAR Values for 10 ZIP Codes", 
      align = "r") %>%
  kable_minimal("hover")
Check SPAI and SPAR Values for 10 ZIP Codes
ZIP SPAI SPAR
15001 0.0000000 0.0000000
15003 0.0000000 0.0000000
15005 0.0000000 0.0000000
15007 0.2966425 0.3886845
15012 0.0000000 0.0000000
15014 0.1448114 0.1897434
15015 0.2966425 0.3886845
15017 0.0600278 0.0786532
15018 0.0000000 0.0000000
15024 0.1448114 0.1897434

Map SPAI and SPAR values

The easiest way to understand these values, and to check if they intuitively make sense, is to plot the values on a map and see if population centers and locations of service providers seem to impact spatial access values. The code below creates an interactive map that can display either SPAI or SPAR values for Allegheny County ZIP codes (using the Census Bureau’s ZCTA shapefiles). Although the values are different between the two variables, both retain the same rank of ZIP codes for spatial access.

# Create SPAI and SPAR values shapefile
SPAR_map <- left_join((populations %>%
                         mutate(ZIP = as.character(ZIP))),
                      SPAI,
                      by = "ZIP")


# Define color palettes for map
spar_palette <- colorNumeric(palette = "viridis",
                             domain = SPAR_map$SPAR)

spai_palette <- colorNumeric(palette = "magma",
                             domain = SPAR_map$SPAI)

# Create interactive SPAR and SPAI map
SPAR_map %>%
  leaflet() %>%
  addProviderTiles(provider = "CartoDB.Positron") %>%
  
  # Add SPAR values
  addPolygons(stroke = TRUE,
              weight = 1,
              color = "white",
              smoothFactor = 0,
              fillOpacity = 0.7,
              fillColor = ~spar_palette(SPAR),
              popup = paste0("<b>ZIP Code: ",
                             SPAR_map$ZIP,
                             "</b><br>SPAR :",
                             round(SPAR_map$SPAR,2),
                             "<br>Allegheny County Population: ",
                             prettyNum(SPAR_map$Allegheny_Population,
                                       big.mark = ",")),
              group = "SPAR") %>%
  
  # Add SPAI values
  addPolygons(stroke = TRUE,
              weight = 1,
              color = "white",
              smoothFactor = 0,
              fillOpacity = 0.7,
              fillColor = ~spai_palette(SPAI),
              popup = paste0("<b>ZIP Code: ",
                             SPAR_map$ZIP,
                             "</b><br>SPAI :",
                             round(SPAR_map$SPAI,2),
                             "<br>Allegheny County Population: ",
                             prettyNum(SPAR_map$Allegheny_Population,
                                       big.mark = ",")),
              group = "SPAI") %>%
  
  # Add Pittsburgh UPMC hospital locations
  addCircleMarkers(lng = upmc_hospitals$Longitude,
                   lat = upmc_hospitals$Latitude,
                   radius = 5,
                   stroke = TRUE,
                   color = "black",
                   fillOpacity = .75,
                   fillColor = "red",
                   popup = paste0("<b>",
                                  upmc_hospitals$Location,
                                  "</b><br>",
                                  upmc_hospitals$Address),
                   group = c("Hospitals")) %>%
  
  # Add legends
  addLegend("topright",
            pal = spar_palette,
            values = ~SPAR,
            opacity = 1,
            group = "SPAR") %>%
  addLegend("topright",
            pal = spai_palette,
            values = ~SPAI,
            opacity = 1,
            group = "SPAI") %>%
  
  # Add layers control (and default "SPAI" to off)
  addLayersControl(position = "topleft",
                   options = layersControlOptions(collapsed = FALSE,
                                                  autoZIndex = FALSE),
                   baseGroups = c("SPAR","SPAI"),
                   overlayGroups = c("Hospitals")) %>%
  hideGroup("SPAI")

Note: ZCTA shapefiles from 2010, which were used in place of ZIP codes and are only updated once every 10 years, can contain holes in physical space covered. This is why there are several small spots in the map above that appear black but are surrounded by SPAI or SPAR values. There is a comment about this under “Key Differences […]” on this page.

6. Details to consider

The enhanced two-step floating catchment area analysis is a useful way to quantify spatial access to healthcare providers. However, there are several details that can be modified during the calculations that can have a big impact on the final results.

Decisions on these choices should be made with an informed understanding of how the relationship between that specific type of healthcare provider and the specific population in the study interact with each other. This information should be useful for determinig how big of a catchment area is used, how many subzones are in the catchment area, what the time ranges of each subzone are, and how weights are determined for calculating SPAI. These choices can have a meaningful impact on the final results, so they should be based on an understanding of the scenario being measured.