Package 'spStack'

spStack: Bayesian Geostatistics Using Predictive Stacking

Description

This package delivers functions to fit Bayesian hierarchical spatial process models for point-referenced Gaussian, Poisson, binomial, and binary data using stacking of predictive densities. It involves sampling from analytically available posterior distributions conditional upon some candidate values of the spatial process parameters for both Gaussian response model as well as non-Gaussian responses, and, subsequently assimilate inference from these individual posterior distributions using Bayesian predictive stacking. Our algorithm is highly parallelizable and hence, much faster than traditional Markov chain Monte Carlo algorithms while delivering competitive predictive performance.

In context of inference for spatial point-referenced data, Bayesian hierarchical models involve latent spatial processes characterized by spatial process parameters, which besides lacking substantive relevance in scientific contexts, are also weakly identified and hence, impedes convergence of MCMC algorithms. This motivates us to build methodology that involves fast sampling from posterior distributions conditioned on a grid of the weakly identified model parameters and combine the inference by stacking of predictive densities (Yao et. al 2018). We exploit the Bayesian conjugate linear modeling framework for the Gaussian case (Zhang, Tang and Banerjee 2025) and the generalized conjugate multivariate distribution theory (Pan, Zhang, Bradley and Banerjee 2025) to analytically derive the individual posterior distributions.

Details

Package:	spStack
Type:	Package
Version:	1.1.0
License:	GPL-3

Accepts a formula, e.g., y~x1+x2, for most regression models accompanied by candidate values of spatial process parameters, and returns posterior samples of the regression coefficients and the latent spatial random effects. Posterior inference or prediction of any quantity of interest proceed from these samples. Main functions are -
spLMexact()
spGLMexact()
spLMstack()
spGLMstack()

Author(s)

Maintainer: Soumyakanti Pan [email protected] (ORCID)

Authors:

• Soumyakanti Pan [email protected] (ORCID)

• Sudipto Banerjee [email protected]

References

Zhang L, Tang W, Banerjee S (2025). "Bayesian Geostatistics Using Predictive Stacking." Journal of the American Statistical Association, In press. doi:10.1080/01621459.2025.2566449.

Pan S, Zhang L, Bradley JR, Banerjee S (2025). "Bayesian Inference for Spatial-temporal Non-Gaussian Data Using Predictive Stacking." Bayesian Analysis, In Press. doi:10.1214/25-BA1582.

Yao Y, Vehtari A, Simpson D, Gelman A (2018). "Using Stacking to Average Bayesian Predictive Distributions (with Discussion)." Bayesian Analysis, 13(3), 917-1007. doi:10.1214/17-BA1091.

See Also

Useful links:

• https://span-18.github.io/spStack-dev/

• Report bugs at https://github.com/SPan-18/spStack-dev/issues

---

Create a collection of candidate models for stacking

Description

Creates an object of class 'candidateModels' that contains a list of candidate models for stacking. The function takes a list of candidate values for each model parameter and returns a list of possible combinations of these values based on either simple aggregation or Cartesian product of indivdual candidate values.

Usage

candidateModels(params_list, aggregation = "simple")

Arguments

params_list	a list of candidate values for each model parameter. See examples for details.
aggregation	a character string specifying the type of aggregation to be used. Options are 'simple' and 'cartesian'. Default is 'simple'.

Value

an object of class 'candidateModels'

Author(s)

Soumyakanti Pan [email protected],
Sudipto Banerjee [email protected]

See Also

stvcGLMstack()

Examples

m1 <- candidateModels(list(phi_s = c(1, 1), phi_t = c(1, 2)), "simple")
m1
m2 <- candidateModels(list(phi_s = c(1, 1), phi_t = c(1, 2)), "cartesian")
m2
m3 <- candidateModels(list(phi_s = list(c(1, 1), c(1, 2)),
                          phi_t = list(c(1, 3), c(2, 3)),
                          boundary = c(0.5, 0.75)),
                      "simple")

---

Different Cholesky factor updates

Description

Provides functions that implements different types of updates of a Cholesky factor that includes rank-one update, single row/column deletion update and a block deletion update.

Usage

cholUpdateRankOne(A, v, alpha, beta, lower = TRUE)

cholUpdateDel(A, del.index, lower = TRUE)

cholUpdateDelBlock(A, del.start, del.end, lower = TRUE)

Arguments

A	an $n\times n$ triangular matrix
v	an $n\times 1$ matrix/vector
alpha	scalar; if not supplied, default is 1
beta	scalar; if not supplied, default is 1
lower	logical; if A is lower-triangular or not
del.index	an integer from 1 to $n$ indicating the row/column to be deleted
del.start	an integer from 1 to $n$ indicating the first row/column of a block to be deleted, must be at least 1 less than del.end
del.end	an integer from 1 to $n$ indicating the last row/column of a block to be deleted, must be at least 1 more than del.start

Details

Suppose $B = AA^\top$ is a $n \times n$ matrix with $A$ being its lower-triangular Cholesky factor. Then rank-one update corresponds to finding the Cholesky factor of the matrix $C = \alpha B + \beta vv^\top$ for some $\alpha,\beta\in\mathbb{R}$ given $A$ (see, Krause and Igel 2015). Similarly, single row/column deletion update corresponds to finding the Cholesky factor of the $(n-1)\times(n-1)$ matrix $B_i$ which is obtained by removing the $i$-th row and column of $B$, given $A$ for some $i - 1, \ldots, n$. Lastly, block deletion corresponds to finding the Cholesky factor of the $(n-n_k)\times(n-n_k)$ matrix $B_{I}$ for a subset $I$ of $\{1, \ldots, n\}$ containing $n_k$ consecutive indices, given the factor $A$.

Value

An $m \times m$ lower-triangular matrix with $m = n$ in case of cholUpdateRankOne(), $m = n - 1$ in case of cholUpdateDel(), and, $m = n - n_k$ in case of cholUpdateDelBlock() where $n_k$ is the size of the block removed.

Author(s)

Soumyakanti Pan [email protected],
Sudipto Banerjee [email protected]

References

