Spatial-Temporal Regression Models • spStack

In this article, we discuss the following functions -

These functions can be used to fit Gaussian and non-Gaussian spatial-temporal point-referenced data with spatially-temporally varying coefficients.

library(patchwork)
set.seed(1729)

Data

The package lazy-loads the synthetic dataset simSpaceTime, observed at 500 space-time coordinates $\ell = (s, t)$ with $s$ in the unit square and $t$ in the unit interval. It has two covariates x1 and x2, a varying intercept z1_true that is a wave travelling across space over time, $z_1(\ell) = \sin(2\pi s_1 - \pi t)$ , and a varying slope of x1, z2_true, $z_2(\ell) = \cos(2\pi s_2) \cos(\pi t)$ . The Gaussian response y_gauss and the Poisson response y_pois share them. Both surfaces are deterministic, not draws from Gaussian processes. See ?simSpaceTime for the code that generated the data. We use the first 200 space-time coordinates for the following analysis.

library(spStack)
data("simSpaceTime")
n_train <- 200
dat <- simSpaceTime[1:n_train, ]
head(dat)

##          s1        s2  t_coords          x1         x2   y_gauss y_pois
## 1 0.2216173 0.9919157 0.8458822 -0.16471416  2.0451334 -1.032669      4
## 2 0.5937945 0.3285576 0.7965761 -1.47214370 -1.1420950 -4.159786     24
## 3 0.5511327 0.1808101 0.9182483  0.53238120 -0.1035955  5.339375      9
## 4 0.2295121 0.7660980 0.8099342  1.64998720  1.3164616  7.439179      0
## 5 0.9524661 0.9817298 0.7379233 -0.03672285 -0.9846361  2.723525      3
## 6 0.6469231 0.4993637 0.4131629 -0.34204648 -1.0552852  1.432997     11
##      z1_true     z2_true
## 1 -0.9535938 -0.88391758
## 2  0.9419539  0.38028823
## 3  0.5464392 -0.40735390
## 4 -0.8923008 -0.08350197
## 5 -0.5009308 -0.67530256
## 6  0.3661277 -0.26943332

Formula for varying coefficients model

We define the spatially-temporally varying coefficients model using a formula, similar to that in the widely used lm() function in the stats package. Suppose $\ell = (s, t)$ refers to a space-time coordinate. See “Technical Overview” for more details. The formula y_gauss ~ x1 + x2 + (x1) corresponds to the spatial-temporal linear model $y(\ell) = \beta_0 + \beta_1 x_1(\ell) + \beta_2 x_2(\ell) + z_1(\ell) + x_1(\ell) z_2(\ell) + \epsilon(\ell)\;,$ where y_gauss corresponds to the response variable $y(\ell)$ , which is regressed on the predictors x1 and x2 given by $x_1(\ell)$ and $x_2(\ell)$ . The model variables specified outside the parentheses correspond to predictors with fixed effects, and the model inside the parentheses correspond to variables with spatial-temporal varying coefficients. The intercept is automatically considered within both the fixed and varying coefficient components of the model, and hence y_gauss ~ x1 + x2 + (x1) is functionally equivalent to y_gauss ~ 1 + x1 + x2 + (1 + x1). For now, we only support the cor.fn="gneiting-decay" covariogram. To implement a model with just a spatial-temporal random effect, one may specify the formula y_gauss ~ x1 + x2 + (1).

Bayesian Gaussian spatially-temporally varying coefficient models

For a Gaussian response, the model is $y(\ell) = x(\ell)^{ \scriptstyle \top }\beta + \tilde{x}(\ell)^{ \scriptstyle \top }z(\ell) + \epsilon(\ell), \quad \epsilon(\ell) \sim \mathrm{N}(0, \sigma^2),$ where each of the $r$ varying coefficients $z_j$ is an independent Gaussian process with variance $\sigma^2_{z_j} = \sigma^2/\delta^2_j$ and a Gneiting correlation function with decay parameters $\phi_{s,j}$ and $\phi_{t,j}$ . Given these and the noise-to-spatial variance ratios $\delta^2_j$ , the joint posterior distribution is available in closed form and stvcLMexact() samples from it exactly. With process.type = "independent", each process has its own $(\phi_s, \phi_t, \delta^2)$ , supplied as vectors of length $r$ ; with process.type = "independent.shared", the processes share one.

