cf_glm: Coarse-to-fine spatial generalized linear mixed models...

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

Coarse-to-fine spatial generalized linear mixed models (CF-GLMMs)

Description

Scalable prediction, regression, and multiscale analysis via CF-GLMMs.

Usage

cf_glm(
  y,
  x = NULL,
  coords,
  offset = NULL,
  x0 = NULL,
  coords0 = NULL,
  offset0 = NULL,
  mod_hv,
  robust_se = TRUE,
  se_type = c("prediction", "mean"),
  se_method = c("opt", "classic"),
  keep_scales = TRUE
)

Arguments

y	Vector of response variables (N x 1), including continuous, count, and binary responses following an exponential family distribution.
x	Matrix of covariates (N x K).
coords	Matrix of 2-dimensional point coordinates (N x 2).
offset	Optional. Vector of offset variables (N x 1) included in the linear predictor, consistent with glm.
x0	Optional. Matrix of covariates at prediction sites (N0 x K).
coords0	Optional. Matrix of 2-dimensional point coordinates at prediction sites (N0 x 2).
offset0	Optional. Vector of offset variables at prediction sites (N0 x 1)
mod_hv	Output object of the cf_glm_hv function.
robust_se	If TRUE (default), coefficient standard errors and predictive uncertainty are computed using a cluster-robust sandwich estimator accounting for local spatial correlation. Set FALSE to use naive SEs (not recommended).
se_type	Type of predictive uncertainty in pred/pred_q. "prediction" (default) returns the OBSERVATION predictive for a new data point, holdout-calibrated on the cf_glm_hv validation samples (Gaussian: mean uncertainty + residual variance, split-conformal SD scaling; Poisson: negative-binomial count predictive; binomial: temperature -calibrated probability with pred_sd = sqrt(p(1-p))). Negative binomial (negbin): negative-binomial count predictive with the fitted dispersion. Other families use a moment-matched observation predictive with the holdout-estimated dispersion (Gamma: gamma; inverse.gaussian: inverse Gaussian; quasipoisson: negative binomial; quasibinomial: beta for proportions, as binomial for 0/1 data; otherwise normal), with the mean-uncertainty scale calibrated to 95% holdout coverage. The mean/signal versions are kept in pred_signal/pred_q_signal. "mean" returns the signal (mean) uncertainty only (previous behaviour). See other$calibration for the fitted calibration.
se_method	Cluster-robust coefficient-SE estimator (used when robust_se = TRUE). "opt" (default) splits the sandwich meat into a field-removed observation-noise part and a field part that adds the calibrated field variance back with a within-block exp(-d/h) correlation (h = median committed bandwidth); this is near-nominal. A refit-free leverage leave-one-out ceiling then caps the field term, preventing over-coverage for count (Poisson) responses while leaving already-calibrated families unchanged. "classic" keeps the realised field inside the working residual (the previous behaviour), which is valid but conservative.
keep_scales	If TRUE (default), the scale-wise processes Z, Z_sd, Z0 and Z0_sd are kept in the output. They hold one column per selected scale for every sample (and prediction) site, which makes them the largest part of the fitted object for large data. With FALSE they are dropped (NULL); predictions, standard errors, sd_summary and the maps of spCFmap for the total prediction are unchanged, but sp_scalewise needs them.

Details

The link-scale spatial-process predictive variance is bounded stage by stage as in cf_lm: \min(\tau\, pv_r, \kappa s_r^2) with the stage caps rescaled to sum to the marginal field variance, and the holdout factor \tau solving the working-weighted moment equation. For the binomial family the field variance is left uncapped.

Value

A list with the following elements:

beta

Regression coefficients, their standard errors, and the lower and upper limits of the 95 percent confidence intervals.

sd_summary

Standard deviation of the regression term (xb), spatial process (spatial_scale1, spatial_scale2,...), additional learning, and residuals.

e_summary

Holdout validation accuracy evaluated on the validation samples: R-squared (validation_Pseudo-R2), root mean squared error (validation_RMSE), and mean absolute error (validation_MAE).

pred

Predictive means and standard deviations (sample sites). The spatial-process contribution to the predictive SD is rescaled by a holdout-calibrated factor (stored as other$tau) estimated on the validation samples.

pred0

Predictive means and standard deviations (prediction sites).

pred_q

Predictive quantiles on the response scale at the sample sites. A data frame whose columns q0.005, q0.025, q0.05, q0.1, ..., q0.9, q0.95, q0.975, q0.995 give the corresponding quantile levels, obtained by Gaussian approximation on the link scale followed by inverse-link transformation (with se_type = "prediction", from the calibrated observation predictive). Not stored in the object: mod$pred_q computes it on access, at the 15 levels listed above; predict.cf_glm gives them at other levels and at new sites.

pred0_q

Predictive quantiles on the response scale at the prediction sites. Column structure is identical to pred_q. NULL when prediction sites are not supplied.

bands

