Section 0: Setup & Install

Run the following cell to install and load the necessary spatial packages for this workshop. Note: Google Colab does not pre-install these R spatial packages, so this might take a few minutes.

install.packages(c('sf', 'spdep', 'tmap', 'MASS', 'spmoran', 'ggplot2', 'dplyr', 'RColorBrewer'))
library(sf); library(spdep); library(tmap); library(MASS); library(ggplot2); library(dplyr)
tmap_mode('plot')

Section 1: Load & Explore Data

Download the Miami-Dade COVID-19 dataset from OSF and read it into an sf (Simple Features) object. We use GEOM_POSSIBLE_NAMES=geom so sf automatically identifies the WKT geometry column.

download.file('https://files.osf.io/v1/resources/tjcny/providers/osfstorage/685441a5f40bd63220a1c9a3?view_only=e598d81c1c5447c2a4df998a25f58825', 'MiamiDade_UDB_CBG_covid.csv')

# Read data and filter out block groups with 0 population to avoid calculation errors
data_sf <- st_read('MiamiDade_UDB_CBG_covid.csv', options = 'GEOM_POSSIBLE_NAMES=geom')
data_sf <- data_sf %>% filter(total_pop > 0)

Let’s look at the structure and dimensions of our spatial dataframe:

dim(data_sf)
names(data_sf)
summary(data_sf[, c('total_pop', 'median_household_income', 'pct_below_poverty', 'pct_hispanic_or_latino')])

Basic Exploratory Data Analysis

Let’s visualize the statistical distribution of cases before we map them.

ggplot(data_sf, aes(x = case_period3)) +
  geom_histogram(fill = 'steelblue', color = 'black', bins=30) +
  theme_minimal() +
  labs(title = 'Histogram of COVID-19 Cases (Period 3: Sep-Dec 2020)')

# Boxplot of rates across a few periods
boxplot(data_sf$rate_period1, data_sf$rate_period3, data_sf$rate_period5, names=c('Period 1', 'Period 3', 'Period 5'), main='Infection Rates by Pandemic Phase')

Section 2: Thematic Mapping

Let’s map the COVID-19 smoothed rates for Period 3 (Fall 2020) along with some key predictors.

Discussion: Do you visually notice spatial correlations? Are the highest COVID-19 rates in the same neighborhoods as the lowest incomes?

tm_shape(data_sf) +
  tm_polygons('smoothed_rate_period3', style='quantile', palette='Reds', title='Smoothed Rate P3') +
  tm_layout(legend.outside = TRUE, title = 'COVID-19 Rates - Fall 2020')

# YOUR TURN: Try changing variables to other demographics like 'pct_black_or_african_american' or 'pct_below_poverty'
tm_shape(data_sf) +
  tm_polygons(c('pct_hispanic_or_latino', 'median_household_income'), style='quantile', title=c('% Hispanic/Latino', 'Median HH Income')) +
  tm_layout(legend.outside = TRUE)

Section 3: Spatial Weights Matrix

To run spatial statistics, we need to define ‘neighbors’. We’ll use Queen contiguity (sharing borders or corners).

nb <- poly2nb(data_sf, queen=TRUE)
summary(nb)

Let’s visualize the connectivity graph overlaid on the census block groups.

coords <- st_coordinates(st_centroid(st_geometry(data_sf)))
plot(st_geometry(data_sf), border='lightgrey')
plot(nb, coords, add=TRUE, col='red', lwd=0.5)
title('Queen Contiguity Neighborhood Graph')

Row-standardize the weights. We use zero.policy=TRUE to gracefully handle ‘islands’ that have no neighbors.

listw <- nb2listw(nb, style='W', zero.policy=TRUE)
summary(listw, zero.policy=TRUE)

Section 4: Global Spatial Autocorrelation

Calculate Global Moran’s I to test for spatial clustering of COVID-19 rates across the entire county.

moran_test <- moran.test(data_sf$smoothed_rate_period3, listw, zero.policy=TRUE)
print(moran_test)

The Moran Scatterplot plots the standardized variable against its spatial lag (the average of its neighbors). The slope of the regression line is Moran’s I.

moran.plot(data_sf$smoothed_rate_period3, listw, zero.policy=TRUE, pch=20, col='darkblue')

# General G statistic to test clustering of High vs Low values
globalG.test(data_sf$smoothed_rate_period3, listw, zero.policy=TRUE)

Section 5: Local Spatial Autocorrelation — LISA

While Global Moran’s I tells us if there is clustering, Local Indicators of Spatial Association (LISA) tell us where the clusters are.

High-High: High rate surrounded by high rates (Hotspots)
Low-Low: Low rate surrounded by low rates (Coldspots)
High-Low / Low-High: Spatial outliers

lisa <- localmoran(data_sf$smoothed_rate_period3, listw, zero.policy=TRUE)
signif <- lisa[, 5] < 0.05
z_val <- scale(data_sf$smoothed_rate_period3)
lag_z_val <- lag.listw(listw, z_val, zero.policy=TRUE)

