Getis-Ord Gi* Hot Spot Analysis

Spatiotemporal Trends in COVID-19 Cases

This notebook is designed for the RCMI 2026 workshop on spatial epidemiology. It focuses on identifying and analyzing hot and cold spots using the Getis-Ord Gi* statistic across multiple time periods.

Section 0: Setup & Data Loading

First, we will install and load the necessary packages.

# Install packages (Colab doesn't have these pre-installed for R)
install.packages(c("sf", "spdep", "tmap", "dplyr", "ggplot2", "RColorBrewer", "tidyr"), repos="http://cran.us.r-project.org")

# Load libraries
library(sf)
library(spdep)
library(tmap)
library(dplyr)
library(ggplot2)
library(RColorBrewer)
library(tidyr)

# Download dataset
if(!file.exists('MiamiDade_UDB_CBG_covid.csv')) {
  download.file('https://files.osf.io/v1/resources/tjcny/providers/osfstorage/685441a5f40bd63220a1c9a3?view_only=e598d81c1c5447c2a4df998a25f58825', 'MiamiDade_UDB_CBG_covid.csv')
}

# Read CSV and convert to spatial (sf) object
df <- read.csv('MiamiDade_UDB_CBG_covid.csv')
cbg_sf <- st_as_sf(df, wkt = "geom", crs = 4326)

# Drop empty geometries or 0 population
cbg_sf <- cbg_sf %>% filter(!st_is_empty(geom), total_pop > 0)

print(paste("Number of valid CBGs:", nrow(cbg_sf)))

Section 1: What is Gi*?

The Getis-Ord Gi* statistic identifies statistically significant spatial clusters of high values (hot spots) and low values (cold spots). It looks at each feature within the context of neighboring features.

Formula Overview: The local sum for a feature and its neighbors is compared proportionally to the sum of all features; when the local sum is very different from the expected local sum, and that difference is too large to be the result of random chance, a statistically significant z-score results.

• Positive Gi* z-score → hot spot (high values surrounded by high values)
• Negative Gi* z-score → cold spot (low values surrounded by low values)

Significance is typically assessed at 90%, 95%, or 99% confidence levels.

Section 2: Building Spatial Weights

To compute Gi*, we first need to define “neighbors.” We’ll use Queen contiguity (sharing boundaries or vertices).

**Crucial Note for Gi*:** Unlike the \(G_i\) statistic, the \(G_i^*\) statistic includes the feature itself in its own neighbor list. This means we must add self-loops to our neighbor graph.

# 1. Queen contiguity neighbors
nb <- poly2nb(cbg_sf, queen = TRUE)

# 2. Include self (crucial for Gi*)
nb_self <- include.self(nb)

# 3. Create spatial weights
lw_self <- nb2listw(nb_self, style = "W", zero.policy = TRUE)

# Let's verify that self is included
print(nb[[1]])
print(nb_self[[1]])

Section 3: Computing Gi* for Period 3

Let’s compute the Gi* statistic for the raw COVID-19 rates in Period 3 (Sep-Dec 2020, fall wave).

# Compute Gi* (using localG from spdep)
# localG returns z-scores
cbg_sf$gi_star_z_p3 <- localG(cbg_sf$rate_period3, lw_self)
cbg_sf$gi_star_z_p3 <- as.numeric(cbg_sf$gi_star_z_p3) # convert to standard numeric

# Classify based on z-score
# Z-score thresholds for significance (approx):
# 99%: > 2.58 or < -2.58
# 95%: > 1.96 or < -1.96
# 90%: > 1.645 or < -1.645

cbg_sf <- cbg_sf %>%
  mutate(
    gi_class_p3 = case_when(
      gi_star_z_p3 > 2.58 ~ "Hot Spot (99%)",
      gi_star_z_p3 > 1.96 ~ "Hot Spot (95%)",
      gi_star_z_p3 > 1.645 ~ "Hot Spot (90%)",
      gi_star_z_p3 < -2.58 ~ "Cold Spot (99%)",
      gi_star_z_p3 < -1.96 ~ "Cold Spot (95%)",
      gi_star_z_p3 < -1.645 ~ "Cold Spot (90%)",
      TRUE ~ "Not Significant"
    )
  )

