Package {sfclust}

Title:	Bayesian Spatial Functional Clustering
Version:	1.1.0
Description:	Bayesian clustering of spatial regions with similar functional shapes using spanning trees and latent Gaussian models. The method enforces spatial contiguity within clusters and supports a wide range of latent Gaussian models, including non-Gaussian likelihoods, via the R-INLA framework. The algorithm is based on Zhong, R., Chacón-Montalván, E. A., and Moraga, P. (2026) <doi:10.1002/sim.70597>, extending the approach of Zhang, B., Sang, H., Luo, Z. T., and Huang, H. (2023) <doi:10.1214/22-AOAS1643>. The package includes tools for model fitting, convergence diagnostics, visualization, and summarization of clustering results.
License:	MIT + file LICENSE
URL:	https://erickchacon.github.io/sfclust/, https://github.com/ErickChacon/sfclust
BugReports:	https://github.com/ErickChacon/sfclust/issues
Encoding:	UTF-8
Imports:	igraph, sf, SparseM, stars, dplyr, methods, Matrix, ggplot2, patchwork
Suggests:	ggraph, here, INLA, knitr, rmarkdown, testthat (≥ 3.0.0)
Additional_repositories:	https://inla.r-inla-download.org/R/stable
Depends:	R (≥ 4.1.0)
Config/testthat/edition:	3
VignetteBuilder:	knitr
LazyData:	true
Config/Needs/website:	rmarkdown, fda, maps
Config/roxygen2/version:	8.0.0
NeedsCompilation:	no
Packaged:	2026-07-28 15:23:11 UTC; rstudio
Author:	Erick A. Chacón-Montalván [aut, cre], Ruiman Zhong [aut], Paula Moraga [aut]
Maintainer:	Erick A. Chacón-Montalván <erick.chaconmontalvan@wur.nl>
Repository:	CRAN
Date/Publication:	2026-07-28 22:40:03 UTC

sfclust: Bayesian Spatial Functional Clustering

Description

The sfclust package implements the clustering algorithm proposed in "Bayesian Spatial Functional Data Clustering: Applications in Disease Surveillance", available at doi:10.1002/sim.70597. The package provides tools for performing Bayesian spatial functional clustering, as well as methods for diagnostic analysis, visualization, and summarization of results.

Author(s)

Ruiman Zhong ruiman.zhong@kaust.edu.sa, Erick A. Chacón-Montalván erick.chaconmontalvan@wur.nl, Paula Moraga paula.moraga@kaust.edu.sa

See Also

Useful links:

• https://erickchacon.github.io/sfclust/

• https://github.com/ErickChacon/sfclust

• Report bugs at https://github.com/ErickChacon/sfclust/issues

Spatio-temporal NDWI2 raster data from El Chaparrillo

Description

A stars raster object containing the Normalized Difference Water Index 2 (NDWI2) derived from Sentinel-2 satellite imagery over El Chaparrillo, an agricultural research center in Spain, across 86 time points.

Usage

data(chapa)

Format

A stars object with:

ndwi2

NDWI2 values (numeric) derived from Sentinel-2 bands. Non-NA values define the study area mask.

dimensions

Three dimensions: x (26 columns), y (27 rows), time (86 dates from 2025-05-15 to 2025-12-01). CRS: WGS 84 / UTM zone 30N.

Examples

library(sfclust)

data(chapa)
chapa
plot(chapa["ndwi2"])

Prepare data in long format

Description

Convert spatio-temporal data to long format with a flat spatial index (ids) and ordered observation indices (⁠id_<dimname>⁠ for each non-spatial dimension). This is a pure converter: all rows are returned including cells that are all-NA. No filtering is applied. Use filter_df() downstream to restrict to valid cells.

Usage

data_all(x, spnames = "geometry")

Arguments

Value

A long-format data frame with columns id (flat array position in the original stars object, column-major over all dimensions), ids (flat spatial index, column-major over spnames, spnames[1] varies fastest), ⁠id_<dimname>⁠ for each functional dimension, and all variables from x. All rows are returned; ids is not remapped to ⁠1..n_valid⁠.

Examples

library(sfclust)
library(stars)

dims <- st_dimensions(
  geometry = st_sfc(lapply(1:5, function(i) st_point(c(i, i)))),
  time = seq(as.Date("2024-01-01"), by = "1 day", length.out = 3)
)
stdata <- st_as_stars(cases = array(1:15, dim = c(5, 3)), dimensions = dims)

data_all(stdata)

Filter a long-format data frame to valid spatial cells

Description