quad <- rep('Not Significant', length(z_val))
quad[z_val > 0 & lag_z_val > 0 & signif] <- 'High-High'
quad[z_val < 0 & lag_z_val < 0 & signif] <- 'Low-Low'
quad[z_val > 0 & lag_z_val < 0 & signif] <- 'High-Low'
quad[z_val < 0 & lag_z_val > 0 & signif] <- 'Low-High'

data_sf$lisa_cluster <- factor(quad, levels=c('High-High', 'Low-Low', 'High-Low', 'Low-High', 'Not Significant'))

Map the LISA clusters to locate precise neighborhood dynamics.

tm_shape(data_sf) +
  tm_polygons('lisa_cluster', palette=c('red', 'blue', 'pink', 'lightblue', 'grey90')) +
  tm_layout(legend.outside = TRUE, title='LISA Clusters Period 3')

Section 6: Hot Spot Analysis — Getis-Ord Gi*

Getis-Ord Gi* specifically identifies statistically significant Hot and Cold Spots based on z-scores.

# Gi* requires weights that include the self-neighbor
nb_self <- include.self(nb)
listw_self <- nb2listw(nb_self, style='W', zero.policy=TRUE)

gi_star <- localG(data_sf$smoothed_rate_period3, listw_self, zero.policy=TRUE)
data_sf$gi_star_z <- as.numeric(gi_star)

Map the calculated Gi* z-scores (Values > 1.96 or < -1.96 are significant at p < 0.05).

tm_shape(data_sf) +
  tm_polygons('gi_star_z', style='cont', midpoint=0, palette='-RdBu', title='Gi* Z-Score') +
  tm_layout(legend.outside=TRUE, title='Getis-Ord Gi* Z-Scores')

Compare our computed result with the pre-computed Gi* clusters.

Observation: Notice the strong cold spots running through the urban core along the I-95 corridor.

tm_shape(data_sf) +
  tm_polygons('gi_star_cluster_type_period3', palette=c('High'='red', 'Low'='blue', 'Not Significant'='grey90')) +
  tm_layout(legend.outside=TRUE, title='Pre-computed Gi* Clusters')

Section 7: Negative Binomial Regression

We model the case counts adjusting for socioeconomic factors. Since case data is count-based with overdispersion, Negative Binomial is preferred over Poisson. We use offset(log(total_pop)) to properly model rates instead of raw counts.

nbm_model <- glm.nb(
  case_period3 ~ offset(log(total_pop)) + 
    pct_hispanic_or_latino + 
    pct_black_or_african_american + 
    pct_below_poverty + 
    pct_no_health_insurance + 
    median_household_income + 
    pct_no_high_school + 
    pct_2_01_or_more_occupants_per_room + 
    pop_dens, 
  data = data_sf
)
summary(nbm_model)

Exponentiate the coefficients to calculate Incidence Rate Ratios (IRRs) for easier interpretation.

irr_table <- data.frame(
  IRR = exp(coef(nbm_model)),
  Lower_CI = suppressWarnings(exp(confint(nbm_model)[,1])),
  Upper_CI = suppressWarnings(exp(confint(nbm_model)[,2]))
)
print(round(irr_table, 3))

Are the residuals spatially autocorrelated? If so, the standard model violates the assumption of independent errors, leading to biased significance tests.

data_sf$nbm_resid <- residuals(nbm_model, type='deviance')

# Test residuals with Moran's I
moran_test_resid <- moran.test(data_sf$nbm_resid, listw, zero.policy=TRUE)
print(moran_test_resid)

# Map residuals to visually check for spatial patterns
tm_shape(data_sf) +
  tm_polygons('nbm_resid', style='quantile', palette='-RdBu', midpoint=0, title='Deviance Residuals') +
  tm_layout(legend.outside=TRUE, title='NBM Residual Pattern')

Section 8: Spatial Eigenvector Filtering (SEF)

Our residuals showed significant spatial clustering, which means the model violates the assumption of independent errors. This can lead to biased significance tests and unreliable coefficient estimates.

Spatial Eigenvector Filtering addresses this by: 1. Decomposing the spatial weights matrix into eigenvectors that represent distinct spatial patterns 2. Selecting a parsimonious subset of eigenvectors via stepwise AIC (to avoid overfitting) 3. Including them as additional predictors to absorb the spatial structure in the residuals

This is conceptually similar to adding ‘spatial control variables’ — the eigenvectors capture latent spatial processes that our substantive predictors missed.

Step 1: Compute eigenvectors from the centered spatial weights matrix.

# Build binary weights matrix and center it
W <- nb2mat(nb, style='B', zero.policy=TRUE)
n <- nrow(W)
MCM <- diag(n) - matrix(1/n, n, n)  # Centering matrix
MWM <- MCM %*% W %*% MCM

# Eigendecomposition — this may take a moment for large datasets
cat('Computing eigendecomposition of', n, 'x', n, 'matrix...\n')
eig <- eigen(MWM, symmetric=TRUE)