# Set factor levels for plotting
gi_levels <- c("Cold Spot (99%)", "Cold Spot (95%)", "Cold Spot (90%)", 
               "Not Significant", 
               "Hot Spot (90%)", "Hot Spot (95%)", "Hot Spot (99%)")
cbg_sf$gi_class_p3 <- factor(cbg_sf$gi_class_p3, levels = gi_levels)

# Plot Z-scores
tm_shape(cbg_sf) +
  tm_polygons("gi_star_z_p3", 
              style = "cont", 
              midpoint = 0, 
              palette = "-RdBu", 
              title = "Gi* Z-Score (Period 3)", border.col = "transparent")

# Map the classified spots
gi_colors <- c("#4575b4", "#91bfdb", "#e0f3f8", "#ffffbf", "#fee090", "#fc8d59", "#d73027")

tm_shape(cbg_sf) +
  tm_polygons("gi_class_p3", 
              palette = gi_colors, 
              title = "Gi* Classification",
              border.col = "transparent")

Section 4: Comparing with Pre-computed Gi*

The dataset already includes a column gi_star_cluster_type_period3. Let’s compare our manual classification with the pre-computed one.

table(Our_Classification = cbg_sf$gi_class_p3, Pre_computed = cbg_sf$gi_star_cluster_type_period3)

# Notice any differences? They might stem from different spatial weights matrices (e.g., k-nearest neighbors vs queen contiguity) or handling of edge cases (e.g., zero policy).

Section 5: Multi-Period Gi* Analysis

How did hot spots shift over time? Let’s compute Gi* for all 5 periods. - Period 1: Mar-May 2020 - Period 2: Jun-Aug 2020 - Period 3: Sep-Dec 2020 - Period 4: Jan-May 2021 - Period 5: Jun-Sep 2021

# Loop over periods to compute Gi*
for(i in 1:5) {
  rate_col <- paste0("rate_period", i)
  z_col <- paste0("gi_z_p", i)
  class_col <- paste0("gi_class_p", i)
  
  # Compute Z-score
  cbg_sf[[z_col]] <- as.numeric(localG(cbg_sf[[rate_col]], lw_self))
  
  # Classify
  cbg_sf[[class_col]] <- case_when(
    cbg_sf[[z_col]] > 2.58 ~ "Hot Spot (99%)",
    cbg_sf[[z_col]] > 1.96 ~ "Hot Spot (95%)",
    cbg_sf[[z_col]] > 1.645 ~ "Hot Spot (90%)",
    cbg_sf[[z_col]] < -2.58 ~ "Cold Spot (99%)",
    cbg_sf[[z_col]] < -1.96 ~ "Cold Spot (95%)",
    cbg_sf[[z_col]] < -1.645 ~ "Cold Spot (90%)",
    TRUE ~ "Not Significant"
  )
  cbg_sf[[class_col]] <- factor(cbg_sf[[class_col]], levels = gi_levels)
}

# Create faceted maps
maps <- list()
for(i in 1:5) {
  class_col <- paste0("gi_class_p", i)
  maps[[i]] <- tm_shape(cbg_sf) +
    tm_polygons(class_col, palette = gi_colors, border.col = "transparent", legend.show = (i == 5)) +
    tm_layout(title = paste("Period", i), title.position = c("left", "top"))
}

tmap_arrange(maps[[1]], maps[[2]], maps[[3]], maps[[4]], maps[[5]], ncol = 3)

Interpretation

• Early periods (1-2): Hot spots often clustered in specific dense, lower-income areas.
• Later periods (4-5): Patterns shift. As vaccines rolled out (Period 4) and Delta emerged (Period 5), hot spots may shift towards areas with lower vaccination rates or different demographics.
• Persistent cold spots: You may notice persistent cold spots in areas like Liberty City, Little Haiti, and Opa-Locka. However, due to structural barriers, this could represent a “detection paradox” where low testing artificially depresses case counts.