mod_lm <- stvcLMexact(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",
                      sptParams = list(phi_s = c(2, 2), phi_t = c(4, 4)),
                      noise_sp_ratio = c(0.25, 0.5),
                      n.samples = 1000, loopd = TRUE, verbose = TRUE)

## ----------------------------------------
##  Model description
## ----------------------------------------
## Model fit with 200 observations.
## 
## Number of covariates 3 (including intercept).
## Number of covariates with spatially-temporally varying coefficients 2.
## 
## Using the gneiting-decay spatial-temporal correlation function.
## Process type: independent.
## 
## Priors:
##  beta flat.
##  sigma.sq: flat, proportional to 1/sigma.sq.
## 
## Spatial-temporal process parameters:
##  process 1: phi_s = 2.00, phi_t = 4.00, noise-to-spatial variance ratio = 0.25.
##  process 2: phi_s = 2.00, phi_t = 4.00, noise-to-spatial variance ratio = 0.50.
## 
## Number of posterior samples = 1000.
## 
## LOO-PD calculation method = exact.
## ----------------------------------------

The fixed effects are summarized below; y_gauss was simulated with $\beta = (2, 5, -1)^{ \scriptstyle \top }$ . The samples of the process variances are returned as sigmaSq.z.

summary_beta <- t(apply(mod_lm$samples$beta, 1, function(x) quantile(x, c(0.025, 0.5, 0.975))))
rownames(summary_beta) <- mod_lm$X.names
print(summary_beta)

##                   2.5%        50%      97.5%
## (Intercept)  0.6587955  1.5608960  2.4468595
## x1           4.2039317  4.8940221  5.5610880
## x2          -1.0515911 -0.9563878 -0.8644576

Next, we stack candidate models built from candidate values of the process parameters and noise-to-spatial variance ratios. For process.type = "independent", the candidate values are vectors of length $r$ , supplied as list() entries in candidateModels().

mod.list.lm <- 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.25, 0.5), c(0.5, 1))),
                               "cartesian")

mod_lm_stack <- 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.list.lm,
                            n.samples = 1000, verbose = TRUE)

## 
## STACKING WEIGHTS:
## 
##           | phi_s[1] | phi_s[2] | phi_t[1] | phi_t[2] | noise_sp_ratio[1] | noise_sp_ratio[2] | weight |
## +---------+----------+----------+----------+----------+-------------------+-------------------+--------+
## | Model 1 |         2|         2|         4|         4|               0.25|                0.5| 0.567  |
## | Model 2 |         3|         3|         4|         4|               0.25|                0.5| 0.000  |
## | Model 3 |         2|         2|         8|         8|               0.25|                0.5| 0.116  |
## | Model 4 |         3|         3|         8|         8|               0.25|                0.5| 0.000  |
## | Model 5 |         2|         2|         4|         4|               0.50|                1.0| 0.317  |
## | Model 6 |         3|         3|         4|         4|               0.50|                1.0| 0.000  |
## | Model 7 |         2|         2|         8|         8|               0.50|                1.0| 0.000  |
## | Model 8 |         3|         3|         8|         8|               0.50|                1.0| 0.000  |
## +---------+----------+----------+----------+----------+-------------------+-------------------+--------+

post_lm <- stackedSampler(mod_lm_stack)

The samples of $z$ are stacked by process: rows $1, \ldots, n$ hold $z_1$ and rows $n+1, \ldots, 2n$ hold $z_2$ . We compare their stacked posteriors with the true values.

post_z1_summ <- t(apply(post_lm$z[1:n_train, ], 1,
                        function(x) quantile(x, c(0.025, 0.5, 0.975))))
post_z2_summ <- t(apply(post_lm$z[n_train + 1:n_train, ], 1,
                        function(x) quantile(x, c(0.025, 0.5, 0.975))))

z1_combn <- data.frame(z = dat$z1_true, zL = post_z1_summ[, 1],
                       zM = post_z1_summ[, 2], zU = post_z1_summ[, 3])
z2_combn <- data.frame(z = dat$z2_true, zL = post_z2_summ[, 1],
                       zM = post_z2_summ[, 2], zU = post_z2_summ[, 3])

