Using GAM for hotspot mapping

Question

I’m going to preface this with, I’m no Generalized Additive Model or spatial statistics expert … I know enough to be “dangerous” … maybe that’s a bad thing 😉

Today I had a little bit of a thought experiment while thinking about Hotspot mapping. I’ve been wading pretty deeply into the spatio-temporal Generalized Additive Model (GAM) world recently and have been wondering if a GAM framework could be used to identify hotspots. Generally, in spatial statistics hotspots are defined as an area where values are greater than the mean. Typically, this is performed using the Getis-Ord Gi statistic which generally is a modified z-score to account for nearest neighbor distances in a spatial network. The Getis-Ord statistic can also be used to identify cold spots (values less than the mean).

If we have a generalized GAM model like

yi ~ s(yr) + s(x,y) + ti(x,y,yr)

where we have a time (yr), location (x,y), and an interaction term (x,y,yr) between location and time to evaluate changes in space and time. Given this model, the location effects plot shows spatially how spatially the data varies. Couldn’t this be used to identify locations that are (statistically) significantly greater or less than the mean (or some scaled value)? I know when attempting to detect significant changes in a time series from a GAM you need to use the derivative to determine if its significantly different from 0 (...i think) based on Dr Simpson's post. Would something similar have to be done with the s(x,y) term?

---

Here is some basic R-code with a reproducible example to get the ball rolling. In the code, there is a section that includes a Getis-Ord Gi statistic calculation. The GAM section of the code is relatively short as I'm not sure what is the next best step...as seen in the effects plots there are some negative and positive regions that correspond pretty closely to the Getis-Ord statistics.

## Spatial Data/Analysis
library(raster)
library(sf)
library(spdep)
library(spatstat)

## go-to GAM magic
library(mgcv)
library(gratia)

## Make some data 
## https://www.r-bloggers.com/2020/02/spatial-predictions-with-gams-and-rasters/

rbase <- raster(extent(c(0,1,0,1)), nrow = 50, ncol = 50)
rx1 <- rx2 <- rbase
rx1[] <- xFromCell(rbase, 1:ncell(rbase))
rx2[] <- yFromCell(rbase, 1:ncell(rbase))
rtrue <- 6*rx1  - 7*rx1^2- 4*rx2
par(mfrow = c(2,2))
plot(rx1, col = RColorBrewer::brewer.pal(11, "RdBu"), main = "x1")
plot(rx2, col = RColorBrewer::brewer.pal(11, "RdBu"), main = "x2")
plot(rtrue, col = RColorBrewer::brewer.pal(11, "RdBu"), main = "True values")

set.seed(42)
site_means <- data.frame(x = runif(50), y = runif(50)) |>
  mutate(site = 1:50,
         x1 = extract(rx1, cbind(x,y)),
         x2 = extract(rx2, cbind(x,y)),
         eta = extract(rtrue, cbind(x,y)),
         z = rnorm(50, sd = 1),
         yhat = eta + z)

dat <- merge(site_means,expand.grid(site = 1:50, yr = 1:5),"site")|>
  mutate(b = rnorm(250, mean = yhat, sd = 0.5))
dat$b<-dat$b+5 #(added 5 to have positive values ... I know I could just do a rlnorm)

## Getis Ord analysis
## a ref https://swampthingecology.org/blog/hot-spot-analysis-geospatial-data-analysis-in-rstats.-part-3/
# turn data into a spatial point data type
# the spstats package needs spatial data to do the nearest neighbor analysis

p.sf <- st_as_sf(dat, coords = c("y", "x"), crs = 4326) 

# Find distance range
ptdist=pointDistance(p.sf)

min.dist<-min(ptdist); # Minimum
mean.dist<-mean(ptdist); # Mean

nb<-dnearneigh(st_coordinates(p.sf),min.dist,mean.dist)

nb_lw<-nb2listw(nb,style="B")

## Global GI
globalG.test(p.sf$b,nb_lw,alternative="two.sided")

# local G
nb_lw<-nb2listw(nb)
local_g<-localG(p.sf$b,nb_lw)

# convert to matrix
local_g=as.matrix(local_g)

# column-bind the local_g data
p.sf<-cbind(p.sf,local_g)

# two side ρ-value
p.sf$pval<- 2*pnorm(-abs(p.sf$local_g))

# should identify several significant hotspots (points)

## simple GAM
m1=gam(b~
         s(yr,k=5,bs="tp")+
         s(x,y,bs="ds")+
         ti(x,y,yr,d=c(2,1),bs=c("ds","tp")),
       data=dat)
summary(m1)

draw(m1)

Answer 1

I think you can do this in principle, but it will depend on how you want to calculate the mean. For example, if you have data over space and time, you may want to compare against the mean across all timepoints. Or you may want to compare within timepoints, i.e. which spatial areas have the highest or lowest relative expectations within a given year? Either way, I believe you can use the {marginaleffects} package to help facilitate these calculations. This package makes it very easy to set up prediction grids and run hypothesis tests using model predictions and associated uncertainties.

Below I give an example of how this could be done by simply creating a 25 x 25 grid of points to make predictions over for each year in the simulated data. These are fed to the predictions() function to calculate a mean expectation over all timepoints, and then the hypotheses() function is used to identify locations in space and time that differ significantly (based on a p-value of 0.05) from the mean after adjusting for multiple comparisons. Doing it this way is important because, as your model illustrates, you can easily have spatial and/or temporal terms that are spread across multiple smooth functions in a GAM in {mgcv}. So you cannot get a full sense of what spatial or temporal patterns the model expects to see by simply looking at the smooths in isolation. Predictions are key.