# Keep candidate eigenvectors with positive eigenvalues
# (these represent positive spatial autocorrelation patterns)
pos_idx <- which(eig$values > 0.001)
n_candidates <- min(50, length(pos_idx))
EV <- eig$vectors[, pos_idx[1:n_candidates]]
colnames(EV) <- paste0('EV', 1:n_candidates)

cat(sprintf('Candidate pool: %d eigenvectors with positive spatial autocorrelation\n', n_candidates))

# Add candidates to the dataset
data_sf <- cbind(data_sf, EV)

Step 2: Use stepwise AIC selection to find the most parsimonious set of eigenvectors.

Rather than blindly including all eigenvectors (which would overfit), we let AIC decide which ones meaningfully improve the model. AIC penalizes each additional parameter, so only eigenvectors that capture substantial spatial structure will be retained.

# Start from the full model with all 50 candidates and let stepwise AIC prune
base_form <- 'case_period3 ~ offset(log(total_pop)) + pct_hispanic_or_latino + pct_black_or_african_american + pct_below_poverty + pct_no_health_insurance + median_household_income + pct_no_high_school + pct_2_01_or_more_occupants_per_room + pop_dens'
ev_terms <- paste0('EV', 1:n_candidates, collapse=' + ')
full_form <- as.formula(paste(base_form, '+', ev_terms))

cat('Fitting full model with all candidate eigenvectors...\n')
full_model <- glm.nb(full_form, data=data_sf, control=glm.control(maxit=100))

cat('Running stepwise AIC selection (this may take a minute)...\n')
sef_model <- step(full_model, direction='both', trace=0)

# How many eigenvectors were retained?
retained_evs <- grep('^EV', names(coef(sef_model)), value=TRUE)
cat(sprintf('\nStepwise retained %d of %d candidate eigenvectors: %s\n', 
    length(retained_evs), n_candidates, paste(retained_evs, collapse=', ')))

cat(sprintf('\nBaseline NBM AIC: %.1f\n', AIC(nbm_model)))
cat(sprintf('SEF Model AIC:    %.1f  (improvement: %.1f)\n', AIC(sef_model), AIC(nbm_model) - AIC(sef_model)))

Step 3: Check if SEF reduced the spatial autocorrelation in the residuals.

data_sf$sef_resid <- residuals(sef_model, type='deviance')

# Compare Moran's I before and after SEF
moran_after <- moran.test(data_sf$sef_resid, listw, zero.policy=TRUE)

cat('=== Before SEF ===\n')
cat(sprintf('Moran I = %.4f, p = %s\n\n', moran_test_resid$estimate[1], format.pval(moran_test_resid$p.value, digits=3)))
cat('=== After SEF ===\n')
cat(sprintf('Moran I = %.4f, p = %s\n\n', moran_after$estimate[1], format.pval(moran_after$p.value, digits=3)))
cat(sprintf('Reduction in Moran I: %.1f%%\n', 
    (1 - moran_after$estimate[1]/moran_test_resid$estimate[1]) * 100))

Step 4: Map the residuals side by side.

m1 <- tm_shape(data_sf) +
  tm_polygons('nbm_resid', style='quantile', palette='-RdBu', midpoint=0, title='Residuals') +
  tm_layout(legend.outside=TRUE, title='Before SEF')

m2 <- tm_shape(data_sf) +
  tm_polygons('sef_resid', style='quantile', palette='-RdBu', midpoint=0, title='Residuals') +
  tm_layout(legend.outside=TRUE, title='After SEF')

tmap_arrange(m1, m2, ncol=2)

A note on statistical vs. practical significance

You will likely notice that the Moran’s I p-value is still statistically significant after SEF. Is that a problem?

Not necessarily. With ~1,600 spatial units, even very small Moran’s I values (e.g., 0.03) will be statistically significant due to high statistical power. What matters more is:

The magnitude of reduction — a 50%+ drop in Moran’s I is substantial
Model parsimony — we achieved this with only ~8-15 eigenvectors, not hundreds
AIC improvement — lower AIC confirms the spatial terms genuinely improve fit
Visual comparison — the residual maps should show noticeably less spatial clustering

Key takeaway: With large spatial datasets, focus on the magnitude of spatial autocorrelation (Moran’s I value), not just the p-value. This distinction between statistical and practical significance becomes increasingly important as sample sizes grow.

Section 9: Summary & Discussion

Key Findings

Spatial autocorrelation in residuals confirms that standard regression models miss important spatial processes in disease data.
SEF with stepwise AIC selection provides a principled, parsimonious way to account for spatial dependence without overfitting.
Interpreting spatial patterns in residuals can reveal structural inequities — areas where the model consistently over- or under-predicts may reflect unmeasured access barriers.

The Bigger Picture

Spatial analysis is not just a statistical correction — it is a lens for understanding where and why health outcomes vary geographically. Methods like hot spot analysis, LISA, and spatial regression help us move from ‘what happened’ to ‘where should we intervene’.

Thank you for completing the workshop! Check out more interactive tools at GatorLab-Geo.