Coarse-to-fine spatial GLMM: Application examples

In spCF: Coarse-to-Fine Spatial and Spatio-Temporal Modeling

```{=html}

```r
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

This vignette demonstrates how to estimate spatial generalized linear mixed models (GLMMs) using the spCF package for smoothing/denoising, prediction, and multiscale analysis. The spatial process in the GLMMs is trained via a coarse-to-fine approach that sequentially learns spatial patterns from coarser to finer scales. This spatial GLMM framework is particularly suitable for moderate-to-large samples. See @murakami2026cfsm-glm for further details.

Example 1: Count data modeling for disease mapping

Data and setup

Let us load the required packages

library(spCF)
library(sf)
library(CARBayesdata)

The pollutionhealthdata dataset included in the CARBayesdata package is used here. This dataset comprises panel data on respiratory hospitalisations across 271 zones in Greater Glasgow. For simplicity, the center coordinates of each zone are used for spatial modeling.

In this example, the observed number of hospitalisations in 2011 (observed) is modeled for smoothing and predicting the latent risk of respiratory disease (i.e., disease mapping).

data(pollutionhealthdata)
dat      <- pollutionhealthdata[pollutionhealthdata$year==2011,]

The response variable (observed), covariates (pm10, jsa, price), and an offset variable (expected) are specified as follows:

y        <- dat[,"observed"]             # count data
x        <- dat[,c("pm10","jsa","price")]# covariates
offset   <- log(dat[,"expected"])        # offset variable

The center coordinates are extracted as:

data(GGHB.IZ)                            # polygons of the 271 zones
coords   <- st_coordinates(st_centroid(GGHB.IZ))# coordinates

observed is spatially plotted as follows:

GGHB.IZ$y <- y
plot(GGHB.IZ[,"y"],lwd=0.01,axes=TRUE, key.pos=4,nbreaks=50)

Coarse-to-fine spatial GLMM (CF-GLM)

In CF-GLM, the spatial process is defined as a sum of scalewise processes, where the number of spatial scales, R, is optimized via holdout validation. A smaller R, corresponding to early stopping, allows the spatial process to capture only coarse-scale patterns, whereas a larger R enables the spatial process to represent finer-scale patterns. The cf_glm_hv function performs holdout validation as follows:

mod_hv   <- cf_glm_hv(y = y, x = x, offset=offset, 
                      coords = coords, family=poisson())

As shown in the output, the deviance loss for the validation samples gradually decreases as learning proceeds from the coarsest scale (Scale 1) to finer scales.

After holdout validation, the full model is trained using the cf_glm function:

mod      <- cf_glm(y = y, x = x, offset=offset, 
                   coords = coords, mod_hv = mod_hv)

The estimated regression coefficients, standard deviations of the scalewise components, and error statistics are displayed as follows:

mod

Disease mapping

Predictive values

The predictive means, representing the denoised risk of respiratory disease, and their standard deviations (SDs) are mapped as follows:

# Predictive mean
GGHB.IZ$pred   <- mod$pred$pred
plot(GGHB.IZ[,c("pred")],lwd=0.01,axes=TRUE, key.pos=4,nbreaks=50)

# Predictive SD
GGHB.IZ$pred_sd<- mod$pred$pred_sd
plot(GGHB.IZ[,c("pred_sd")],pal = function(n) hcl.colors(n, "Viridis"), lwd=0.01, 
     axes=TRUE, key.pos=4)

The result suggests that both disease risk and their uncertainty increase in the central area.

Multiscale spatial pattern/feature extraction

The sp_scalewise function extracts scalewise spatial processes with bandwidth values falling within a pre-specified range. For example, the following commands extract the large- and small-scale processes, corresponding to bandwidth ranges of 4000+ and 0–4000, respectively.

mod_s1      <- sp_scalewise(mod,bw_range=c(4000,Inf)) # Large scale (4000 <= bandwidth)
mod_s2      <- sp_scalewise(mod,bw_range=c(0,4000))   # Small scale (bandwidth <= 4000)

The extracted scalewise processes are mapped as follows:

GGHB.IZ$z1  <- mod_s1$pred$pred
GGHB.IZ$z2  <- mod_s2$pred$pred
plot(GGHB.IZ[,c("z1","z2")],lwd=0.01,axes=TRUE,key.pos=4, nbreaks=50)

The sp_scalewise function is useful for multiscale spatial pattern analysis (or feature extraction), which is commonly conducted in ecological, epidemiological, and environmental studies.

Example 2: Binary data modeling and prediction

Data and setup

Let us load the required packages

library(spCF)
library(sp)
library(sf)

The meuse dataset from the sp package consists of observations at 155 sample sites in a floodplain along the River Meuse. In this example, we consider a binary response variable (flood), which takes the value 1 if the flooding frequency class is once every two years (i.e., flood-prone area) and 0 otherwise. The binary response is predicted over 3,103 regularly spaced grid cells (meuse.grid). The covariate considered is the distance to the river (dist).