Keeps only rows whose spatial index (ids) appears in valid_ids and adds a sid column (1..n_valid) giving each valid cell's rank among the valid cells. sid matches the indexing of membership vectors returned by genclust(). If valid_ids is NULL, all rows are kept and sid is assigned by rank of ids.

Usage

filter_df(df, valid_ids = NULL)

Arguments

Value

A filtered data frame with rows corresponding to valid spatial cells only, with an additional sid column (integer, 1..n_valid).

Fitted Values for sfclust Objects

Description

This function calculates the fitted values for a specific clustering sample in an sfclust object, based on the estimated models for each cluster. The fitted values are computed using the membership assignments and model parameters associated with the selected clustering sample.

Usage

## S3 method for class 'sfclust'
fitted(object, sample = object$clust$id, sort = FALSE, aggregate = FALSE, ...)

Arguments

Value

A data frame with fitted values and cluster assignments, keyed by id. For sfclust_stars objects, a stars object is returned instead.

Examples

if (requireNamespace("INLA", quietly = TRUE)) {
library(sfclust)

data(stgaus)
result <- sfclust(stgaus, formula = y ~ f(id_time, model = "rw1"), niter = 10,
  nmessage = 1)

# Estimated values ordering clusters by size
df_est <- fitted(result, sort = TRUE)

# Estimated values aggregated by cluster
df_est <- fitted(result, aggregate = TRUE)

# Estimated values using a particular clustering sample
df_est <- fitted(result, sample = 3)
}

Generate initial cluster assignments

Description

Creates an undirected graph from spatial data or a weighted adjacency matrix, computes its minimum spanning tree (MST), and partitions it into nclust clusters. Accepts weighted matrix or Matrix objects, and stars objects with either vector geometry (contiguity via sf::st_touches()) or raster dimensions (4-connected grid adjacency for non-NA pixels).

Usage

genclust(x, ...)

## Default S3 method:
genclust(x, ...)

## S3 method for class 'matrix'
genclust(x, nclust = 10, weights = NULL, ...)

## S3 method for class 'Matrix'
genclust(x, nclust = 10, weights = NULL, ...)

## S3 method for class 'stars'
genclust(x, nclust = 10, spnames = NULL, response = NULL, weights = NULL, ...)

Arguments

Value

A list with:

• graph: undirected igraph object representing spatial contiguity.

• mst: minimum spanning tree of graph.

• membership: integer vector of cluster assignments (length = number of valid spatial units).

• valid_ids: integer vector of flat spatial positions included in the graph (only present for stars input; NULL for matrix/Matrix input).

Examples

library(sfclust)
library(sf)
library(stars)

# stars object with vector geometry
geom <- st_make_grid(cellsize = c(1, 1), offset = c(0, 0), n = c(3, 2))
x <- st_as_stars(st_sf(z = 1:6, geometry = geom))
clust <- genclust(x, nclust = 3)
plot(st_sf(geom, cluster = factor(clust$membership)))

# stars raster input
x <- st_as_stars(cluster = matrix(1:35, 5))
clust <- genclust(x, nclust = 4)
x$cluster <- clust$membership
plot(x, col = rainbow(4))

# matrix input
x <- matrix(c(0,1,0,1, 1,0,1,0, 0,1,0,1, 1,0,1,0), nrow = 4)
clust <- genclust(x, nclust = 2)
clust$membership

Fit models and compute the log marginal likelihood for all clusters

Description

Fit the specified INLA model to each cluster and compute the log marginal likelihood for each cluster specified in the membership vector.

Usage

log_mlik_all(
  membership,
  data,
  correction = TRUE,
  detailed = FALSE,
  inla_args = NULL
)

Arguments

Value

A numeric vector containing the log marginal likelihood for each cluster or the the fitted INLA model for each cluster when detailed = TRUE.

Plot function for sfclust objects

Description

Visualizes fitted cluster functions (plot 1) and a log marginal likelihood traceplot (plot 2). For sfclust_stars objects, also includes a spatial map of cluster assignments (plot 1 in that method, shifting the others to 2 and 3).

Usage

## S3 method for class 'sfclust'
plot(
  x,
  sample = x$clust$id,
  which = 1:2,
  clusters = NULL,
  sort = FALSE,
  legend = FALSE,
  fnames = NULL,
  ...
)

Arguments

Value

A composed patchwork object.

Plot functional shapes of cluster linear predictors

Description

Plots the estimated mean functional shape (linear predictor or inverse-link scale) for each cluster in a given MCMC sample of an sfclust object.

Usage

plot_clusters_fitted(
  x,
  sample = x$clust$id,
  clusters = NULL,
  sort = FALSE,
  legend = FALSE,
  inv_link = TRUE,
  fnames = NULL,
  ...
)

Arguments

Value

