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

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)
#> 
#> STACKING WEIGHTS:
#> 
#>           | phi_s | phi_t | boundary | weight |
#> +---------+-------+-------+----------+--------+
#> | Model 1 |      2|      4|      0.75| 0.713  |
#> | Model 2 |      3|      4|      0.75| 0.000  |
#> | Model 3 |      4|      4|      0.75| 0.000  |
#> | Model 4 |      2|      8|      0.75| 0.287  |
#> | Model 5 |      3|      8|      0.75| 0.000  |
#> | Model 6 |      4|      8|      0.75| 0.000  |
#> +---------+-------+-------+----------+--------+
#> 
post_samps <- stackedSampler(mod1)