Spatial and Temporal Models
Gilles Colling
2026-10-13
Source:vignettes/spatial-temporal.Rmd
spatial-temporal.RmdInstallation
install.packages("pak") # development version
pak::pak("gcol33/tulpaRatio")Overview
This vignette covers spatial and temporal extensions in tulpaRatio. These allow borrowing strength across neighboring locations or time points, improving estimates in data-sparse regions.
Spatial Models
tulpaRatio supports two spatial priors for areal (discrete region) data:
| Function | Prior | Use Case |
|---|---|---|
spatial_car() |
Intrinsic CAR (ICAR) | Simple spatial smoothing |
spatial_bym2() |
BYM2 | Interpretable spatial fraction |
Both require an adjacency matrix defining which regions are neighbors.
Creating Adjacency Matrices
From Coordinates (Grid Data)
# Camera trap grid
n <- 25
coords <- expand.grid(x = 1:5, y = 1:5)
# Queen contiguity (8 neighbors)
adj <- matrix(0, n, n)
for (i in 1:n) {
for (j in 1:n) {
if (i != j) {
dist <- sqrt((coords$x[i] - coords$x[j])^2 +
(coords$y[i] - coords$y[j])^2)
if (dist <= sqrt(2) + 0.01) adj[i, j] <- 1
}
}
}From sf Objects
library(sf)
#> Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.4.0; sf_use_s2() is TRUE
# North Carolina counties, shipped with sf
nc <- st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE)
nc$county <- seq_len(nrow(nc))
# Queen contiguity: counties sharing a border or a corner
adj_nc <- 1 * st_touches(nc, sparse = FALSE)ICAR Model: spatial_car()
The intrinsic conditional autoregressive (ICAR) prior encourages neighboring regions to have similar effects:
where denotes neighbors of region .
# Four surveys per grid cell; the spatial field enters catch and effort
set.seed(1)
phi <- 0.4 * sin(coords$x) + 0.4 * cos(coords$y)
df <- data.frame(region = rep(1:n, each = 4), x = rnorm(4 * n))
df$effort <- rgamma(nrow(df), shape = 5, rate = 1 / exp(phi[df$region]))
df$count <- rpois(nrow(df), exp(3 + 0.3 * df$x + phi[df$region]))
fit <- tratio(
count | effort ~ x,
spatial = spatial_car(
adjacency = adj,
level = "group",
group_var = "region"
),
data = df,
family = ratiod_poisson_gamma(),
mode = "hmc",
control = list(seed = 1)
)
summary(fit)
#> ratio model summary
#> ===================
#>
#> Inference: Exact (Tier 1) via hmc
#> Family: poisson_gamma
#>
#> Fixed effects (numerator):
#> parameter mean sd 2.5% 97.5% rhat ess_bulk
#> num:(Intercept) 2.927 0.025 2.878 2.974 1.000 2087
#> num:x 0.292 0.026 0.241 0.342 1.002 1755
#>
#> Fixed effects (denominator):
#> parameter mean sd 2.5% 97.5% rhat ess_bulk
#> denom:(Intercept) 1.453 0.047 1.361 1.549 1.000 1604
#> denom:x 0.043 0.055 -0.065 0.150 1.003 1301
#>
#> Variance components:
#> parameter mean sd 2.5% 97.5% rhat ess_bulk
#> shape 4.425 0.608 3.324 5.701 1.002 1474
#> tau_spatial 3.238 1.005 1.602 5.497 1.003 1059
#>
#> Diagnostics:
#> Divergences: 0
#> Avg. acceptance: 0.887
#> Chains: 4
#> Iterations: 2000
#> Warmup: 1000BYM2 Model: spatial_bym2()
The BYM2 reparameterization (Riebler et al., 2016) separates spatial structure from unstructured heterogeneity:
where:
- is the scaled ICAR component
- is unstructured noise
- is the proportion of variance that is spatial
# Sudden infant deaths among live births, 1974-78
nc_df <- st_drop_geometry(nc)
nc_df$nonwhite <- nc_df$NWBIR74 / nc_df$BIR74
fit_bym2 <- tratio(
SID74 | BIR74 ~ nonwhite,
spatial = spatial_bym2(
adjacency = adj_nc,
level = "group",
group_var = "county"
),
data = nc_df,
family = ratiod_binomial(),
mode = "hmc",
control = list(seed = 1)
)
summary(fit_bym2)
#> ratio model summary
#> ===================
#>
#> Inference: Exact (Tier 1) via hmc
#> Family: binomial_fixed
#>
#> Fixed effects (numerator):
#> parameter mean sd 2.5% 97.5% rhat ess_bulk
#> num:(Intercept) -6.893 0.126 -7.141 -6.649 1.013 548
#> num:nonwhite 2.002 0.329 1.354 2.635 1.016 339
#>
#> Variance components:
#> parameter mean sd 2.5% 97.5% rhat ess_bulk
#> rho_spatial 0.457 0.271 0.026 0.953 1.004 487
#> sigma_s_spatial 0.178 0.086 0.038 0.366 1.002 577
#> sigma_spatial 0.275 0.072 0.138 0.420 1.004 740
#> sigma_u_spatial 0.191 0.072 0.055 0.332 1.006 570
#>
#> Diagnostics:
#> Divergences: 0
#> Avg. acceptance: 0.916
#> Chains: 4
#> Iterations: 2000
#> Warmup: 1000
#> R-hat > 1.01 for 2 parameter(s), largest 1.016 (num:nonwhite).
#> Bulk ESS < 400 for 1 parameter(s), smallest 339 (num:nonwhite).
#> Run longer chains or reparameterise before using these estimates.The parameter is directly interpretable: means 80% of the random effect variance is spatially structured.
Extracting Spatial Effects
# County-level posterior summaries
county_rate <- as.data.frame(ratio(fit_bym2, by = "county", summary = TRUE))
# Map visualization
library(ggplot2)
nc_map <- merge(nc, county_rate, by.x = "county", by.y = "group")
ggplot(nc_map) +
geom_sf(aes(fill = mean)) +
scale_fill_viridis_c(name = "SIDS rate")
Temporal Models
tulpaRatio supports temporal priors for time series structure:
| Function | Prior | Use Case |
|---|---|---|
temporal_rw1(), temporal_rw2()
|
Random walk | Smooth trends |
temporal_ar1() |
AR(1) | Mean-reverting dynamics |
Random Walk: temporal_rw1(),
temporal_rw2()
First-order random walk (RW1):
Second-order random walk (RW2) for smoother trends:
# Twelve sites surveyed once a year for twelve years, in a random month
set.seed(2)
df_t <- expand.grid(site = factor(1:12), year = 2010:2021)
df_t$month <- sample(1:12, nrow(df_t), replace = TRUE)
df_t$x <- rnorm(nrow(df_t))
gamma_t <- 0.3 * sin((df_t$year - 2010) / 2) + 0.2 * cos(2 * pi * df_t$month / 12) +
rnorm(12, 0, 0.3)[df_t$site]
df_t$effort <- rgamma(nrow(df_t), shape = 5, rate = 1 / exp(gamma_t))
df_t$count <- rpois(nrow(df_t), exp(3 + 0.2 * df_t$x + gamma_t))
# First-order random walk
fit_rw1 <- tratio(
count | effort ~ x + (1 | site),
temporal = temporal_rw1(time_var = "year"),
data = df_t,
family = ratiod_poisson_gamma(),
control = list(seed = 1)
)
# Second-order for smoother trends
fit_rw2 <- tratio(
count | effort ~ x,
temporal = temporal_rw2(time_var = "year"),
data = df_t,
family = ratiod_poisson_gamma(),
control = list(seed = 1)
)Cyclic Random Walk
For seasonal patterns that wrap around (e.g., month 12 neighbors month 1):
fit <- tratio(
count | effort ~ x,
temporal = temporal_rw1(
time_var = "month",
cyclic = TRUE
),
data = df_t,
family = ratiod_poisson_gamma(),
control = list(seed = 1)
)AR(1): temporal_ar1()
Autoregressive order 1:
The correlation parameter controls persistence.
fit <- tratio(
count | effort ~ x,
temporal = temporal_ar1(
time_var = "year"
),
data = df_t,
family = ratiod_poisson_gamma(),
control = list(seed = 1)
)Group-Specific Temporal Effects
By default, temporal effects are shared across all groups. For
group-specific trends, group_var gives each group its own
walk. Each walk sums to zero over time, so it carries the group’s trend,
and a group intercept carries the group’s level:
# Each site follows its own random-walk trend
df_g <- df_t
site_trend <- apply(matrix(rnorm(12 * 12, 0, 0.15), 12), 2, cumsum)
eta_g <- rnorm(12, 0, 0.3)[df_g$site] +
site_trend[cbind(df_g$year - 2009, as.integer(df_g$site))]
df_g$effort <- rgamma(nrow(df_g), shape = 5, rate = 1 / exp(eta_g))
df_g$count <- rpois(nrow(df_g), exp(3 + eta_g))
fit <- tratio(
count | effort ~ 1 + (1 | site),
temporal = temporal_rw1(
time_var = "year",
group_var = "site"
),
data = df_g,
family = ratiod_poisson_gamma(),
control = list(seed = 1)
)
#> Warning: 1 divergent transition(s) after warmup; raise `adapt_delta` or
#> reparameterise.Extracting Temporal Effects
# Year-level effects
temporal_effects <- temporal(fit_rw1, summary = TRUE)
temporal_effects
#> time_idx time mean sd q2.5 q50 q97.5
#> 1 1 2010 -0.06178548 0.05178706 -0.16689648 -0.06240901 0.038603395
#> 2 2 2011 0.14243254 0.04260754 0.05924424 0.14169815 0.227371519
#> 3 3 2012 0.15647954 0.04406724 0.07049364 0.15650034 0.242298943
#> 4 4 2013 0.23149480 0.04029654 0.15359599 0.22977105 0.311901868
#> 5 5 2014 0.29412744 0.04154833 0.21304317 0.29472991 0.375692127
#> 6 6 2015 0.18911816 0.04121639 0.10960488 0.18883433 0.267593275
#> 7 7 2016 0.06790121 0.04542897 -0.02017534 0.06699919 0.159292552
#> 8 8 2017 -0.08568612 0.04604968 -0.17636793 -0.08607211 0.004921489
#> 9 9 2018 -0.17981113 0.04713092 -0.27205409 -0.17919335 -0.086293365
#> 10 10 2019 -0.25760129 0.04964322 -0.35906271 -0.25630280 -0.165184251
#> 11 11 2020 -0.25406534 0.04960702 -0.34767276 -0.25441512 -0.154887377
#> 12 12 2021 -0.24260434 0.05135490 -0.34242301 -0.24178457 -0.140325552
# Plot trend
plot(temporal(fit_rw1))
Combining Spatial and Temporal
For spatiotemporal data, include both:
# The 25 grid cells surveyed twice a year for four years
set.seed(3)
df_st <- expand.grid(region = 1:n, year = 2015:2018, survey = 1:2)
df_st$x <- rnorm(nrow(df_st))
df_st$lon <- coords$x[df_st$region]
df_st$lat <- coords$y[df_st$region]
delta <- matrix(rnorm(n * 4, sd = 0.5), n, 4)
# Spatial field, a random walk over the years, and a region-by-year term
year_eff <- cumsum(rnorm(4, 0, 0.3))
eta_st <- phi[df_st$region] + year_eff[df_st$year - 2014] +
delta[cbind(df_st$region, df_st$year - 2014)]
df_st$effort <- rgamma(nrow(df_st), shape = 5, rate = 1 / exp(eta_st))
df_st$count <- rpois(nrow(df_st), exp(3 + 0.3 * df_st$x + eta_st))
fit <- tratio(
count | effort ~ x,
spatial = spatial_bym2(adj, level = "group", group_var = "region"),
temporal = temporal_rw1(time_var = "year"),
data = df_st,
family = ratiod_poisson_gamma(),
mode = "hmc",
control = list(seed = 1)
)This fits additive spatial and temporal effects. For interaction effects, see the next section.
Spatiotemporal Interaction
When spatial patterns change over time (or temporal trends differ across space), additive models may be insufficient. tulpaRatio supports spatiotemporal interaction effects following Knorr-Held (2000).
Interaction Types
| Type | Description | Parameters |
|---|---|---|
| Type I | Unstructured (IID) | |
| Type II | Structured time at each location | |
| Type III | Structured space at each time | |
| Type IV | Fully structured (Kronecker) | |
| Separable | For GP spatial/temporal |
Type I: Unstructured Interaction
The simplest form assumes independent random effects for each space-time combination:
fit_type1 <- tratio(
count | effort ~ x,
data = df_st,
family = ratiod_poisson_gamma(),
spatiotemporal = spatiotemporal(
spatial = spatial_car(adj, level = "group", group_var = "region"),
temporal = temporal_rw1("year"),
type = "I"
),
control = list(seed = 1)
)Type II: Structured Time at Each Location
Each location has its own temporal random walk:
This captures location-specific temporal trends.
fit_type2 <- tratio(
count | effort ~ x,
data = df_st,
family = ratiod_poisson_gamma(),
spatiotemporal = spatiotemporal(
spatial = spatial_car(adj, level = "group", group_var = "region"),
temporal = temporal_rw1("year"),
type = "II"
),
control = list(seed = 1)
)Type III: Structured Space at Each Time Point
Each time point has its own spatial field:
This captures time-specific spatial patterns.
fit_type3 <- tratio(
count | effort ~ x,
data = df_st,
family = ratiod_poisson_gamma(),
spatiotemporal = spatiotemporal(
spatial = spatial_car(adj, level = "group", group_var = "region"),
temporal = temporal_rw1("year"),
type = "III"
),
control = list(seed = 1)
)Type IV: Fully Structured (Kronecker)
Precision is the Kronecker product of spatial and temporal precision:
This is the most constrained and often most appropriate model.
fit_type4 <- tratio(
count | effort ~ x,
data = df_st,
family = ratiod_poisson_gamma(),
spatiotemporal = spatiotemporal(
spatial = spatial_car(adj, level = "group", group_var = "region"),
temporal = temporal_rw1("year"),
type = "IV"
),
control = list(seed = 1)
)Non-Separable Spatiotemporal GP
For continuous space-time data, type = "nonsep_gp" puts
a covariance on the space-time grid rather than a GMRF stencil on it.
The interaction carries its own coords, so an areal main
effect pairs with a continuous interaction over the same units:
fit <- tratio(
count | effort ~ x,
data = df_st,
family = ratiod_poisson_gamma(),
spatiotemporal = spatiotemporal(
spatial = spatial_car(adj, level = "group", group_var = "region"),
temporal = temporal_rw1("year"),
type = "nonsep_gp",
nonsep_type = "gneiting",
coords = ~ lon + lat,
parameterization = "collapsed"
),
control = list(iter = 1000, chains = 2, seed = 1)
)A GP interaction has one value per space-time cell, and with the
field sampled alongside its precision and two ranges NUTS needs a small
step size and deep trees. parameterization = "collapsed"
integrates the interaction out with a Laplace approximation inside HMC,
so the sampler moves only its three hyperparameters beside the other
blocks, and the interaction is drawn from its Gaussian approximation at
each saved iteration. On this model the two chains above take about five
minutes on one core each and reach rhat of 1.00 for every
hyperparameter, where the centred interaction takes about 390 leapfrog
steps per draw.
Non-separability types:
-
"product": the separable covariance,C_s(h) C_t(u).type = "separable"is this arm. -
"gneiting": Gneiting (2002) class
An additive covariance C_s(h) + C_t(u) is not among
them. On a complete S x T grid its rank is
S + T - 1, and the directions it drops are exactly the
space-time interactions the block exists to carry.
Both ranges are estimated on their own log scale, alongside the
interaction’s precision. The field is approximated by an NNGP over the
grid, with nn neighbours per cell chosen on a joint
space-time metric.
Extracting Spatiotemporal Effects
# Summary format
st_effects <- spatiotemporal_effects(fit_type4, format = "summary")
head(st_effects)
#> Spatiotemporal interaction posterior summary: 25 locations x 4 times, 4000 posterior draws
#> s t mean q2.5 q50 q97.5 sd
#> 1 1 1 -7.053193e-01 -0.9862563 -0.709673205 -0.4178550 0.14563858
#> 2 2 1 1.074187e+00 0.9386545 1.075267924 1.2030313 0.06779440
#> 3 3 1 -4.182706e-01 -0.6765575 -0.416632197 -0.1626701 0.13352410
#> 4 4 1 4.928829e-02 -0.1480003 0.046848155 0.2477211 0.09918574
#> 5 5 1 -9.941303e-05 -0.2158890 -0.001129724 0.2120820 0.10742076
#> 6 6 1 -3.938569e-01 -0.5985219 -0.390020116 -0.2009649 0.09851789
# Array format (S x T x draws)
st_array <- spatiotemporal_effects(fit_type4, format = "array")
dim(st_array)
#> [1] 25 4 4000
# Visualization
plot(st_effects, type = "heatmap")
plot(st_effects, type = "time_series")
Model Selection
Compare interaction types using LOO-CV:
ratiod_compare(fit_type1, fit_type4)
#> Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.
#> Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.
#> model elpd_diff se_diff p_worse diag_diff diag_elpd
#> fit_type1 0.0 0.0 NA 174 k_psis > 0.7
#> fit_type4 -107.4 15.2 1.00 185 k_psis > 0.7
#>
#> Diagnostic flags present.
#> See ?`loo-glossary` (sections `diag_diff` and `diag_elpd`)
#> or https://mc-stan.org/loo/reference/loo-glossary.html.Each space-time cell carries its own interaction parameter, informed here by two surveys, so leaving one survey out shifts that parameter strongly and PSIS-LOO flags most observations with Pareto ; the ranking is a guide, and K-fold cross-validation over held-out cells gives the reliable comparison.
Prior Specification
Spatial Priors
The sigma prior is the half-Cauchy on the total standard
deviation of a BYM2 field (and on random-effect standard deviations);
its scale sets how large a spatial effect is plausible before the data
are seen. An ICAR field’s precision takes a Gamma(1, 0.01) prior.
fit <- tratio(
count | effort ~ x,
data = df,
family = ratiod_poisson_gamma(),
spatial = spatial_bym2(adj, level = "group", group_var = "region"),
priors = ratiod_priors(sigma = prior_half_cauchy(0.5)),
control = list(seed = 1)
)Temporal Priors
fit <- tratio(
count | effort ~ x,
data = df_t,
family = ratiod_poisson_gamma(),
temporal = temporal_ar1(time_var = "year"),
priors = ratiod_priors(
rho_temporal = prior_beta(5, 2) # Beta on (rho + 1) / 2, mode near 0.6
),
control = list(seed = 1)
)See Also
-
vignette("getting-started")- Basic usage -
vignette("workflows")- Complete analysis examples -
?spatial_car,?spatial_bym2- Function documentation -
?temporal_rw1,?temporal_rw2,?temporal_ar1- Function documentation