Bandwidth values for each scale. The i-th bandwidth corresponds to the i-th column of the Z matrix.

Z

Predictive mean of the spatial process at each scale (sample sites; list).

Z_sd

Predictive standard deviation of the spatial process at each scale (sample sites; list).

Z0

Predictive mean of the spatial process at each scale (prediction sites; list).

Z0_sd

Predictive standard deviation of the spatial process at each scale (prediction sites; list). Z, Z_sd, Z0 and Z0_sd are NULL when keep_scales = FALSE.

other

Other internally used output objects.

Author(s)

Daisuke Murakami

References

Murakami, D., Comber, A., Yoshida, T., Tsutsumida, N., Brunsdon, C., & Nakaya, T. (2025). Coarse-to-fine spatial GLMMs for scalable prediction and multiscale analysis. *ArXiv preprint*, 2605.01157. https://doi.org/10.48550/arXiv.2605.01157

See Also

cf_glm_hv, sp_scalewise

Examples

################ Example 1: Count data modeling/Disease mapping/smoothing
set.seed(1234)
require( CARBayesdata )
require( sf )
data(pollutionhealthdata)
data(GGHB.IZ)

### Data
dat      <- pollutionhealthdata[pollutionhealthdata$year==2011,]
y        <- dat[,"observed"]             # count data
x        <- dat[,c("pm10","jsa","price")]
offset   <- log(dat[,"expected"])
coords   <- st_coordinates(st_centroid(GGHB.IZ))

### Holdout validation optimizing the number of spatial scales
mod_hv   <- cf_glm_hv(y = y, x = x, offset=offset, coords = coords, family=poisson())

### Spatial modeling and prediction
mod      <- cf_glm(y = y, x = x, coords = coords, mod_hv = mod_hv)
mod

### Mapping predictive mean and standard deviations (SD)
GGHB.IZ$y      <- y
GGHB.IZ$pred   <- mod$pred$pred
GGHB.IZ$pred_sd<- mod$pred$pred_sd
plot(GGHB.IZ[,c("pred")],lwd=0.2,axes=TRUE, key.pos=4,nbreaks=50)   # Predictive mean
plot(GGHB.IZ[,c("pred_sd")],lwd=0.2,axes=TRUE, key.pos=4,nbreaks=50)# Predictive SD

### Multiscale spatial pattern/feature extraction
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)
GGHB.IZ$z1  <- mod_s1$pred$pred
GGHB.IZ$z2  <- mod_s2$pred$pred
plot(GGHB.IZ[,c("z1","z2")],lwd=0.2,axes=TRUE,key.pos=4, nbreaks=50)# Extracted features



################ Example 2: Binary data modeling/spatial prediction
set.seed(1234)
require(sp); require(sf)
data(meuse)
data(meuse.grid)

### Data
y        <- ifelse(meuse$ffreq==1, 1, 0 )# binary data
coords   <- meuse[,c("x","y")]
x        <- meuse[,"dist"]

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

### Holdout validation optimizing the number of spatial scales
mod_hv   <- cf_glm_hv(y = y, x = x, coords = coords, family=binomial())

### Spatial modeling and prediction
mod      <- cf_glm(y = y, x=x, coords = coords, x0=x0, coords0 = coords0,
                   mod_hv = mod_hv)
mod

### Mapping predictive mean and standard deviations (SD)
meuse.grid$pred   <- mod$pred0$pred
meuse.grid$pred_sd<- mod$pred0$pred_sd
meuse.grid_sf     <- st_as_sf(meuse.grid, coords = c("x","y"))
plot(meuse.grid_sf[,"pred"], pch = 15, cex = 0.8, nbreaks = 20)   # Predictive mean
plot(meuse.grid_sf[,"pred_sd"], pch = 15, cex = 0.8, nbreaks = 20)# Predictive SD

### Multiscale spatial pattern/feature extraction
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)
meuse.grid_sf$z1    <- mod_s1$pred0$pred
meuse.grid_sf$z2    <- mod_s2$pred0$pred
plot(meuse.grid_sf[,c("z1","z2")], pch = 15,
     cex = 0.5, nbreaks = 20,axes=TRUE) # Predictive means

### The same fit, explored interactively over a basemap
# spCFmap(mod, crs = 28992)   # crs = the system the coordinates are in

### Prediction with predict(): the model can be fitted WITHOUT prediction
### sites (no x0, coords0) and used to predict at any sites later; the
### training data are not needed then. For a binary response,
### se_type = "mean" gives the quantiles of the probability (those of a
### single 0/1 observation are degenerate).
mod_f    <- cf_glm(y = y, x = x, coords = coords, mod_hv = mod_hv)  # no x0, coords0
p        <- predict(mod_f, x0 = x0, coords0 = coords0,
                    probs = c(0.025, 0.975), se_type = "mean")
head(p)
all.equal(p$pred, mod$pred0_signal$pred)   # same as cf_glm(..., coords0 = coords0)

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