A ggplot2 object.

Plot a spatial map of cluster assignments

Description

Produces a ggplot2 map of spatial regions colored by their cluster assignment for a given MCMC sample of an sfclust_stars object.

Usage

plot_clusters_map(
  x,
  sample = x$clust$id,
  clusters = NULL,
  sort = FALSE,
  legend = FALSE,
  geom_before = NULL,
  ...
)

Arguments

Value

A ggplot2 object.

Plot observed time series by cluster

Description

Plots individual region time series faceted by cluster, overlaid with the cluster mean, for a given variable in an sfclust object.

Usage

plot_clusters_series(x, var, clusters = NULL, sort = FALSE, fnames = NULL, ...)

Arguments

Value

A ggplot2 object with one facet per cluster.

Plot log marginal likelihood convergence trace

Description

Plots the log marginal likelihood across MCMC samples for an sfclust object, highlighting the selected sample.

Usage

plot_log_mlik(x, sample = x$clust$id, ...)

Arguments

Value

A ggplot2 object.

Print method for sfclust objects

Description

Prints details of an sfclust object, including the (i) within-cluster formula; (ii) hyperparameters used for the MCMC sample such as the number of clusters penalty (q) and the movement probabilities (move_prob); (iii) the number of movement type counts during the MCMC sampling; and (iv) the log marginal likelihood of the model of the last clustering sample.

Usage

## S3 method for class 'sfclust'
print(x, ...)

Arguments

Value

Invisibly returns the input sfclust object x. The function also prints a summary of:

• the within-cluster model formula,

• clustering hyperparameters,

• movement counts from the MCMC sampler,

• and the log marginal likelihood of the selected sample.

Bayesian spatial functional clustering

Description

sfclust() is the main user-facing function for Bayesian spatial functional clustering via reversible-jump MCMC. It dispatches on the class of the first argument:

Usage

sfclust(x, ...)

## Default S3 method:
sfclust(x, ...)

## S3 method for class 'data.frame'
sfclust(
  x,
  adjacency,
  graphdata = NULL,
  fnames = NULL,
  nclust = 10,
  move_prob = c(0.425, 0.425, 0.1, 0.05),
  logpen = log(1 - 0.5),
  correction = TRUE,
  niter = 100,
  burnin = 0,
  thin = 1,
  nmessage = 10,
  path_save = NULL,
  nsave = nmessage,
  ...
)

## S3 method for class 'stars'
sfclust(
  x,
  nclust = 10,
  graphdata = NULL,
  spnames = NULL,
  move_prob = c(0.425, 0.425, 0.1, 0.05),
  logpen = log(1 - 0.5),
  correction = TRUE,
  niter = 100,
  burnin = 0,
  thin = 1,
  nmessage = 10,
  path_save = NULL,
  nsave = nmessage,
  ...
)

Arguments

Details

• sfclust.data.frame(): core interface — takes a pre-built long-format data frame and a weighted adjacency matrix. Use this when working with any data format after converting it yourself.

• sfclust.stars(): stars wrapper — takes a stars spatio-temporal object, converts it to long format, builds the spatial graph, and calls the core algorithm.

This implementation draws inspiration from the methods described in the paper: "Bayesian Clustering of Spatial Functional Data with Application to a Human Mobility Study During COVID-19" by Bohai Zhang, Huiyan Sang, Zhao Tang Luo, and Hui Huang, published in The Annals of Applied Statistics, 2023. For further details on the methodology, please refer to:

• The paper: doi:10.1214/22-AOAS1643

• Supplementary material: doi:10.1214/22-AOAS1643SUPPB

The MCMC algorithm in this implementation is largely based on the supplementary material provided in the paper. However, we have generalized the computation of the marginal likelihood ratio by leveraging INLA (Integrated Nested Laplace Approximation). This generalization enables integration over all parameters and hyperparameters, allowing for inference within a broader family of distribution functions and model terms, thereby extending the scope and flexibility of the original approach. Further details of our approach can be found in our paper "Bayesian spatial functional data clustering: applications in disease surveillance" by Ruiman Zhong, Erick A. Chacón-Montalván, Paula Moraga:

• The paper: doi:10.1002/sim.70597

Value

An sfclust object (from sfclust.data.frame) or an sfclust_stars object inheriting from sfclust (from sfclust.stars). Both contain:

• samples: MCMC trace with membership, log_mlike, and move_counts.

• clust: selected clustering with id, membership, and fitted models.

sfclust_stars additionally carries input_args with stars (structural shell of the input), spnames, and fnames, used by spatial plot methods.

Author(s)