### Data at samples sites
data(meuse)
flood    <- ifelse(meuse$ffreq==1, 1, 0 )# Binary response variable
coords   <- meuse[,c("x","y")]           # Coordinates
x        <- meuse[,"dist"]               # Covariate

### Data at prediction sites
data(meuse.grid)
coords0  <- meuse.grid[,c("x","y")]      # Coordinates
x0       <- meuse.grid[,"dist"]          # Covariate

flood is spatially plotted as follows:

obs_s    <- st_as_sf( data.frame(coords, flood), coords= c("x","y"), crs=28992)
plot(obs_s[,"flood"], pch = 20, key.pos=4, axes=TRUE)

Coarse-to-fine spatial GLMM (CF-GLM)

The code implementing the CF-GLM is essentially the same as in the previous example, except that family is set to binomial() for modeling binary data:

set.seed(1234) # For this vignette, training samples are fixed
mod_hv   <- cf_glm_hv(y = flood, x = x, coords = coords, family=binomial())

In the subsequent full model training using the cf_glm function, the covariates (x0) and coordinates (coords0) at the prediction sites are also specified to enable spatial prediction:

mod      <- cf_glm(y = flood, x=x, coords = coords, 
                   x0=x0, coords0 = coords0, mod_hv = mod_hv)

The estimated regression coefficients, standard deviations of the scalewise components, and error statistics are displayed as follows:

mod

Mapping

Predictive values

The predictive values and their standard deviations at the grid cells are mapped as follows:

### Convert gridded points to gridded polygons (for clear visualization)
meuse.grid_sp             <- meuse.grid
coordinates(meuse.grid_sp)<- c("x", "y")
gridded(meuse.grid_sp)    <- TRUE
meuse.grid_sf             <- st_as_sf(as(meuse.grid_sp, "SpatialPolygons"))
st_crs(meuse.grid_sf)     <- 28992

### Mapping predictive mean and standard deviations
meuse.grid_sf$pred        <- mod$pred0$pred   # Predictive mean
meuse.grid_sf$pred_sd     <- mod$pred0$pred_sd# Predictive standard deviations
plot(meuse.grid_sf[,"pred"], border = NA, nbreaks = 20, key.pos=4,axes=TRUE)
plot(meuse.grid_sf[,"pred_sd"], pal = function(n) hcl.colors(n, "Viridis"),
     border = NA,key.pos=4,axes=TRUE)

Prediction at new sites with predict()

As in cf_lm, the model can be fitted without prediction sites (no x0, coords0), and predict() then predicts at any sites later from the fitted object alone, with the same result as giving those sites as coords0 in cf_glm():

mod_f <- cf_glm(y = flood, x = x, coords = coords, mod_hv = mod_hv)   # no x0, coords0
p     <- predict(mod_f, x0 = x0, coords0 = coords0, probs = c(0.025, 0.5, 0.975),
                 se_type = "mean")

Applied to mod above, which avoids refitting here, it reproduces mod$pred0. For a binary response the interval of a single 0/1 observation is degenerate, so the quantiles of the predicted probability are asked for with se_type = "mean":

p <- predict(mod, x0 = x0, coords0 = coords0, probs = c(0.025, 0.5, 0.975),
             se_type = "mean")
head(p)

Multiscale spatial pattern/feature extraction

Let us extract the large- and small-scale processes, corresponding to bandwidth ranges of 1000+ and 0–1000, respectively.

mod_s1<- sp_scalewise(mod,bw_range=c(1000,Inf)) # Large scale (1000 <= bandwidth)
mod_s2<- sp_scalewise(mod,bw_range=c(0,1000))   # Small scale (0 <= bandwidth <= 1000)

The extracted scalewise processes are mapped as follows:

meuse.grid_sf$z1    <- mod_s1$pred0$pred
meuse.grid_sf$z2    <- mod_s2$pred0$pred
plot(meuse.grid_sf[,c("z1","z2")], border=NA, nbreaks=20, key.pos=4, axes=TRUE)

Interactive mapping

The maps above are static. spCFmap() opens a Shiny application that draws the same results over a basemap. Without arguments it launches the full application, whose bundled demo data include meuse and the Greater Glasgow health data used in the two examples above:

spCFmap()

Passing a fitted model maps that object directly, with crs naming the coordinate reference system of the coordinates given to the model. For the binary example above the coordinates are in the Dutch RD grid (EPSG:28992), and for the count example in the British National Grid (EPSG:27700):

spCFmap(mod, crs = 28992)   # the flood-probability model of Example 2

A model fitted without prediction sites -- as in Example 1, where cf_glm was called without x0/coords0 -- is mapped at its sample sites instead, so both examples can be explored this way.

Reference

spCF documentation built on Oct. 5, 2026, 5:07 p.m.