library(ggplot2)
plot_z1_summ <- ggplot(data = z1_combn, aes(x = z)) +
  geom_errorbar(aes(ymin = zL, ymax = zU), alpha = 0.5, color = "skyblue") +
  geom_point(aes(y = zM), size = 0.5, color = "darkblue", alpha = 0.5) +
  geom_abline(slope = 1, intercept = 0, color = "red", linetype = "solid") +
  xlab("True z1") + ylab("Stacked posterior of z1") + theme_bw() +
  theme(panel.grid = element_blank(), aspect.ratio = 1)

plot_z2_summ <- ggplot(data = z2_combn, aes(x = z)) +
  geom_errorbar(aes(ymin = zL, ymax = zU), alpha = 0.5, color = "skyblue") +
  geom_point(aes(y = zM), size = 0.5, color = "darkblue", alpha = 0.5) +
  geom_abline(slope = 1, intercept = 0, color = "red", linetype = "solid") +
  xlab("True z2") + ylab("Stacked posterior of z2") + theme_bw() +
  theme(panel.grid = element_blank(), aspect.ratio = 1)

plot_z1_summ + plot_z2_summ

Stacked posterior of the varying coefficients against their true values.

Bayesian non-Gaussian spatially-temporally varying coefficient models

We now analyze the Poisson counts y_pois. Given family = "poisson", the formula y_pois ~ x1 + x2 + (x1) corresponds to the spatial-temporal generalized linear model $y(\ell) \sim \mathsf{Poisson}(\lambda(\ell)), \quad \log \lambda(\ell) = \beta_0 + \beta_1 x_1(\ell) + \beta_2 x_2(\ell) + z_1(\ell) + x_1(\ell) z_2(\ell)\;.$ The spatially-temporally varying coefficients $z(\ell) = (z_1(\ell), z_2(\ell))^{{ \scriptstyle \top }}$ is a multivariate Gaussian process, and we pursue the following specifications for $z(\ell)$ - independent process, independent process with shared parameters, and a multivariate process. See “Technical Overview” for more details.

Using fixed hyperparameters

We use the function stvcGLMexact() to fit spatially-temporally varying coefficient generalized linear models. In the following code snippets, we demonstrate the uasge of the argument process.Type to implement different variations of spatial-temporal process specifications for the varying coefficients.

Independent processes

In this case, since there are two independent processes $z_1(\ell)$ and $z_2(\ell)$ the candidate values of the spatial-temporal process parameters sptParams is a list with tags phi_s and phi_t, with each tag being of length 2. Here, the scale parameter $\sigma = (\sigma^2_{z1}, \sigma^2_{z2})^{{ \scriptstyle \top }}$ has dimension 2.

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",
                     priors = list(nu.beta = 5, nu.z = 5),
                     sptParams = list(phi_s = c(3, 3), phi_t = c(4, 4)),
                     verbose = FALSE, n.samples = 500)

## Some priors were not supplied. Using defaults.

Posterior samples of the scale parameters can be recovered by running recoverGLMscale() on mod1.

mod1 <- recoverGLMscale(mod1)

We visualize the posterior distributions of the scale parameters as follows.

post_scale_df <- data.frame(value = sqrt(c(mod1$samples$z.scale[1, ], mod1$samples$z.scale[2, ])),
                            group = factor(rep(c("sigma.z1", "sigma.z2"),
                                    each = length(mod1$samples$z.scale[1, ]))))
library(ggplot2)
ggplot(post_scale_df, aes(x = value)) +
  geom_density(fill = "lightblue", alpha = 0.6) +
  facet_wrap(~ group, scales = "free") + labs(x = "", y = "Density") +
  theme_bw() + theme(panel.background = element_blank(),
                     panel.grid = element_blank(), aspect.ratio = 1)

Posterior distributions of the scale parameters.

Independent shared processes

In this case, the processes $z_1(\ell)$ and $z_2(\ell)$ are independent but share a common covariance matrix. Hence, sptParams is a list with tags phi_s and phi_t, with each tag being of length 1. Here, the scale parameter $\sigma = \sigma_z^2$ is 1-dimensional.

mod2 <- 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.shared",
                     priors = list(nu.beta = 5, nu.z = 5),
                     sptParams = list(phi_s = 3, phi_t = 4),
                     verbose = FALSE, n.samples = 500)

## Some priors were not supplied. Using defaults.