Ruiman Zhong ruiman.zhong@kaust.edu.sa, Erick A. Chacón-Montalván erick.chaconmontalvan@wur.nl, Paula Moraga paula.moraga@kaust.edu.sa

Examples

if (requireNamespace("INLA", quietly = TRUE)) {
library(sfclust)

# Stars interface: Gaussian model
data(stgaus)
result <- sfclust(stgaus, formula = y ~ f(id_time, model = "rw1"),
  niter = 10, nmessage = 1)

print(result)
summary(result, sort = TRUE)
fitted(result, sort = TRUE)

plot(result)
plot_clusters_series(result, var = y)

result2 <- update(result, niter = 2, nmessage = 1)
plot(result2)

# Stars interface: Binomial model
data(stbinom)
result <- sfclust(stbinom, formula = cases ~ poly(id_time, 2) + f(id),
  family = "binomial", Ntrials = population, niter = 10, nmessage = 1)

print(result)
summary(result, sort = TRUE)
fitted(result, sort = TRUE)

plot(result)
plot_clusters_series(result, var = cases/population)

result2 <- update(result, niter = 2, nmessage = 1)
plot(result2)

# data.frame interface: Poisson model
 ns <- 6L; nt <- 4L
 set.seed(4)
 df  <- data.frame(
   id       = seq_len(ns * nt),
   ids      = rep(seq_len(ns), nt),
   id_time  = rep(seq_len(nt), each = ns),
   expected = rep(10L, ns * nt)
 )
 df <- transform(df,
   y = rpois(ns * nt, expected * exp(0.5 * id_time * rep(c(-1, 1), each = ns / 2, nt)))
 )
 adj <- Matrix::sparseMatrix(
   i = c(1, 2, 4, 5, 1, 2, 3), j = c(2, 3, 5, 6, 4, 5, 6),
   x = 1L, dims = c(ns, ns), symmetric = TRUE
 )

 result <- sfclust(df, adjacency = adj, nclust = 3, fnames = "id_time",
   formula = y ~ 1 + id_time, family = "poisson", E = expected,
   niter = 3, burnin = 0, thin = 1, nmessage = 1)

 print(result)
 summary(result, sort = TRUE)
 fitted(result, sort = TRUE)

 plot(result)
 plot_clusters_series(result, var = y)

 result2 <- update(result, niter = 2, nmessage = 1)
 plot(result2)
}

Spatio-temporal Binomial data

Description

A simulated stars object containing binomial response data with a functional clustering pattern defined by polynomial fixed effects. This dataset includes the variables cases and population observed across 100 simulated spatial regions over 91 time points.

Usage

data(stbinom)

Format

A stars object with:

cases

Number of observed cases (integer)

population

Population at risk (integer)

dimensions

Two dimensions: geometry (spatial features) and time (daily observations)

Examples

library(sfclust)

data(stbinom)
stbinom
plot(stbinom["cases"])

Spatio-temporal Gaussian data

Description

A simulated stars object containing Gaussian response data with a functional clustering pattern uging random walk processes. This dataset includes the response variable y observed across 100 simulated spatial regions over 91 time points.

Usage

data(stgaus)

Format

A stars object with:

y

Response variable

Examples

library(sfclust)

data(stgaus)
stgaus
plot(stgaus["y"])

Summary method for sfclust objects

Description

This function summarizes the cluster assignments from the desired clustering sample.

Usage

## S3 method for class 'sfclust'
summary(object, sample = object$clust$id, sort = FALSE, ...)

Arguments

Value

Invisibly returns a table with the number of regions in each cluster for the selected sample. The function also prints a summary that includes:

• the within-cluster model formula,

• the total number of MCMC clustering samples,

• the cluster membership counts for the specified sample (optionally sorted),

• and the log marginal likelihood of the selected clustering sample.

Update MCMC Clustering Procedure

Description

This function continues the MCMC sampling of a sfclust object based on previous results or update the model fitting for a specified sample clustering if the argument sample is provided.

Usage

## S3 method for class 'sfclust'
update(
  object,
  niter = NULL,
  burnin = 0,
  thin = 1,
  nmessage = 10,
  sample = NULL,
  path_save = NULL,
  nsave = nmessage,
  ...
)

Arguments

Details

This function takes the last state of the Markov chain from a previous sfclust_fit / sfclust execution and uses it as the starting point for additional MCMC iterations. If niter = NULL (the default), no additional iterations are run; instead, the within-cluster INLA models are fitted for the clustering given by sample (or the current clust$id if sample is NULL). This is useful when loading a checkpoint saved during a run, where the models were not yet stored: result <- update(result).

Value

An updated sfclust object with (i) new clustering samples if sample is not specified, or (ii) updated within-cluster model results if sample is given.