Oswin Krause and Christian Igel. 2015. "A More Efficient Rank-one Covariance Matrix Update for Evolution Strategies". In Proceedings of the 2015 ACM Conference on Foundations of Genetic Algorithms XIII (FOGA '15). Association for Computing Machinery, New York, NY, USA, 129-136. doi:10.1145/2725494.2725496.

Examples

n <- 10
A <- matrix(rnorm(n^2), n, n)
A <- crossprod(A)
cholA <- chol(A)

## Rank-1 update
v <- 1:n
APlusvvT <- A + tcrossprod(v)
cholA1 <- t(chol(APlusvvT))
cholA2 <- cholUpdateRankOne(cholA, v, lower = FALSE)
print(all(abs(cholA1 - cholA2) < 1E-9))

## Single Row-deletion update
ind <- 2
A1 <- A[-ind, -ind]
cholA1 <- t(chol(A1))
cholA2 <- cholUpdateDel(cholA, del.index = ind, lower = FALSE)
print(all(abs(cholA1 - cholA2) < 1E-9))

## Block-deletion update
start_ind <- 2
end_ind <- 6
del_ind <- c(start_ind:end_ind)
A1 <- A[-del_ind, -del_ind]
cholA1 <- t(chol(A1))
cholA2 <- cholUpdateDelBlock(cholA, start_ind, end_ind, lower = FALSE)
print(all(abs(cholA1 - cholA2) < 1E-9))

---

Optimal stacking weights

Description

Obtains optimal stacking weights given leave-one-out predictive densities for each candidate model.

Usage

get_stacking_weights(log_loopd, solver = NULL, verbose = TRUE)

Arguments

log_loopd	an $n \times M$ matrix with $i$-th row containing the leave-one-out predictive densities for the $i$-th data point for the $M$ candidate models.
solver	specifies the solver to use for obtaining optimal weights. Default is "CLARABEL". Internally calls CVXR::psolve().
verbose	if TRUE, prints output of optimization routine.

Details

The weights maximize the log score of the stacked leave-one-out predictive densities (Yao et al. 2018) over the simplex, using the CVXR solvers in the order given above. If none of them reaches an optimal solution, loo::stacking_weights() is used as a fallback when the package loo is installed; otherwise the weights are returned as NA with status "failed".

Value

A list with elements:

weights

optimal stacking weights as a numeric vector of length $M$ (NA if no solver succeeded, see Details).

status

solver status, returns "optimal" if solver succeeded, and "failed" if no solver succeeded.

solver

name of the solver used ("none" if no solver succeeded).

details

a list with the installed CVXR solvers (installed), the requested solver(s) (requested) and those of them not installed (missing.requested), the order in which the solvers were tried (search.order), a data frame of the attempts with the status or error of each solver (attempts), and whether the fallback loo::stacking_weights() was used (fallback).

Author(s)

Soumyakanti Pan [email protected],
Sudipto Banerjee [email protected]

References

Yao Y, Vehtari A, Simpson D, Gelman A (2018). "Using Stacking to Average Bayesian Predictive Distributions (with Discussion)." Bayesian Analysis, 13(3), 917-1007. doi:10.1214/17-BA1091.

See Also

CVXR::psolve(), spLMstack(), spGLMstack()

Examples

set.seed(1234)
data(simSpatial)
dat <- simSpatial[1:100, ]

cand.mod <- candidateModels(list(phi = c(2, 4, 6),
                                 nu = c(1, 1.5, 1.75),
                                 noise_sp_ratio = c(0.1)), "cartesian")

mod1 <- spLMstack(y_gauss ~ x1 + x2, data = dat,
                  coords = as.matrix(dat[, c("s1", "s2")]),
                  cor.fn = "matern",
                  candidate.models = cand.mod,
                  n.samples = 1000, loopd.method = "exact",
                  parallel = FALSE, verbose = TRUE)

loopd_mat <- do.call('cbind', mod1$loopd)
w_hat <- get_stacking_weights(loopd_mat)
print(round(w_hat$weights, 4))
print(w_hat$solver)
print(w_hat$status)

---

Calculate distance matrix

Description

Computes the inter-site Euclidean distance matrix for one or two sets of points.

Usage

iDist(coords.1, coords.2, ...)

Arguments

coords.1	an $n\times p$ matrix with each row corresponding to a point in $p$-dimensional space.
coords.2	an $m\times p$ matrix with each row corresponding to a point in $p$ dimensional space. If this is missing then coords.1 is used.
...	currently no additional arguments.

Value

The $n\times n$ or $n\times m$ inter-site Euclidean distance matrix.

Author(s)

Soumyakanti Pan [email protected],
Sudipto Banerjee [email protected]

Examples

n <- 10
p1 <- cbind(runif(n),runif(n))
m <- 5
p2 <- cbind(runif(m),runif(m))
D <- iDist(p1, p2)

---

Prediction of latent process at new spatial or temporal locations

Description

A function to sample from the posterior predictive distribution of the latent spatial or spatial-temporal process.

Usage

posteriorPredict(mod_out, coords_new, covars_new, joint = FALSE, nBinom_new)

Arguments

mod_out	an object returned by any model fit under fixed hyperparameters or using predictive stacking, i.e., spLMexact(), spLMstack(), spGLMexact(), spGLMstack(), svcLMexact(), svcLMstack(), svcGLMexact(), svcGLMstack(), stvcGLMexact(), stvcGLMstack(), stvcLMexact() or stvcLMstack().
coords_new	new spatial coordinates (an $n_{new} \times 2$ matrix) or, for the spatial-temporal models, a list with tags sp (an $n_{new} \times 2$ matrix) and time (a vector or one-column matrix) at which the latent process, the mean, and the response is to be predicted.
covars_new	new covariates at the new coordinates: a matrix or, for the varying coefficients models (svc and stvc), a list with tags fixed (covariates with fixed effects) and vc (covariates with spatially or spatially-temporally varying coefficients). See examples for the structure of this list.
joint	a logical value indicating whether to return the joint posterior predictive samples of the latent process at the new locations or times. Defaults to FALSE.
nBinom_new	a vector of the number of trials for each new prediction location or time. Only required if the model family is "binomial". Defaults to a vector of ones, indicating one trial for each new prediction.

Value

A modified object with the class name preceeded by the identifier pp separated by a .. For example, if input is of class spLMstack, then the output of this prediction function would be pp.spLMstack. The entry with the tag samples is updated and will include samples from the posterior predictive distribution of the latent process, the mean, and the response at the new locations or times. An entry with the tag prediction is added and contains the new coordinates and covariates, and whether the joint posterior predictive samples were requested.

Author(s)

Soumyakanti Pan [email protected],
Sudipto Banerjee [email protected]

See Also

spLMexact(), spLMstack(), spGLMexact(), spGLMstack(), svcLMexact(), svcLMstack(), svcGLMexact(), svcGLMstack(), stvcGLMexact(), stvcGLMstack(), stvcLMexact(), stvcLMstack()

Examples

set.seed(1234)
# training and test data sizes
n_train <- 100
n_pred <- 10

# Example 1: Spatial linear model
# split the data into training and prediction sets
data(simSpatial)
dat_train <- simSpatial[1:n_train, ]
dat_pred <- simSpatial[n_train + 1:n_pred, ]

# fit a spatial linear model using predictive stacking
cand.mod <- candidateModels(list(phi = c(2, 4, 6), nu = c(1, 1.5, 1.75),
                                 noise_sp_ratio = c(0.1, 0.25)), "cartesian")
mod1 <- spLMstack(y_gauss ~ x1 + x2, data = dat_train,
                  coords = as.matrix(dat_train[, c("s1", "s2")]),
                  cor.fn = "matern",
                  candidate.models = cand.mod,
                  n.samples = 1000, loopd.method = "exact",
                  parallel = FALSE, verbose = TRUE)

# prepare new coordinates and covariates for prediction
sp_pred <- as.matrix(dat_pred[, c("s1", "s2")])
X_new <- cbind(1, dat_pred$x1, dat_pred$x2)

# carry out posterior prediction
mod.pred <- posteriorPredict(mod1, coords_new = sp_pred, covars_new = X_new,
                             joint = TRUE)

# sample from the stacked posterior and posterior predictive distribution
post_samps <- stackedSampler(mod.pred)

# compare the predicted spatial effects with the truth
z_pred <- apply(post_samps$z.pred, 1, median)
cor(z_pred, dat_pred$z_true)

---

Recover posterior samples of scale parameters of spatial/spatial-temporal generalized linear models

Description

A function to recover posterior samples of scale parameters that were marginalized out during model fit. This is only applicable for spatial or, spatial-temporal generalized linear models. This function applies on outputs of functions that fits a spatial/spatial-temporal generalized linear model, such as spGLMexact(), spGLMstack(), svcGLMexact(), svcGLMstack(), stvcGLMexact(), and stvcGLMstack().

Usage

recoverGLMscale(mod_out)

Arguments

mod_out

an object returned by a fitting a spatial or spatial-temporal GLM.

Value

An object of the same class as input, and updates the list tagged samples with the posterior samples of the scale parameters. The new tags are sigmasq.beta and z.scale.

Author(s)

Soumyakanti Pan [email protected],
Sudipto Banerjee [email protected]

See Also

spGLMexact(), spGLMstack(), svcGLMexact(), svcGLMstack(), stvcGLMexact(), stvcGLMstack()

Examples

set.seed(1234)
data(simSpatial)
dat <- simSpatial[1:100, ]
cand.mod <- candidateModels(list(phi = c(2, 4, 8), nu = c(1.5, 1.75),
                                 boundary = c(0.5)), "cartesian")

mod1 <- spGLMstack(y_pois ~ x1 + x2, data = dat, family = "poisson",
                   coords = as.matrix(dat[, c("s1", "s2")]), cor.fn = "matern",
                   candidate.models = cand.mod,
                   n.samples = 100,
                   loopd.controls = list(method = "CV", CV.K = 10, nMC = 500),
                   verbose = TRUE)

# Recover posterior samples of scale parameters
mod1.1 <- recoverGLMscale(mod1)

# sample from the stacked posterior distribution
post_samps <- stackedSampler(mod1.1)

---

Synthetic point-referenced spatial-temporal data

Description

Dataset of size 500 with spatial coordinates sampled uniformly from the unit square and temporal coordinates sampled uniformly from the unit interval, two covariates, an intercept and a slope of x1 that vary over space and time, and Gaussian and Poisson responses that share them.

Usage

data(simSpaceTime)

Format

a data.frame object with 500 rows and columns

⁠s1, s2⁠

2-D coordinates in the unit square.

t_coords

temporal coordinates in the unit interval.

⁠x1, x2⁠

covariates sampled from the standard normal distribution.

y_gauss

Gaussian response.

y_pois

Poisson count response.

z1_true

true spatial-temporal effect associated with the intercept.

z2_true

true spatial-temporal effect associated with x1.

Details

With $\ell = (s, t)$, the varying intercept is a wave travelling across space over time and the varying slope of $x_1$ changes with $s_2$ and $t$,

$z_1(\ell) = \sin(2 \pi s_1 - \pi t), \quad z_2(\ell) = \cos(2 \pi s_2) \cos(\pi t),$

deterministic surfaces that are not draws from Gaussian processes. The responses are generated as

$\begin{aligned} y_{\mathrm{gauss}}(\ell) &\sim N(2 + 5 x_1(\ell) - x_2(\ell) + z_1(\ell) + x_1(\ell) z_2(\ell), 0.5^2),\\ y_{\mathrm{pois}}(\ell) &\sim \mathrm{Poisson}(\exp\{2 - 0.5 x_1(\ell) + 0.3 x_2(\ell) + z_1(\ell) + x_1(\ell) z_2(\ell)\}). \end{aligned}$

This data can be generated with the code given in the example.

Author(s)

Soumyakanti Pan [email protected]

See Also

simSpatial

Examples

set.seed(1726)
n <- 500
s1 <- runif(n)
s2 <- runif(n)
t_coords <- runif(n)
x1 <- rnorm(n)
x2 <- rnorm(n)
z1 <- sin(2 * pi * s1 - pi * t_coords)
z2 <- cos(2 * pi * s2) * cos(pi * t_coords)
dat <- data.frame(
  s1 = s1, s2 = s2, t_coords = t_coords, x1 = x1, x2 = x2,
  y_gauss = rnorm(n, 2 + 5 * x1 - x2 + z1 + x1 * z2, sd = 0.5),
  y_pois = rpois(n, exp(2 - 0.5 * x1 + 0.3 * x2 + z1 + x1 * z2)),
  z1_true = z1, z2_true = z2
)
all.equal(dat, simSpaceTime)

---

Synthetic point-referenced spatial data

Description

Dataset of size 500 with spatial coordinates sampled uniformly from the unit square, two covariates, a spatial effect with a checkerboard pattern and one response for each of the Gaussian, Poisson, binomial and binary families, plus a Gaussian and a Poisson response with spatially varying coefficients (a varying intercept with a ring pattern and a varying slope of x1). All the responses share the covariates and the spatial effects, so that every spatial model in the package can be fitted to the same data.

Usage

data(simSpatial)

Format

a data.frame object with 500 rows and 13 columns

⁠s1, s2⁠

2-D coordinates in the unit square.

⁠x1, x2⁠

covariates sampled from the standard normal distribution.

y_gauss

Gaussian response.

y_pois

Poisson count response.

y_binom

binomial count response, out of n_trials trials.

n_trials

number of trials of the binomial response (5 to 20).

y_bin

binary response.

y_svc

Gaussian response with the spatially varying intercept z2 and the spatially varying slope 5 + z of x1.

y_svc_pois

Poisson count response with the spatially varying intercept z2 and the spatially varying slope -0.5 + z of x1.

z_true

true spatial effect (checkerboard) that generated the data; the spatially varying part of the slope of x1 in y_svc and y_svc_pois.

z2_true

true spatially varying intercept (rings) of y_svc and y_svc_pois.

Details

The spatial effect is the checkerboard

$z(s) = 1.5 \sin(2 \pi s_1) \sin(2 \pi s_2),$

and the spatially varying intercept is the rings

$z_2(s) = \cos(4 \pi \lVert s - s_0 \rVert) - c, \quad s_0 = (0.5, 0.5),$

where the constant $c \approx 0.129$ makes the mean of $z_2$ over the unit square zero (the mean of $z$ is zero by symmetry, and the two surfaces are orthogonal on the unit square). Both are deterministic surfaces, not draws from a Gaussian process. The responses are generated as

$\begin{aligned} y_{\mathrm{gauss}}(s) &\sim N(2 + 5 x_1(s) - x_2(s) + z(s), 0.5^2),\\ y_{\mathrm{pois}}(s) &\sim \mathrm{Poisson}(\exp\{2 - 0.5 x_1(s) + 0.3 x_2(s) + z(s)\}),\\ y_{\mathrm{binom}}(s) &\sim \mathrm{Binomial}(m(s), \mathrm{ilogit}\{0.5 - 0.5 x_1(s) + 0.5 x_2(s) + z(s)\}),\\ y_{\mathrm{bin}}(s) &\sim \mathrm{Bernoulli}(\mathrm{ilogit}\{0.5 - 0.5 x_1(s) + 0.5 x_2(s) + z(s)\}),\\ y_{\mathrm{svc}}(s) &\sim N(2 + z_2(s) + \{5 + z(s)\} x_1(s) - x_2(s), 0.5^2),\\ y_{\mathrm{svc\_pois}}(s) &\sim \mathrm{Poisson}(\exp\{2 + z_2(s) + \{-0.5 + z(s)\} x_1(s) + 0.3 x_2(s)\}), \end{aligned}$

where $m(s)$ is the number of trials. This data can be generated with the code given in the example.

Author(s)

Soumyakanti Pan [email protected]

See Also

simSpaceTime

Examples

set.seed(1729)
n <- 500
s1 <- runif(n)
s2 <- runif(n)
x1 <- rnorm(n)
x2 <- rnorm(n)
z <- 1.5 * sin(2 * pi * s1) * sin(2 * pi * s2)
ring <- function(s1, s2) cos(4 * pi * sqrt((s1 - 0.5)^2 + (s2 - 0.5)^2))
g <- (seq_len(1000) - 0.5) / 1000
z2 <- ring(s1, s2) - mean(outer(g, g, ring))
n_trials <- sample(5:20, n, replace = TRUE)
dat <- data.frame(
  s1 = s1, s2 = s2, x1 = x1, x2 = x2,
  y_gauss = rnorm(n, 2 + 5 * x1 - x2 + z, sd = 0.5),
  y_pois = rpois(n, exp(2 - 0.5 * x1 + 0.3 * x2 + z)),
  y_binom = rbinom(n, n_trials, plogis(0.5 - 0.5 * x1 + 0.5 * x2 + z)),
  n_trials = n_trials,
  y_bin = rbinom(n, 1, plogis(0.5 - 0.5 * x1 + 0.5 * x2 + z)),
  y_svc = rnorm(n, 2 + z2 + (5 + z) * x1 - x2, sd = 0.5),
  y_svc_pois = rpois(n, exp(2 + z2 + (-0.5 + z) * x1 + 0.3 * x2)),
  z_true = z,
  z2_true = z2
)
all.equal(dat, simSpatial)

---

Univariate Bayesian spatial generalized linear model

Description

Fits a Bayesian spatial generalized linear model with fixed values of spatial process parameters and some auxiliary model parameters. The output contains posterior samples of the fixed effects, spatial random effects and, if required, finds leave-one-out predictive densities.

Usage

spGLMexact(
  formula,
  data = parent.frame(),
  family,
  coords,
  cor.fn,
  priors,
  spParams,
  boundary = 0.5,
  n.samples,
  loopd = FALSE,
  loopd.method = "exact",
  CV.K = 10,
  loopd.nMC = 500,
  verbose = TRUE,
  ...
)

Arguments

formula	a symbolic description of the regression model to be fit. See example below.
data	an optional data frame containing the variables in the model. If not found in data, the variables are taken from environment(formula), typically the environment from which spGLMexact is called.
family	Specifies the distribution of the response as a member of the exponential family. Supported options are 'poisson', 'binomial' and 'binary'.
coords	an $n \times 2$ matrix of the observation coordinates in $\mathbb{R}^2$ (e.g., easting and northing).
cor.fn	a quoted keyword that specifies the correlation function used to model the spatial dependence structure among the observations. Supported covariance model key words are: 'exponential' and 'matern'. See below for details.
priors	(optional) a list with each tag corresponding to a hyperparameter name and containing hyperprior details. Valid tags include V.beta, nu.beta, nu.z and sigmaSq.xi. Values of nu.beta and nu.z must be at least 2.1. If not supplied, uses defaults.
spParams	fixed values of spatial process parameters.
boundary	Specifies the boundary adjustment parameter. Must be a real number between 0 and 1. Default is 0.5.
n.samples	number of posterior samples to be generated.
loopd	logical. If loopd=TRUE, returns leave-one-out predictive densities, using method as given by loopd.method. Default is FALSE.
loopd.method	character. Ignored if loopd=FALSE. If loopd=TRUE, valid inputs are 'exact', 'CV' and 'PSIS'. The option 'exact' corresponds to exact leave-one-out predictive densities which requires computation almost equivalent to fitting the model $n$ times. The options 'CV' and 'PSIS' are faster and they implement $K$-fold cross validation and Pareto-smoothed importance sampling to find approximate leave-one-out predictive densities (Vehtari et al. 2017).
CV.K	An integer between 10 and 20. Considered only if loopd.method='CV'. Default is 10 (as recommended in Vehtari et. al 2017).
loopd.nMC	Number of Monte Carlo samples to be used to evaluate leave-one-out predictive densities when loopd.method is set to either 'exact' or 'CV'.
verbose	logical. If verbose = TRUE, prints model description.
...	currently no additional argument.

Details

With this function, we fit a Bayesian hierarchical spatial generalized linear model by sampling exactly from the joint posterior distribution utilizing the generalized conjugate multivariate distribution theory (Bradley and Clinch 2024). Suppose $\chi = (s_1, \ldots, s_n)$ denotes the $n$ spatial locations the response $y$ is observed. Let $y(s)$ be the outcome at location $s$ endowed with a probability law from the natural exponential family, which we denote by

$y(s) \sim \mathrm{EF}(x(s)^\top \beta + z(s); b, \psi)$

for some positive parameter $b > 0$ and unit log partition function $\psi$. We consider the following response models based on the input supplied to the argument family.

'poisson'

It considers point-referenced Poisson responses $y(s) \sim \mathrm{Poisson}(e^{x(s)^\top \beta + z(s)})$. Here, $b = 1$ and $\psi(t) = e^t$.

'binomial'

It considers point-referenced binomial counts $y(s) \sim \mathrm{Binomial}(m(s), \pi(s))$ where, $m(s)$ denotes the total number of trials and probability of success $\pi(s) = \mathrm{ilogit}(x(s)^\top \beta + z(s))$ at location $s$. Here, $b = m(s)$ and $\psi(t) = \log(1+e^t)$.

'binary'

It considers point-referenced binary data (0 or, 1) i.e., $y(s) \sim \mathrm{Bernoulli}(\pi(s))$, where probability of success $\pi(s) = \mathrm{ilogit}(x(s)^\top \beta + z(s))$ at location $s$. Here, $b = 1$ and $\psi(t) = \log(1 + e^t)$.

The hierarchical model is given as

$\begin{aligned} y(s_i) &\mid \beta, z, \xi \sim EF(x(s_i)^\top \beta + z(s_i) + \xi_i - \mu_i; b_i, \psi_y), i = 1, \ldots, n\\ \xi &\mid \beta, z, \sigma^2_\xi, \alpha_\epsilon \sim \mathrm{GCM_c}(\cdots),\\ \beta &\mid \sigma^2_\beta \sim N(0, \sigma^2_\beta V_\beta), \quad \sigma^2_\beta \sim \mathrm{IG}(\nu_\beta/2, \nu_\beta/2)\\ z &\mid \sigma^2_z \sim N(0, \sigma^2_z R(\chi; \phi, \nu)), \quad \sigma^2_z \sim \mathrm{IG}(\nu_z/2, \nu_z/2), \end{aligned}$

where $\mu = (\mu_1, \ldots, \mu_n)^\top$ denotes the discrepancy parameter. We fix the spatial process parameters $\phi$ and $\nu$ and the hyperparameters $V_\beta$, $\nu_\beta$, $\nu_z$ and $\sigma^2_\xi$. The term $\xi$ is known as the fine-scale variation term which is given a conditional generalized conjugate multivariate distribution as prior. For more details, see Pan et al. 2025. Default values for $V_\beta$, $\nu_\beta$, $\nu_z$, $\sigma^2_\xi$ are diagonal with each diagonal element 100, 2.1, 2.1 and 0.1 respectively.

Value

An object of class spGLMexact, which is a list with the following tags -

priors

details of the priors used, containing the values of the boundary adjustment parameter (boundary), the variance parameter of the fine-scale variation term (simasq.xi) and others.

samples

a list of length 3, containing posterior samples of fixed effects (beta), spatial effects (z) and the fine-scale variation term (xi).

loopd

If loopd=TRUE, contains leave-one-out predictive densities.

model.params

Values of the fixed parameters that includes phi (spatial decay), nu (spatial smoothness; NA for the exponential correlation function).

diagnostics

a list of diagnostics. Element numerical is a data frame with one row and columns min.pivot (the smallest relative Cholesky pivot of the $n \times n$ correlation matrix; values below 1e-8 indicate a nearly singular correlation matrix), min.cor and max.cor (the correlations of the two farthest-apart and of the two closest locations; values of min.cor above 0.95 suggest an effective range far exceeding the extent of the data, values of max.cor below 0.05 nearly uncorrelated locations), obtained from quantities the fit computes anyway. If verbose = TRUE, a "Diagnostics" section is printed when any threshold is crossed.

The return object might include additional data that can be used for subsequent prediction and/or model fit evaluation.

Author(s)

Soumyakanti Pan [email protected]

References

Bradley JR, Clinch M (2024). "Generating Independent Replicates Directly from the Posterior Distribution for a Class of Spatial Hierarchical Models." Journal of Computational and Graphical Statistics, 0(0), 1-17. doi:10.1080/10618600.2024.2365728.

Pan S, Zhang L, Bradley JR, Banerjee S (2025). "Bayesian Inference for Spatial-temporal Non-Gaussian Data Using Predictive Stacking." Bayesian Analysis, In Press. doi:10.1214/25-BA1582.

Vehtari A, Gelman A, Gabry J (2017). "Practical Bayesian Model Evaluation Using Leave-One-out Cross-Validation and WAIC." Statistics and Computing, 27(5), 1413-1432. ISSN 0960-3174. doi:10.1007/s11222-016-9696-4.

See Also

spLMexact()

Examples

data(simSpatial)
dat <- simSpatial[1:100, ]
coords <- as.matrix(dat[, c("s1", "s2")])

# Example 1: spatial Poisson count data
mod1 <- spGLMexact(y_pois ~ x1 + x2, data = dat, family = "poisson",
                   coords = coords, cor.fn = "matern",
                   spParams = list(phi = 4, nu = 1.75),
                   n.samples = 100, verbose = TRUE)
post_beta <- mod1$samples$beta
print(t(apply(post_beta, 1, function(x) quantile(x, c(0.025, 0.5, 0.975)))))

# Example 2: spatial binomial count data
mod2 <- spGLMexact(cbind(y_binom, n_trials) ~ x1 + x2, data = dat,
                   family = "binomial", coords = coords, cor.fn = "matern",
                   spParams = list(phi = 4, nu = 1.75),
                   n.samples = 100, verbose = FALSE)

# Example 3: spatial binary data
mod3 <- spGLMexact(y_bin ~ x1 + x2, data = dat, family = "binary",
                   coords = coords, cor.fn = "matern",
                   spParams = list(phi = 4, nu = 1.75),
                   n.samples = 100, verbose = FALSE)

---

Bayesian spatial generalized linear model using predictive stacking

Description

Fits Bayesian spatial generalized linear model on a collection of candidate models constructed based on some candidate values of some model parameters specified by the user and subsequently combines inference by stacking predictive densities. See Pan, Zhang, Bradley, and Banerjee (2025) for more details.

Usage

spGLMstack(
  formula,
  data = parent.frame(),
  family,
  coords,
  cor.fn,
  priors,
  candidate.models,
  n.samples,
  loopd.controls,
  parallel = FALSE,
  solver = NULL,
  verbose = TRUE,
  ...
)

Arguments

formula	a symbolic description of the regression model to be fit. See example below.
data	an optional data frame containing the variables in the model. If not found in data, the variables are taken from environment(formula), typically the environment from which spLMstack is called.
family	Specifies the distribution of the response as a member of the exponential family. Supported options are 'poisson', 'binomial' and 'binary'.
coords	an $n \times 2$ matrix of the observation coordinates in $\mathbb{R}^2$ (e.g., easting and northing).
cor.fn	a quoted keyword that specifies the correlation function used to model the spatial dependence structure among the observations. Supported covariance model key words are: 'exponential' and 'matern'. See below for details.
priors	(optional) a list with each tag corresponding to a parameter name and containing prior details. Valid tags include V.beta, nu.beta, nu.z and sigmaSq.xi.
candidate.models	an object of class candidateModels containing a list of candidate models for stacking. See candidateModels() for details.
n.samples	number of posterior samples to be generated.
loopd.controls	a list with details on how leave-one-out predictive densities (LOO-PD) are to be calculated. Valid tags include method, CV.K, nMC and CV.update. The tag method can be either 'exact' or 'CV'. If sample size is more than 100, then the default is 'CV' with CV.K equal to its default value 10 (Gelman et al. 2024). The tag nMC decides how many Monte Carlo samples will be used to evaluate the leave-one-out predictive densities, which must be at least 500 (default). The tag CV.update is an advanced option, used only if method = 'CV', that decides how the pre-processing of the model fit on each fold is obtained, and should be changed with care as the faster choice depends on the BLAS library that R is linked with (see sessionInfo()). CV.update = 'update' obtains it by deletion updates of the full-data Cholesky factors, which is the faster choice with the reference BLAS that R ships with. CV.update = 'direct' recomputes it on each fold, which is faster only with an optimized BLAS such as OpenBLAS, Intel MKL or Apple Accelerate (vecLib), and is slower otherwise. The default CV.update = 'auto' uses 'direct' if such an optimized BLAS is detected from the library paths reported by R and 'update' otherwise; set it explicitly if the BLAS is not detected correctly (for example, an optimized BLAS installed in place of Rblas.dll on Windows). Both choices give the same results up to floating-point rounding, so only the run time is affected.
parallel	logical. If parallel=FALSE, the parallelization plan, if set up by the user, is ignored. If parallel=TRUE, the function inherits the parallelization plan that is set by the user via the function future::plan() only. Depending on the parallel backend available, users may choose their own plan. More details are available at https://cran.R-project.org/package=future.
solver	(optional) Specifies the name of the solver that will be used to obtain optimal stacking weights for each candidate model. Default order is c("CLARABEL", "ECOS", "SCS"). Users can use other solvers supported by the CVXR-package package.
verbose	logical. If TRUE, prints model-specific optimal stacking weights.
...	currently no additional argument.

Details

Instead of assigning a prior on the process parameters $\phi$ and $\nu$, the boundary adjustment parameter $\epsilon$, we consider a set of candidate models based on some candidate values of these parameters supplied by the user. Suppose the set of candidate models is $\mathcal{M} = \{M_1, \ldots, M_G\}$. Then for each $g = 1, \ldots, G$, we sample from the posterior distribution $p(\sigma^2, \beta, z \mid y, M_g)$ under the model $M_g$ and find leave-one-out predictive densities $p(y_i \mid y_{-i}, M_g)$. Then we solve the optimization problem

$\begin{aligned} \max_{w_1, \ldots, w_G}& \, \frac{1}{n} \sum_{i = 1}^n \log \sum_{g = 1}^G w_g p(y_i \mid y_{-i}, M_g) \\ \text{subject to} & \quad w_g \geq 0, \sum_{g = 1}^G w_g = 1 \end{aligned}$

to find the optimal stacking weights $\hat{w}_1, \ldots, \hat{w}_G$.

Value

An object of class spGLMstack, which is a list including the following tags -

family

the distribution of the responses as indicated in the function call

samples

a list of length equal to total number of candidate models with each entry corresponding to a list of length 3, containing posterior samples of fixed effects (beta), spatial effects (z) and fine-scale variation term (xi) for that particular model.

loopd

a list of length equal to total number of candidate models with each entry containing leave-one-out predictive densities under that particular model.

loopd.method

a list containing details of the algorithm used for calculation of leave-one-out predictive densities. For K-fold cross-validation, its tag cv.update records the pre-processing method ('update' or 'direct') that was used.

n.models

number of candidate models that are fit.

model.params

a list with one element per candidate model, each a named list of its parameters: phi, nu (NA for the exponential correlation function) and boundary (boundary adjustment parameter).

stacking.summary

a matrix with one row per candidate model, containing its parameters and its optimal stacking weight, for display.

stacking.weights

a numeric vector of length equal to the number of candidate models storing the optimal stacking weights.

run.time

a proc_time object with runtime details.

diagnostics

a list of diagnostics. Element numerical is a data frame with one row per candidate model and columns min.pivot (the smallest relative Cholesky pivot of the $n \times n$ correlation matrix; values below 1e-8 indicate a nearly singular correlation matrix), min.cor and max.cor (the correlations of the two farthest-apart and of the two closest locations; values of min.cor above 0.95 suggest an effective range far exceeding the extent of the data, values of max.cor below 0.05 nearly uncorrelated locations), obtained from quantities the fit computes anyway. Element solver describes the optimization for the stacking weights: the solver used (used) and its status (status), the installed and requested solvers, the search order, the attempts with their status, and whether the fallback loo::stacking_weights() was used. If verbose = TRUE, a "Diagnostics" section is printed if there is an issue: numerical flags of candidate models with stacking weight above 0.05 (extreme candidates with negligible weight are expected in a stacking grid and are only counted), and solver problems (a requested solver not installed, an inaccurate solution, or the fallback).

The return object might include additional data that is useful for subsequent prediction, model fit evaluation and other utilities.

Author(s)

Soumyakanti Pan [email protected],
Sudipto Banerjee [email protected]

References

Pan S, Zhang L, Bradley JR, Banerjee S (2025). "Bayesian Inference for Spatial-temporal Non-Gaussian Data Using Predictive Stacking." Bayesian Analysis, In Press. doi:10.1214/25-BA1582.

Vehtari A, Simpson D, Gelman A, Yao Y, Gabry J (2024). "Pareto Smoothed Importance Sampling." Journal of Machine Learning Research, 25(72), 1-58. URL https://jmlr.org/papers/v25/19-556.html.

See Also

spGLMexact(), spLMstack()

Examples

set.seed(1234)
data(simSpatial)
dat <- simSpatial[1:150, ]
cand.mod <- candidateModels(list(phi = c(2, 4, 8),
                                 nu = c(1.5, 1.75),
                                 boundary = c(0.75)), "cartesian")

mod1 <- spGLMstack(y_pois ~ x1 + x2, data = dat, family = "poisson",
                   coords = as.matrix(dat[, c("s1", "s2")]), cor.fn = "matern",
                   candidate.models = cand.mod,
                   n.samples = 500,
                   loopd.controls = list(method = "CV", CV.K = 10, nMC = 500),
                   parallel = FALSE, verbose = TRUE)

post_samps <- stackedSampler(mod1)
post_beta <- post_samps$beta
print(t(apply(post_beta, 1, function(x) quantile(x, c(0.025, 0.5, 0.975)))))

# compare the posterior medians of the spatial effects with the truth
cor(apply(post_samps$z, 1, median), dat$z_true)

---

Univariate Bayesian spatial linear model

Description

Fits a Bayesian spatial linear model with spatial process parameters and the noise-to-spatial variance ratio fixed to a value supplied by the user. The output contains posterior samples of the fixed effects, variance parameter, spatial random effects and, if required, leave-one-out predictive densities.

Usage

spLMexact(
  formula,
  data = parent.frame(),
  coords,
  cor.fn,
  priors = "flat",
  spParams,
  noise_sp_ratio,
  n.samples,
  loopd = FALSE,
  loopd.method = "exact",
  verbose = TRUE,
  ...
)

Arguments

formula	a symbolic description of the regression model to be fit. See example below.
data	an optional data frame containing the variables in the model. If not found in data, the variables are taken from environment(formula), typically the environment from which spLMexact is called.
coords	an $n \times 2$ matrix of the observation coordinates in $\mathbb{R}^2$ (e.g., easting and northing).
cor.fn	a quoted keyword that specifies the correlation function used to model the spatial dependence structure among the observations. Supported covariance model key words are: 'exponential' and 'matern'. See below for details.
priors	either "flat" (default), which assigns the prior $p(\beta, \sigma^2) \propto 1/\sigma^2$, or a list with tags beta.norm (a list containing $\mu_\beta$ and $V_\beta$) and/or sigma.sq.ig (a vector containing $a_\sigma$ and $b_\sigma$). A component not supplied in the list receives its flat prior, $p(\beta) \propto 1$ or $p(\sigma^2) \propto 1/\sigma^2$.
spParams	fixed value of spatial process parameters.
noise_sp_ratio	noise-to-spatial variance ratio.
n.samples	number of posterior samples to be generated.
loopd	logical. If loopd=TRUE, returns leave-one-out predictive densities, using method as given by loopd.method. Default is FALSE.
loopd.method	character. Ignored if loopd=FALSE. If loopd=TRUE, valid inputs are 'exact' and 'PSIS'. The option 'exact' corresponds to exact leave-one-out predictive densities which requires computation almost equivalent to fitting the model $n$ times. The option 'PSIS' is faster and finds approximate leave-one-out predictive densities using Pareto-smoothed importance sampling (Gelman et al. 2024).
verbose	logical. If verbose = TRUE, prints model description.
...	currently no additional argument.

Details

Suppose $\chi = (s_1, \ldots, s_n)$ denotes the $n$ spatial locations the response $y$ is observed. With this function, we fit a conjugate Bayesian hierarchical spatial model

$\begin{aligned} y \mid z, \beta, \sigma^2 &\sim N(X\beta + z, \sigma^2 I_n), \quad z \mid \sigma^2_z \sim N(0, \sigma^2_z R(\chi; \phi, \nu)), \\ \beta \mid \sigma^2 &\sim N(\mu_\beta, \sigma^2 V_\beta), \quad \sigma^2 \sim \mathrm{IG}(a_\sigma, b_\sigma) \end{aligned}$

where we fix the noise-to-spatial variance ratio $\delta^2 = \sigma^2 / \sigma^2_z$, the spatial process parameters $\phi$ and $\nu$, and the hyperparameters $\mu_\beta$, $V_\beta$, $a_\sigma$ and $b_\sigma$. If priors = "flat", we instead assign the prior $p(\beta, \sigma^2) \propto 1/\sigma^2$. We utilize a composition sampling strategy to sample the model parameters from their joint posterior distribution which can be written as

$p(\sigma^2, \beta, z \mid y) = p(\sigma^2 \mid y) \times p(\beta \mid \sigma^2, y) \times p(z \mid \beta, \sigma^2, y).$

We proceed by first sampling $\sigma^2$ from its marginal posterior, then given the samples of $\sigma^2$, we sample $\beta$ and subsequently, we sample $z$ conditioned on the posterior samples of $\beta$ and $\sigma^2$ (Banerjee 2020). Posterior samples of the spatial variance are obtained as $\sigma^2_z = \sigma^2 / \delta^2$.

Value

An object of class spLMexact, which is a list with the following tags -

samples

a list of length 4, containing posterior samples of fixed effects (beta), measurement error variance (sigmaSq), spatial variance (sigmaSq.z), and spatial effects (z).

loopd

If loopd=TRUE, contains leave-one-out predictive densities.

model.params

Values of the fixed parameters that includes phi (spatial decay), nu (spatial smoothness; NA for the exponential correlation function) and noise_sp_ratio (noise-to-spatial variance ratio).

diagnostics

a list of fit diagnostics, obtained from quantities the fit computes anyway. Element numerical is a data frame with one row and columns min.pivot (the smallest relative Cholesky pivot of the $n \times n$ factorizations; values below 1e-8 indicate a nearly singular covariance matrix), min.cor and max.cor (the correlations of the two farthest-apart and of the two closest locations; values of min.cor above 0.95 suggest an effective range far exceeding the extent of the data, values of max.cor below 0.05 nearly uncorrelated locations). If loopd.method = 'PSIS', element pareto is a list with the Pareto $k$ diagnostic values of the leave-one-out predictive densities (k), the threshold above which they are unreliable (threshold, Vehtari et al. 2024) and the number of values above it (n.high). If verbose = TRUE, a "Diagnostics" section is printed when any threshold is crossed.

The return object might include additional data used for subsequent prediction and/or model fit evaluation.

Author(s)

Soumyakanti Pan [email protected],
Sudipto Banerjee [email protected]

References

Banerjee S (2020). "Modeling massive spatial datasets using a conjugate Bayesian linear modeling framework." Spatial Statistics, 37, 100417. ISSN 2211-6753. doi:10.1016/j.spasta.2020.100417.

Vehtari A, Simpson D, Gelman A, Yao Y, Gabry J (2024). "Pareto Smoothed Importance Sampling." Journal of Machine Learning Research, 25(72), 1-58. URL https://jmlr.org/papers/v25/19-556.html.

See Also

spLMstack()

Examples

data(simSpatial)
dat <- simSpatial[1:100, ]

# setup prior list
muBeta <- c(0, 0, 0)
VBeta <- diag(100, 3)
sigmaSqIGa <- 2
sigmaSqIGb <- 0.1
prior_list <- list(beta.norm = list(muBeta, VBeta),
                   sigma.sq.ig = c(sigmaSqIGa, sigmaSqIGb))

mod1 <- spLMexact(y_gauss ~ x1 + x2, data = dat,
                  coords = as.matrix(dat[, c("s1", "s2")]),
                  cor.fn = "matern",
                  priors = prior_list,
                  spParams = list(phi = 4, nu = 1.75),
                  noise_sp_ratio = 0.1,
                  n.samples = 100,
                  loopd = TRUE, loopd.method = "exact")

post_beta <- mod1$samples$beta
print(t(apply(post_beta, 1, function(x) quantile(x, c(0.025, 0.5, 0.975)))))

# compare the posterior medians of the spatial effects with the truth
cor(apply(mod1$samples$z, 1, median), dat$z_true)

---

Bayesian spatial linear model using predictive stacking

Description

Fits Bayesian spatial linear model on a collection of candidate models constructed based on some candidate values of some model parameters specified by the user and subsequently combines inference by stacking predictive densities. See Zhang, Tang and Banerjee (2025) for more details.

Usage

spLMstack(
  formula,
  data = parent.frame(),
  coords,
  cor.fn,
  priors = "flat",
  candidate.models,
  n.samples,
  loopd.method,
  parallel = FALSE,
  solver = NULL,
  verbose = TRUE,
  ...
)

Arguments

formula	a symbolic description of the regression model to be fit. See example below.
data	an optional data frame containing the variables in the model. If not found in data, the variables are taken from environment(formula), typically the environment from which spLMstack is called.
coords	an $n \times 2$ matrix of the observation coordinates in $\mathbb{R}^2$ (e.g., easting and northing).
cor.fn	a quoted keyword that specifies the correlation function used to model the spatial dependence structure among the observations. Supported covariance model key words are: 'exponential' and 'matern'. See below for details.
priors	either "flat" (default), which assigns the prior $p(\beta, \sigma^2) \propto 1/\sigma^2$, or a list with tags beta.norm (a list containing $\mu_\beta$ and $V_\beta$) and/or sigma.sq.ig (a vector containing $a_\sigma$ and $b_\sigma$). A component not supplied in the list receives its flat prior, $p(\beta) \propto 1$ or $p(\sigma^2) \propto 1/\sigma^2$.
candidate.models	an object of class candidateModels containing a list of candidate models for stacking. See candidateModels() for details.
n.samples	number of posterior samples to be generated.
loopd.method	character. Valid inputs are 'exact' and 'PSIS'. The option 'exact' corresponds to exact leave-one-out predictive densities. The option 'PSIS' is faster, as it finds approximate leave-one-out predictive densities using Pareto-smoothed importance sampling (Gelman et al. 2024).
parallel	logical. If parallel=FALSE, the parallelization plan, if set up by the user, is ignored. If parallel=TRUE, the function inherits the parallelization plan that is set by the user via the function future::plan() only. Depending on the parallel backend available, users may choose their own plan. More details are available at https://cran.R-project.org/package=future.
solver	(optional) Specifies the name of the solver that will be used to obtain optimal stacking weights for each candidate model. Default order is c("CLARABEL", "ECOS", "SCS"). Users can use other solvers supported by the CVXR-package package.
verbose	logical. If TRUE, prints model-specific optimal stacking weights.
...	currently no additional argument.

Details

Instead of assigning a prior on the process parameters $\phi$ and $\nu$, noise-to-spatial variance ratio $\delta^2$, we consider a set of candidate models based on some candidate values of these parameters supplied by the user. Suppose the set of candidate models is $\mathcal{M} = \{M_1, \ldots, M_G\}$. Then for each $g = 1, \ldots, G$, we sample from the posterior distribution $p(\sigma^2, \beta, z \mid y, M_g)$ under the model $M_g$ and find leave-one-out predictive densities $p(y_i \mid y_{-i}, M_g)$. Then we solve the optimization problem

$\begin{aligned} \max_{w_1, \ldots, w_G}& \, \frac{1}{n} \sum_{i = 1}^n \log \sum_{g = 1}^G w_g p(y_i \mid y_{-i}, M_g) \\ \text{subject to} & \quad w_g \geq 0, \sum_{g = 1}^G w_g = 1 \end{aligned}$

to find the optimal stacking weights $\hat{w}_1, \ldots, \hat{w}_G$.

Value

An object of class spLMstack, which is a list including the following tags -

samples

a list of length equal to total number of candidate models with each entry corresponding to a list of length 4, containing posterior samples of fixed effects (beta), measurement error variance (sigmaSq), spatial variance (sigmaSq.z), and spatial effects (z) for that model.

loopd

a list of length equal to total number of candidate models with each entry containing leave-one-out predictive densities under that particular model.

n.models

number of candidate models that are fit.

model.params

a list with one element per candidate model, each a named list of its parameters: phi, nu (NA for the exponential correlation function) and noise_sp_ratio (noise-to-spatial variance ratio).

stacking.summary

a matrix with one row per candidate model, containing its parameters and its optimal stacking weight, for display.

stacking.weights

a numeric vector of length equal to the number of candidate models storing the optimal stacking weights.

run.time

a proc_time object with runtime details.

diagnostics

a list of diagnostics. Element numerical is a data frame with one row per candidate model and columns min.pivot (the smallest relative Cholesky pivot of the $n \times n$ factorizations; values below 1e-8 indicate a nearly singular covariance matrix), min.cor and max.cor (the correlations of the two farthest-apart and of the two closest locations; values of min.cor above 0.95 suggest an effective range far exceeding the extent of the data, values of max.cor below 0.05 nearly uncorrelated locations), obtained from quantities the fits compute anyway. If loopd.method = 'PSIS', element pareto is a list with the Pareto $k$ diagnostic values of each candidate model (k), the threshold above which they are unreliable (threshold) and the number of values above it for each model (n.high). Element solver describes the optimization for the stacking weights: the solver used (used) and its status (status), the installed and requested solvers, the search order, the attempts with their status, and whether the fallback loo::stacking_weights() was used. If verbose = TRUE, a "Diagnostics" section is printed if there is an issue: numerical flags of candidate models with stacking weight above 0.05 (extreme candidates with negligible weight are expected in a stacking grid and are only counted), Pareto $k$ values above the threshold for any model, and solver problems (a requested solver not installed, an inaccurate solution, or the fallback).

The return object might include additional data that is useful for subsequent prediction, model fit evaluation and other utilities.

Author(s)

Soumyakanti Pan [email protected],
Sudipto Banerjee [email protected]

References

Vehtari A, Simpson D, Gelman A, Yao Y, Gabry J (2024). "Pareto Smoothed Importance Sampling." Journal of Machine Learning Research, 25(72), 1-58. URL https://jmlr.org/papers/v25/19-556.html.

Zhang L, Tang W, Banerjee S (2025). "Bayesian Geostatistics Using Predictive Stacking." Journal of the American Statistical Association, In press. doi:10.1080/01621459.2025.2566449.

See Also

spLMexact(), spGLMstack()

Examples

set.seed(1234)
data(simSpatial)
dat <- simSpatial[1:100, ]

cand.mod <- candidateModels(list(phi = c(2, 4, 6),
                                 nu = c(1, 1.5, 1.75),
                                 noise_sp_ratio = c(0.1, 0.25)),
                            "cartesian")

mod1 <- spLMstack(y_gauss ~ x1 + x2, data = dat,
                  coords = as.matrix(dat[, c("s1", "s2")]),
                  cor.fn = "matern",
                  candidate.models = cand.mod,
                  n.samples = 1000, loopd.method = "exact",
                  parallel = FALSE, verbose = TRUE)

post_samps <- stackedSampler(mod1)
post_beta <- post_samps$beta
print(t(apply(post_beta, 1, function(x) quantile(x, c(0.025, 0.5, 0.975)))))

# compare the posterior medians of the spatial effects with the truth
cor(apply(post_samps$z, 1, median), dat$z_true)

---

Sample from the stacked posterior distribution

Description

A helper function to sample from the stacked posterior distribution to obtain final posterior samples that can be used for subsequent analysis. This function applies on outputs of functions spLMstack(), spGLMstack(), svcLMstack(), svcGLMstack(), stvcGLMstack() and stvcLMstack().

Usage

stackedSampler(mod_out, n.samples)

Arguments

mod_out	an object that is an output of a model fit or a prediction task, i.e., the class should be either spLMstack, pp.spLMstack, spGLMstack, pp.spGLMstack, stvcGLMstack, pp.stvcGLMstack, stvcLMstack, pp.stvcLMstack, svcLMstack, pp.svcLMstack, svcGLMstack or pp.svcGLMstack.
n.samples	(optional) If missing, inherits the number of posterior samples from the original output. Otherwise, it specifies number of posterior samples to draw from the stacked posterior. If it exceeds the number of posterior draws used in the original function, then a message is thrown and the samples are obtained by resampling. We recommended running the original model fit/prediction with enough samples.

Details

After obtaining the optimal stacking weights $\hat{w}_1, \ldots, \hat{w}_G$, posterior inference of quantities of interest subsequently proceed from the stacked posterior,

$\tilde{p}(\cdot \mid y) = \sum_{g = 1}^G \hat{w}_g p(\cdot \mid y, M_g),$

where $\mathcal{M} = \{M_1, \ldots, M_g\}$ is the collection of candidate models.

Value

An object of class stacked_posterior, which is a list that includes the following tags -

beta

samples of the fixed effect from the stacked joint posterior.

z

samples of the spatial random effects from the stacked joint posterior.

The list may also include other scale parameters corresponding to the model.

Author(s)

Soumyakanti Pan [email protected],
Sudipto Banerjee [email protected]

See Also

spLMstack(), spGLMstack(), svcLMstack(), svcGLMstack(), stvcGLMstack(), stvcLMstack()

Examples

set.seed(1234)
data(simSpatial)
dat <- simSpatial[1:100, ]
cand.mod <- candidateModels(list(phi = c(2, 4, 6),
                                 nu = c(1, 1.5, 1.75),
                                 noise_sp_ratio = c(0.1)), "cartesian")

mod1 <- spLMstack(y_gauss ~ x1 + x2, data = dat,
                  coords = as.matrix(dat[, c("s1", "s2")]),
                  cor.fn = "matern",
                  candidate.models = cand.mod,
                  n.samples = 1000, loopd.method = "exact",
                  parallel = FALSE, verbose = TRUE)
print(mod1$diagnostics$solver$used)
print(mod1$diagnostics$solver$status)
print(mod1$run.time)

post_samps <- stackedSampler(mod1)
post_beta <- post_samps$beta
print(t(apply(post_beta, 1, function(x) quantile(x, c(0.025, 0.5, 0.975)))))

---

Bayesian spatially-temporally varying generalized linear model

Description

Fits a Bayesian generalized linear model with spatially-temporally varying coefficients under fixed values of spatial-temporal process parameters and some auxiliary model parameters. The output contains posterior samples of the fixed effects, spatial-temporal random effects and, if required, finds leave-one-out predictive densities.

Usage

stvcGLMexact(
  formula,
  data = parent.frame(),
  family,
  sp_coords,
  time_coords,
  cor.fn,
  process.type,
  sptParams,
  priors,
  boundary = 0.5,
  n.samples,
  loopd = FALSE,
  loopd.method = "exact",
  CV.K = 10,
  loopd.nMC = 500,
  verbose = TRUE,
  ...
)

Arguments

formula	a symbolic description of the regression model to be fit. Variables in parenthesis are assigned spatially-temporally varying coefficients. See examples.
data	an optional data frame containing the variables in the model. If not found in data, the variables are taken from environment(formula), typically the environment from which stvcGLMexact is called.
family	Specifies the distribution of the response as a member of the exponential family. Supported options are 'poisson', 'binomial' and 'binary'.
sp_coords	an $n \times 2$ matrix of the observation spatial coordinates in $\mathbb{R}^2$ (e.g., easting and northing).
time_coords	an $n \times 1$ matrix of the observation temporal coordinates in $\mathcal{T} \subseteq [0, \infty)$.
cor.fn	a quoted keyword that specifies the correlation function used to model the spatial-temporal dependence structure among the observations. Supported covariance model key words are: 'gneiting-decay' (Gneiting and Guttorp 2010). See below for details.
process.type	a quoted keyword specifying the model for the spatial-temporal process. Supported keywords are 'independent' which indicates independent processes for each varying coefficients characterized by different process parameters, independent.shared implies independent processes for the varying coefficients that shares common process parameters, and multivariate implies correlated processes for the varying coefficients modeled by a multivariate Gaussian process with an inverse-Wishart prior on the correlation matrix. The input for sptParams and priors must be given accordingly.
sptParams	fixed values of spatial-temporal process parameters in usually a list of length 2. If cor.fn='gneiting-decay', then it is a list of length 2 with tags phi_s and phi_t. If process.type='independent', then phi_s and phi_t contain fixed values of the $r$ spatial-temporal processes, otherwise they will contain scalars. See examples below.
priors	(optional) a list with each tag corresponding to a hyperparameter name and containing hyperprior details. Valid tags include V.beta, nu.beta, nu.z, sigmaSq.xi and IW.scale. Values of nu.beta and nu.z must be at least 2.1. If not supplied, uses defaults.
boundary	Specifies the boundary adjustment parameter. Must be a real number between 0 and 1. Default is 0.5.
n.samples	number of posterior samples to be generated.
loopd	logical. If loopd=TRUE, returns leave-one-out predictive densities, using method as given by loopd.method. Default is FALSE.
loopd.method	character. Ignored if loopd=FALSE. If loopd=TRUE, valid inputs are 'exact', 'CV'. The option 'exact' corresponds to exact leave-one-out predictive densities which requires computation almost equivalent to fitting the model $n$ times. The options 'CV' is faster as it implements $K$-fold cross validation to find approximate leave-one-out predictive densities (Vehtari et al. 2017).
CV.K	An integer between 10 and 20. Considered only if loopd.method='CV'. Default is 10 (as recommended in Vehtari et. al 2017).
loopd.nMC	Number of Monte Carlo samples to be used to evaluate leave-one-out predictive densities when loopd.method is set to either 'exact' or 'CV'.
verbose	logical. If verbose = TRUE, prints model description.
...	currently no additional argument.

Details

With this function, we fit a Bayesian hierarchical spatially-temporally varying generalized linear model by sampling exactly from the joint posterior distribution utilizing the generalized conjugate multivariate distribution theory (Bradley and Clinch 2024). Suppose $\chi = (\ell_1, \ldots, \ell_n)$ denotes the $n$ spatial-temporal co-ordinates in $\mathcal{L} = \mathcal{S} \times \mathcal{T}$, the response $y$ is observed. Let $y(\ell)$ be the outcome at the co-ordinate $\ell$ endowed with a probability law from the natural exponential family, which we denote by

$y(\ell) \sim \mathrm{EF}(x(\ell)^\top \beta + \tilde{x}(\ell)^\top z(\ell); b(\ell), \psi)$

for some positive parameter $b(\ell) > 0$ and unit log partition function $\psi$. Here, $\tilde{x}(\ell)$ denotes covariates with spatially-temporally varying coefficients We consider the following response models based on the input supplied to the argument family.

'poisson'

It considers point-referenced Poisson responses $y(\ell) \sim \mathrm{Poisson}(e^{x(\ell)^\top \beta + \tilde{x}(\ell)^\top z(\ell)})$. Here, $b(\ell) = 1$ and $\psi(t) = e^t$.

'binomial'

It considers point-referenced binomial counts $y(\ell) \sim \mathrm{Binomial}(m(\ell), \pi(\ell))$ where, $m(\ell)$ denotes the total number of trials and probability of success $\pi(\ell) = \mathrm{ilogit}(x(\ell)^\top \beta + \tilde{x}(\ell)^\top z(\ell))$ at spatial-temporal co-ordinate $\ell$. Here, $b = m(\ell)$ and $\psi(t) = \log(1+e^t)$.

'binary'

It considers point-referenced binary data (0 or, 1) i.e., $y(\ell) \sim \mathrm{Bernoulli}(\pi(\ell))$, where probability of success $\pi(\ell) = \mathrm{ilogit}(x(\ell)^\top \beta + \tilde{x}(\ell)^\top z(\ell))$ at spatial-temporal co-ordinate $\ell$. Here, $b(\ell) = 1$ and $\psi(t) = \log(1 + e^t)$.

The hierarchical model is given as

$\begin{aligned} y(\ell_i) &\mid \beta, z, \xi \sim EF(x(\ell_i)^\top \beta + \tilde{x}(\ell_i)^\top z(s_i) + \xi_i - \mu_i; b_i, \psi_y), i = 1, \ldots, n\\ \xi &\mid \beta, z, \sigma^2_\xi, \alpha_\epsilon \sim \mathrm{GCM_c}(\cdots),\\ \beta &\mid \sigma^2_\beta \sim N(0, \sigma^2_\beta V_\beta), \quad \sigma^2_\beta \sim \mathrm{IG}(\nu_\beta/2, \nu_\beta/2)\\ z_j &\mid \sigma^2_{z_j} \sim N(0, \sigma^2_{z_j} R(\chi; \phi_s, \phi_t)), \quad \sigma^2_{z_j} \sim \mathrm{IG}(\nu_z/2, \nu_z/2), j = 1, \ldots, r \end{aligned}$

where $\mu = (\mu_1, \ldots, \mu_n)^\top$ denotes the discrepancy parameter. We fix the spatial-temporal process parameters $\phi_s$ and $\phi_t$ and the hyperparameters $V_\beta$, $\nu_\beta$, $\nu_z$ and $\sigma^2_\xi$. The term $\xi$ is known as the fine-scale variation term which is given a conditional generalized conjugate multivariate distribution as prior. For more details, see Pan et al. 2025. Default values for $V_\beta$, $\nu_\beta$, $\nu_z$, $\sigma^2_\xi$ are diagonal with each diagonal element 100, 2.1, 2.1 and 0.1 respectively. If process.type = 'multivariate', the $r \times r$ cross-covariance matrix $\Sigma$ of the processes has the prior $\Sigma \sim \mathrm{IW}(\nu_z + 2r, \Psi)$, with $\Psi$ given by IW.scale; its default $\Psi = (\nu_z + r - 1) I_r$ gives $E[\Sigma] = I_r$. The widths of the posterior credible intervals of the processes grow with the prior scale.

Value

An object of class stvcGLMexact, which is a list with the following tags -

priors

details of the priors used, containing the values of the boundary adjustment parameter (boundary), the variance parameter of the fine-scale variation term (simasq.xi) and others.

samples

a list of length 3, containing posterior samples of fixed effects (beta), spatial-temporal effects (z) and the fine-scale variation term (xi). The element with tag z will again be a list of length $r$, each containing posterior samples of the spatial-temporal random effects corresponding to each varying coefficient.

loopd

If loopd=TRUE, contains leave-one-out predictive densities.

model.params

Values of the fixed parameters that includes phi_s (spatial decay), phi_t (temporal smoothness).

diagnostics

a list of diagnostics. Element numerical is a data frame with one row (one row per process, with columns model and process, if process.type = 'independent') and columns min.pivot (the smallest relative Cholesky pivot of the $n \times n$ spatial-temporal correlation matrix; values below 1e-8 indicate a nearly singular correlation matrix), min.cor and max.cor (the correlations of the two farthest-apart and of the two closest space-time locations; values of min.cor above 0.95 suggest an effective range far exceeding the extent of the data, values of max.cor below 0.05 nearly uncorrelated space-time locations), obtained from quantities the fit computes anyway. If verbose = TRUE, a "Diagnostics" section is printed when any threshold is crossed.

The return object might include additional data that can be used for subsequent prediction and/or model fit evaluation.

Author(s)

Soumyakanti Pan [email protected]

References

Bradley JR, Clinch M (2024). "Generating Independent Replicates Directly from the Posterior Distribution for a Class of Spatial Hierarchical Models." Journal of Computational and Graphical Statistics, 0(0), 1-17. doi:10.1080/10618600.2024.2365728.

T. Gneiting and P. Guttorp (2010). "Continuous-parameter spatio-temporal processes." In A.E. Gelfand, P.J. Diggle, M. Fuentes, and P Guttorp, editors, Handbook of Spatial Statistics, Chapman & Hall CRC Handbooks of Modern Statistical Methods, p 427–436. Taylor and Francis.

Pan S, Zhang L, Bradley JR, Banerjee S (2025). "Bayesian Inference for Spatial-temporal Non-Gaussian Data Using Predictive Stacking." Bayesian Analysis, In Press. doi:10.1214/25-BA1582.

Vehtari A, Gelman A, Gabry J (2017). "Practical Bayesian Model Evaluation Using Leave-One-out Cross-Validation and WAIC." Statistics and Computing, 27(5), 1413-1432. ISSN 0960-3174. doi:10.1007/s11222-016-9696-4.

See Also

spGLMexact()

Examples

data(simSpaceTime)
dat <- simSpaceTime[1:100, ]

# Fit a spatial-temporal varying coefficient Poisson GLM
mod1 <- stvcGLMexact(y_pois ~ x1 + x2 + (x1), data = dat, family = "poisson",
                     sp_coords = as.matrix(dat[, c("s1", "s2")]),
                     time_coords = as.matrix(dat[, "t_coords"]),
                     cor.fn = "gneiting-decay",
                     process.type = "independent",
                     sptParams = list(phi_s = c(3, 3), phi_t = c(4, 4)),
                     verbose = FALSE, n.samples = 100)

# rows 1, ..., n of z hold the varying intercept, rows n+1, ..., 2n the
# varying slope of x1
n <- nrow(dat)
z_hat <- apply(mod1$samples$z, 1, median)
cor(z_hat[1:n], dat$z1_true)

---

Bayesian spatially-temporally varying coefficients generalized linear model using predictive stacking

Description

Fits Bayesian spatial-temporal generalized linear model with spatially-temporally varying coefficients on a collection of candidate models constructed based on some candidate values of some model parameters specified by the user and subsequently combines inference by stacking predictive densities. See Pan, Zhang, Bradley, and Banerjee (2025) for more details.

Usage

stvcGLMstack(
  formula,
  data = parent.frame(),
  family,
  sp_coords,
  time_coords,
  cor.fn,
  process.type,
  priors,
  candidate.models,
  n.samples,
  loopd.controls,
  parallel = FALSE,
  solver = NULL,
  verbose = TRUE,
  ...
)

Arguments

formula	a symbolic description of the regression model to be fit. Variables in parenthesis are assigned spatially-temporally varying coefficients. See examples.
data	an optional data frame containing the variables in the model. If not found in data, the variables are taken from environment(formula), typically the environment from which stvcGLMstack is called.
family	Specifies the distribution of the response as a member of the exponential family. Supported options are 'poisson', 'binomial' and 'binary'.
sp_coords	an $n \times 2$ matrix of the observation spatial coordinates in $\mathbb{R}^2$ (e.g., easting and northing).
time_coords	an $n \times 1$ matrix of the observation temporal coordinates in $\mathcal{T} \subseteq [0, \infty)$.
cor.fn	a quoted keyword that specifies the correlation function used to model the spatial-temporal dependence structure among the observations. Supported covariance model key words are: 'gneiting-decay' (Gneiting and Guttorp 2010). See below for details.
process.type	a quoted keyword specifying the model for the spatial-temporal process. Supported keywords are 'independent' which indicates independent processes for each varying coefficients characterized by different process parameters, independent.shared implies independent processes for the varying coefficients that shares common process parameters, and multivariate implies correlated processes for the varying coefficients modeled by a multivariate Gaussian process with an inverse-Wishart prior on the correlation matrix. The input for sptParams and priors must be given accordingly.
priors	(optional) a list with each tag corresponding to a hyperparameter name and containing hyperprior details. Valid tags include V.beta, nu.beta, nu.z, sigmaSq.xi and IW.scale. Values of nu.beta and nu.z must be at least 2.1. If not supplied, uses defaults.
candidate.models	an object of class candidateModels containing a list of candidate models for stacking. See candidateModels() for details.
n.samples	number of samples to be drawn from the posterior distribution.
loopd.controls	a list with details on how leave-one-out predictive densities (LOO-PD) are to be calculated. Valid tags include method, CV.K, nMC and CV.update. The tag method can be either 'exact' or 'CV'. If sample size is more than 100, then the default is 'CV' with CV.K equal to its default value 10 (Gelman et al. 2024). The tag nMC decides how many Monte Carlo samples will be used to evaluate the leave-one-out predictive densities, which must be at least 500 (default). The tag CV.update is an advanced option, used only if method = 'CV', that decides how the pre-processing of the model fit on each fold is obtained, and should be changed with care as the faster choice depends on the BLAS library that R is linked with (see sessionInfo()). CV.update = 'update' obtains it by deletion updates of the full-data Cholesky factors, which is the faster choice with the reference BLAS that R ships with. CV.update = 'direct' recomputes it on each fold, which is faster only with an optimized BLAS such as OpenBLAS, Intel MKL or Apple Accelerate (vecLib), and is slower otherwise. The default CV.update = 'auto' uses 'direct' if such an optimized BLAS is detected from the library paths reported by R and 'update' otherwise; set it explicitly if the BLAS is not detected correctly (for example, an optimized BLAS installed in place of Rblas.dll on Windows). Both choices give the same results up to floating-point rounding, so only the run time is affected.
parallel	logical. If parallel=FALSE, the parallelization plan, if set up by the user, is ignored. If parallel=TRUE, the function inherits the parallelization plan that is set by the user via the function future::plan() only. Depending on the parallel backend available, users may choose their own plan. More details are available at https://cran.R-project.org/package=future.
solver	(optional) Specifies the name of the solver that will be used to obtain optimal stacking weights for each candidate model. Default order is c("CLARABEL", "ECOS", "SCS"). Users can use other solvers supported by the CVXR-package package.
verbose	logical. If TRUE, prints model-specific optimal stacking weights.
...	currently no additional argument.

Value

An object of class stvcGLMstack, which is a list including the following tags -

samples

a list of length equal to total number of candidate models with each entry corresponding to a list of length 3, containing posterior samples of fixed effects (beta), spatial effects (z), and fine scale variation xi for that model.

loopd

a list of length equal to total number of candidate models with each entry containing leave-one-out predictive densities under that particular model.

n.models

number of candidate models that are fit.

model.params

a list with one element per candidate model, each a named list of its parameters: phi_s, phi_t (vectors of length $r$ if process.type = 'independent') and boundary.

stacking.summary

a matrix with one row per candidate model, containing its parameters and its optimal stacking weight, for display.

stacking.weights

a numeric vector of length equal to the number of candidate models storing the optimal stacking weights.

run.time

a proc_time object with runtime details.

diagnostics

a list of diagnostics. Element numerical is a data frame with one row per candidate model (per candidate model and process, with columns model and process, if process.type = 'independent') and columns min.pivot (the smallest relative Cholesky pivot of the $n \times n$ spatial-temporal correlation matrix; values below 1e-8 indicate a nearly singular correlation matrix), min.cor and max.cor (the correlations of the two farthest-apart and of the two closest space-time locations; values of min.cor above 0.95 suggest an effective range far exceeding the extent of the data, values of max.cor below 0.05 nearly uncorrelated space-time locations), obtained from quantities the fit computes anyway. Element solver describes the optimization for the stacking weights: the solver used (used) and its status (status), the installed and requested solvers, the search order, the attempts with their status, and whether the fallback loo::stacking_weights() was used. If verbose = TRUE, a "Diagnostics" section is printed if there is an issue: numerical flags of candidate models with stacking weight above 0.05 (extreme candidates with negligible weight are expected in a stacking grid and are only counted), and solver problems (a requested solver not installed, an inaccurate solution, or the fallback).

This object can be further used to recover posterior samples of the scale parameters in the model, and subsequrently, to make predictions at new locations or times using the function posteriorPredict().

Examples

set.seed(1234)
data(simSpaceTime)
dat <- simSpaceTime[1:100, ]

# create list of candidate models (multivariate)
mod.list2 <- candidateModels(list(phi_s = list(2, 3, 4),
                                  phi_t = list(4, 8),
                                  boundary = c(0.75)), "cartesian")

# fit a spatial-temporal varying coefficient model using predictive stacking
mod1 <- stvcGLMstack(y_pois ~ x1 + x2 + (x1), data = dat, family = "poisson",
                     sp_coords = as.matrix(dat[, c("s1", "s2")]),
                     time_coords = as.matrix(dat[, "t_coords"]),
                     cor.fn = "gneiting-decay",
                     process.type = "multivariate",
                     candidate.models = mod.list2,
                     loopd.controls = list(method = "CV", CV.K = 10, nMC = 500),
                     n.samples = 500)
post_samps <- stackedSampler(mod1)

---