Posterior samples of the scale parameters can be recovered by running recoverGLMscale() on mod2.

mod2 <- recoverGLMscale(mod2)

We visualize the posterior distributions of the scale parameters as follows.

post_scale_df <- data.frame(value = sqrt(mod2$samples$z.scale),
                            group = factor(rep(c("sigma.z"),
                                               each = length(mod2$samples$z.scale))))
ggplot(post_scale_df, aes(x = value)) +
  geom_density(fill = "lightblue", alpha = 0.6) +
  facet_wrap(~ group, scales = "free") + labs(x = "", y = "Density") +
  theme_bw() + theme(panel.background = element_blank(),
                     panel.grid = element_blank(), aspect.ratio = 1)

Posterior distributions of the scale parameters.

Multivariate processes

In this case, $z(\ell) = (z_1(\ell), z_2(\ell))^{{ \scriptstyle \top }}$ is a 2-dimensional Gaussian process with covariance matrix $\Sigma$ . Further, we put an inverse-Wishart prior on $\Sigma$ , which can be specified through the priors argument. If not supplied, uses the default $\mathrm{IW}(\nu_z + 2r, I_r)$ , where $r = 2$ is the dimension of the multivariate process. Here, sptParams is a list with tags phi_s and phi_t, with each tag being of length 1, and the scale parameter $\sigma = \Sigma$ is an $2 \times 2$ matrix.

mod3 <- 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 = "multivariate",
                     priors = list(nu.beta = 5, nu.z = 5),
                     sptParams = list(phi_s = 3, phi_t = 4),
                     verbose = FALSE, n.samples = 500)

## Some priors were not supplied. Using defaults.

Posterior samples of the scale parameters can be recovered by running recoverGLMscale() on mod3.

mod3 <- recoverGLMscale(mod3)

We visualize the posterior distribution of the scale matrix $\Sigma$ as follows.

post_scale_z <- mod3$samples$z.scale

r <- sqrt(dim(post_scale_z)[1])
# Function to get (i,j) index from row number (column-major)
get_indices <- function(k, r) {
  j <- ((k - 1) %/% r) + 1
  i <- ((k - 1) %% r) + 1
  c(i, j)
}

# Generate plots into a matrix
plot_matrix <- matrix(vector("list", r * r), nrow = r, ncol = r)
for (k in 1:(r^2)) {
  ij <- get_indices(k, r)
  i <- ij[1]
  j <- ij[2]

  if (i >= j) {
    df <- data.frame(value = post_scale_z[k, ])
    p <- ggplot(df, aes(x = value)) +
      geom_density(fill = "lightblue", alpha = 0.7) +
      theme_bw(base_size = 9) +
      labs(title = bquote(Sigma[.(i) * .(j)])) +
      theme(axis.title = element_blank(), axis.text = element_text(size = 6),
        plot.title = element_text(size = 9, hjust = 0.5),
        panel.grid = element_blank(), aspect.ratio = 1)
  } else {
    p <- ggplot() + theme_void()
  }

  plot_matrix[j, i] <- list(p)
}

library(patchwork)
# Assemble with patchwork
final_plot <- wrap_plots(plot_matrix, nrow = r)
final_plot

Posterior distributions of elements of the scale matrix.

Using predictive stacking

For implementing predictive stacking for spatially-temporally varying models, we offer a helper function candidateModels() to create a collection of candidate models. The grid of candidate values can be combined either using a Cartesian product or a simple element-by-element combination. We demonstrate stacking based on the multivariate spatial-temporal process model.

Step 1. Create candidate models.

mod.list <- candidateModels(list(
  phi_s = list(2, 3, 4),
  phi_t = list(2, 4, 8),
  boundary = c(0.5, 0.75)), "cartesian")

Step 2. Run stvcGLMstack().

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.list,
                     loopd.controls = list(method = "CV", CV.K = 10, nMC = 500),
                     n.samples = 1000)