Obviously a dedicated user would need to think carefully about what spatial and temporal grids to predict over, what kind of hypothesis to compute, and how to adjust the p-values appropriately, but I think this is easily doable in {marginaleffects}.

# Libraries
library(raster)
library(sf)
library(spdep)
library(spatstat)
library(dplyr)
library(marginaleffects)
library(ggplot2); theme_set(theme_bw())
library(viridis)
library(mgcv)

# Data 
rbase <- raster(extent(c(0,1,0,1)), nrow = 50, ncol = 50)
rx1 <- rx2 <- rbase
rx1[] <- xFromCell(rbase, 1:ncell(rbase))
rx2[] <- yFromCell(rbase, 1:ncell(rbase))
rtrue <- 6*rx1  - 7*rx1^2- 4*rx2

set.seed(42)
site_means <- data.frame(x = runif(50), y = runif(50)) |>
  mutate(site = 1:50,
         x1 = extract(rx1, cbind(x,y)),
         x2 = extract(rx2, cbind(x,y)),
         eta = extract(rtrue, cbind(x,y)),
         z = rnorm(50, sd = 1),
         yhat = eta + z)

dat <- merge(site_means,expand.grid(site = 1:50, yr = 1:5),"site") |>
  mutate(b = rnorm(250, mean = yhat, sd = 0.5))
dat$b<-dat$b+5

# Ord analysis
p.sf <- st_as_sf(dat, coords = c("y", "x"), crs = 4326) 

# Find distance range
ptdist=pointDistance(p.sf)
min.dist<-min(ptdist); # Minimum
mean.dist<-mean(ptdist); # Mean
nb<-dnearneigh(st_coordinates(p.sf),min.dist,mean.dist)
nb_lw<-nb2listw(nb,style="B")

# Global GI
globalG.test(p.sf$b,nb_lw,alternative="two.sided")
#> 
#>  Getis-Ord global G statistic
#> 
#> data:  p.sf$b 
#> weights: nb_lw 
#> 
#> standard deviate = 2.2934, p-value = 0.02183
#> alternative hypothesis: two.sided
#> sample estimates:
#> Global G statistic        Expectation           Variance 
#>       5.235505e-01       5.044177e-01       6.960063e-05

# Local G
nb_lw<-nb2listw(nb)
local_g<-localG(p.sf$b,nb_lw)

# Convert to matrix
local_g=as.matrix(local_g)

# Column-bind the local_g data
p.sf<-cbind(p.sf,local_g)

# Two side ρ-value
p.sf$pval<- 2*pnorm(-abs(p.sf$local_g))

# Plot the observed data over space and time
ggplot(dat, aes(x, y, color = b)) +
  scale_colour_viridis() +
  geom_point() +
  facet_wrap(~yr)

# GAM
m1=gam(b~
         s(yr,k=5,bs="tp")+
         s(x,y,bs="ds")+
         ti(x,y,yr,d=c(2,1),bs=c("ds","tp")),
       data=dat)
summary(m1)
#> 
#> Family: gaussian 
#> Link function: identity 
#> 
#> Formula:
#> b ~ s(yr, k = 5, bs = "tp") + s(x, y, bs = "ds") + ti(x, y, yr, 
#>     d = c(2, 1), bs = c("ds", "tp"))
#> 
#> Parametric coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)  3.73047    0.04852   76.88   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Approximate significance of smooth terms:
#>              edf Ref.df      F p-value    
#> s(yr)       1.00   1.00  0.127   0.722    
#> s(x,y)     28.89  31.42 24.960  <2e-16 ***
#> ti(x,y,yr)  2.00   2.00  0.287   0.751    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> R-sq.(adj) =  0.756   Deviance explained = 78.7%
#> GCV = 0.67781  Scale est. = 0.58865   n = 250

# Get expected values across all times and a grid of spatial locations
# using marginaleffects
even_seq = function(x, length.out = 25){
  minx <- min(x, na.rm = TRUE)
  maxx <- max(x, na.rm = TRUE)
  seq(minx, maxx, length.out = length.out)
}

newdata <- datagrid(model = m1,
                    x = even_seq,
                    y = even_seq,
                    yr = even_seq(dat$yr, length.out = 5))

preds <- predictions(m1, newdata = newdata, type = 'response')

# Calculate mean across all predictions
mean_pred <- mean(preds$estimate)

# Use the delta method to calculate approximate standard errors and then
# conduct an F test against the null hypothesis that all expected values are equal
# to the mean
hyps <- hypotheses(preds,
                   hypothesis = mean_pred) |>
  
  # Run Bonferoni adjustment for multiple comparisons
  mutate(p.value = p.adjust(p.value)) |>
  
  # Calculate significance and classify into hot and cold spots
  mutate(significant = p.value <= 0.05,
         classification = case_when(
           estimate < mean_pred & significant ~ 'coldspot',
           estimate > mean_pred & significant ~ 'hotspot',
           TRUE ~ '-'))

# Plot the classifications
cols <- c('darkred', 'darkblue', 'grey')
names(cols) <- c('hotspot', 'coldspot', '-')
ggplot(hyps, aes(x, y, color = classification)) +
  geom_point() +
  scale_color_manual(values = cols) +
  facet_wrap(~yr)

^{Created on 2024-05-22 with reprex v2.0.2}