Section 6: Characterizing Hot Spots vs Cold Spots

Let’s look at the demographics of these clusters in Period 3.

# Summarize demographics by cluster type
summary_p3 <- cbg_sf %>%
  st_drop_geometry() %>%
  group_by(gi_class_p3) %>%
  summarize(
    count = n(),
    med_income = mean(median_household_income, na.rm = TRUE),
    pct_black = mean(pct_black_or_african_american, na.rm = TRUE),
    pct_hispanic = mean(pct_hispanic_or_latino, na.rm = TRUE),
    pct_uninsured = mean(pct_no_health_insurance, na.rm = TRUE)
  )

print(summary_p3)

# Boxplot comparing % Black by cluster type
ggplot(cbg_sf, aes(x = gi_class_p3, y = pct_black_or_african_american, fill = gi_class_p3)) +
  geom_boxplot() +
  scale_fill_manual(values = gi_colors) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "Percentage of Black/African American Population by Cluster Type (Period 3)",
       x = "Cluster Type", y = "% Black / African American")

Section 7: Emerging Hot Spot Analysis

We can classify areas based on their temporal behavior across the 5 periods. - Persistent Hot Spot: Hot in 4 or 5 periods - Persistent Cold Spot: Cold in 4 or 5 periods

# Count how many times an area was hot or cold
cbg_sf <- cbg_sf %>%
  mutate(
    hot_count = (gi_class_p1 %in% c("Hot Spot (90%)", "Hot Spot (95%)", "Hot Spot (99%)")) +
                (gi_class_p2 %in% c("Hot Spot (90%)", "Hot Spot (95%)", "Hot Spot (99%)")) +
                (gi_class_p3 %in% c("Hot Spot (90%)", "Hot Spot (95%)", "Hot Spot (99%)")) +
                (gi_class_p4 %in% c("Hot Spot (90%)", "Hot Spot (95%)", "Hot Spot (99%)")) +
                (gi_class_p5 %in% c("Hot Spot (90%)", "Hot Spot (95%)", "Hot Spot (99%)")),
                
    cold_count = (gi_class_p1 %in% c("Cold Spot (90%)", "Cold Spot (95%)", "Cold Spot (99%)")) +
                 (gi_class_p2 %in% c("Cold Spot (90%)", "Cold Spot (95%)", "Cold Spot (99%)")) +
                 (gi_class_p3 %in% c("Cold Spot (90%)", "Cold Spot (95%)", "Cold Spot (99%)")) +
                 (gi_class_p4 %in% c("Cold Spot (90%)", "Cold Spot (95%)", "Cold Spot (99%)")) +
                 (gi_class_p5 %in% c("Cold Spot (90%)", "Cold Spot (95%)", "Cold Spot (99%)"))
  ) %>%
  mutate(
    temporal_class = case_when(
      hot_count >= 4 ~ "Persistent Hot Spot",
      cold_count >= 4 ~ "Persistent Cold Spot",
      hot_count > 0 & cold_count > 0 ~ "Oscillating",
      TRUE ~ "Other / Transient"
    )
  )

table(cbg_sf$temporal_class)

tm_shape(cbg_sf) +
  tm_polygons("temporal_class", 
              palette = c("Oscillating" = "purple", "Other / Transient" = "gray90", 
                          "Persistent Cold Spot" = "blue", "Persistent Hot Spot" = "red"),
              border.col = "transparent", title = "Temporal Pattern")

Section 8: Discussion

• The Detection Paradox: Persistent cold spots in predominantly Black/African American neighborhoods (like Liberty City) often reflect severe barriers to testing access, rather than true low transmission. This shows why spatial statistics must be interpreted alongside socioeconomic context.
• Dynamic Shifts: The Oscillating category highlights how the pandemic’s geography was not static. Surges hit different communities at different times.
• Interactive Tools: To explore these metrics further, refer to the workshop’s interactive Shiny applications!