## 
## STACKING WEIGHTS:
## 
##            | phi_s | phi_t | boundary | weight |
## +----------+-------+-------+----------+--------+
## | Model 1  |      2|      2|      0.50| 0.000  |
## | Model 2  |      3|      2|      0.50| 0.000  |
## | Model 3  |      4|      2|      0.50| 0.000  |
## | Model 4  |      2|      4|      0.50| 0.000  |
## | Model 5  |      3|      4|      0.50| 0.000  |
## | Model 6  |      4|      4|      0.50| 0.000  |
## | Model 7  |      2|      8|      0.50| 0.000  |
## | Model 8  |      3|      8|      0.50| 0.000  |
## | Model 9  |      4|      8|      0.50| 0.000  |
## | Model 10 |      2|      2|      0.75| 0.074  |
## | Model 11 |      3|      2|      0.75| 0.000  |
## | Model 12 |      4|      2|      0.75| 0.000  |
## | Model 13 |      2|      4|      0.75| 0.926  |
## | Model 14 |      3|      4|      0.75| 0.000  |
## | Model 15 |      4|      4|      0.75| 0.000  |
## | Model 16 |      2|      8|      0.75| 0.000  |
## | Model 17 |      3|      8|      0.75| 0.000  |
## | Model 18 |      4|      8|      0.75| 0.000  |
## +----------+-------+-------+----------+--------+

Step 3. Recover posterior samples of the scale parameters.

mod1 <- recoverGLMscale(mod1)

Step 4. Sample from the stacked posterior distribution.

post_samps <- stackedSampler(mod1)

Now, we analyze the posterior distribution of the latent process as obtained from the stacked posterior.

post_z <- post_samps$z

post_z1_summ <- t(apply(post_z[1:n_train,], 1,
                        function(x) quantile(x, c(0.025, 0.5, 0.975))))
post_z2_summ <- t(apply(post_z[n_train + 1:n_train,], 1,
                        function(x) quantile(x, c(0.025, 0.5, 0.975))))

z1_combn <- data.frame(z = dat$z1_true, zL = post_z1_summ[, 1],
                       zM = post_z1_summ[, 2], zU = post_z1_summ[, 3])
z2_combn <- data.frame(z = dat$z2_true, zL = post_z2_summ[, 1],
                       zM = post_z2_summ[, 2], zU = post_z2_summ[, 3])

plot_z1_summ <- ggplot(data = z1_combn, aes(x = z)) +
  geom_errorbar(aes(ymin = zL, ymax = zU), alpha = 0.5, color = "skyblue") +
  geom_point(aes(y = zM), size = 0.5, color = "darkblue", alpha = 0.5) +
  geom_abline(slope = 1, intercept = 0, color = "red", linetype = "solid") +
  xlab("True z1") + ylab("Posterior of z1") + theme_bw() +
  theme(panel.grid = element_blank(), aspect.ratio = 1)

plot_z2_summ <- ggplot(data = z2_combn, aes(x = z)) +
  geom_errorbar(aes(ymin = zL, ymax = zU), alpha = 0.5, color = "skyblue") +
  geom_point(aes(y = zM), size = 0.5, color = "darkblue", alpha = 0.5) +
  geom_abline(slope = 1, intercept = 0, color = "red", linetype = "solid") +
  xlab("True z2") + ylab("Posterior of z2") + theme_bw() +
  theme(panel.grid = element_blank(), aspect.ratio = 1)

plot_z1_summ + plot_z2_summ

Next, we analyze the posterior distribution of the scale matrix that models the inter-process dependence structure.

post_scale_z <- post_samps$z.scale
r <- sqrt(dim(post_scale_z)[1])
# Generate plots into a matrix
plot_matrix <- matrix(vector("list", r * r), nrow = r, ncol = r)
for (k in 1:(r^2)) {
  ij <- get_indices(k, r)
  i <- ij[1]
  j <- ij[2]

  if (i >= j) {
    df <- data.frame(value = post_scale_z[k, ])
    p <- ggplot(df, aes(x = value)) +
      geom_density(fill = "lightblue", alpha = 0.7) +
      theme_bw(base_size = 9) +
      labs(title = bquote(Sigma[.(i) * .(j)])) +
      theme(axis.title = element_blank(), axis.text = element_text(size = 6),
        plot.title = element_text(size = 9, hjust = 0.5),
        panel.grid = element_blank(), aspect.ratio = 1)
  } else {
    p <- ggplot() + theme_void()
  }

  plot_matrix[j, i] <- list(p)
}

# Assemble with patchwork
final_plot <- wrap_plots(plot_matrix, nrow = r)
final_plot

Stacked posterior distribution of the elements of the inter-process covariance matrix.