Bayesian spatially-temporally varying coefficients linear model using predictive stacking

Fits Bayesian linear models with spatially-temporally varying coefficients for a Gaussian response on a collection of candidate models, constructed from candidate values of the spatial-temporal process parameters and the noise-to-spatial variance ratios supplied by the user, and combines inference by stacking predictive densities.

Usage

stvcLMstack(
  formula,
  data = parent.frame(),
  sp_coords,
  time_coords,
  cor.fn,
  process.type,
  priors = "flat",
  candidate.models,
  n.samples,
  loopd.method = "exact",
  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 stvcLMstack is called.
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).
process.type: a quoted keyword specifying the model for the spatial-temporal processes of the varying coefficients: 'independent' or 'independent.shared'. See stvcLMexact().
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.
candidate.models: an object of class candidateModels containing a list of candidate models for stacking, each with tags phi_s, phi_t and noise_sp_ratio: vectors of length $r$ if process.type = 'independent' (use list() entries in candidateModels()), otherwise scalars. See candidateModels() for details.
n.samples: number of posterior samples to be generated.
loopd.method: character. Valid inputs are 'exact' (default) and 'PSIS'. The option 'exact' finds the exact leave-one-out predictive densities in closed form. The option 'PSIS' uses Pareto-smoothed importance sampling (Vehtari et al. 2024); with many latent effects its Pareto $k$ diagnostics are often high, so 'exact' is recommended.
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 stvcLMstack, 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), the noise variance (sigmaSq), the process variances (sigmaSq.z) and the spatial-temporal effects (z) for that model, as in stvcLMexact().
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 and noise_sp_ratio (vectors of length $r$ if process.type = 'independent').
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, min.cor and max.cor, as in stvcLMexact(). If loopd.method = 'PSIS', element pareto contains 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.

This object can be used to make predictions at new locations or times with posteriorPredict() and to sample from the stacked posterior with stackedSampler().

Details

Instead of assigning priors on the process parameters $\phi_s$, $\phi_t$ and the noise-to-spatial variance ratios $\delta^2$, we consider a set of candidate models $\mathcal{M} = \{M_1, \ldots, M_G\}$ based on candidate values of these parameters. For each $g$, we sample exactly from the posterior distribution $p(\sigma^2, \beta, z \mid y, M_g)$ (see stvcLMexact()) and find the leave-one-out predictive densities $p(y_i \mid y_{-i}, M_g)$. The stacking weights solve $$ \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} $$ Candidate models that share the process parameters $(\phi_s, \phi_t)$ are fitted together, building and factorizing the spatial-temporal correlation matrices only once.

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 .

Author

Soumyakanti Pan span18@ucla.edu,
Sudipto Banerjee sudipto@ucla.edu

Examples

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

# processes sharing (phi_s, phi_t, noise_sp_ratio): scalar candidates
mod.list <- candidateModels(list(phi_s = c(1, 2, 3), phi_t = c(2, 4, 8),
                                 noise_sp_ratio = c(0.1, 0.25)), "cartesian")
mod1 <- stvcLMstack(y_gauss ~ x1 + x2 + (x1), data = dat,
                    sp_coords = as.matrix(dat[, c("s1", "s2")]),
                    time_coords = as.matrix(dat[, "t_coords"]),
                    cor.fn = "gneiting-decay",
                    process.type = "independent.shared",
                    candidate.models = mod.list,
                    n.samples = 500)
#> 
#> STACKING WEIGHTS:
#> 
#>            | phi_s | phi_t | noise_sp_ratio | weight |
#> +----------+-------+-------+----------------+--------+
#> | Model 1  |      1|      2|            0.10| 0.000  |
#> | Model 2  |      2|      2|            0.10| 0.000  |
#> | Model 3  |      3|      2|            0.10| 0.000  |
#> | Model 4  |      1|      4|            0.10| 0.000  |
#> | Model 5  |      2|      4|            0.10| 0.541  |
#> | Model 6  |      3|      4|            0.10| 0.000  |
#> | Model 7  |      1|      8|            0.10| 0.000  |
#> | Model 8  |      2|      8|            0.10| 0.000  |
#> | Model 9  |      3|      8|            0.10| 0.000  |
#> | Model 10 |      1|      2|            0.25| 0.000  |
#> | Model 11 |      2|      2|            0.25| 0.000  |
#> | Model 12 |      3|      2|            0.25| 0.000  |
#> | Model 13 |      1|      4|            0.25| 0.000  |
#> | Model 14 |      2|      4|            0.25| 0.204  |
#> | Model 15 |      3|      4|            0.25| 0.000  |
#> | Model 16 |      1|      8|            0.25| 0.000  |
#> | Model 17 |      2|      8|            0.25| 0.254  |
#> | Model 18 |      3|      8|            0.25| 0.000  |
#> +----------+-------+-------+----------------+--------+
#> 
post_samps <- stackedSampler(mod1)
n <- nrow(dat)
cor(apply(post_samps$z[1:n, ], 1, median), dat$z1_true)
#> [1] 0.9084215

# independent processes (r = 2): vector-valued candidates via list()
mod.list2 <- candidateModels(list(phi_s = list(c(2, 2), c(3, 3)),
                                  phi_t = list(c(4, 4), c(8, 8)),
                                  noise_sp_ratio = list(c(0.1, 0.25), c(0.25, 0.5))),
                             "cartesian")
mod2 <- stvcLMstack(y_gauss ~ x1 + x2 + (x1), data = dat,
                    sp_coords = as.matrix(dat[, c("s1", "s2")]),
                    time_coords = as.matrix(dat[, "t_coords"]),
                    cor.fn = "gneiting-decay",
                    process.type = "independent",
                    candidate.models = mod.list2,
                    n.samples = 500, verbose = FALSE)