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.
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.
spatial_car()The intrinsic conditional autoregressive (ICAR) prior encourages neighboring regions to have similar effects:
\[ \phi_i | \phi_{-i} \sim \text{Normal}\left(\frac{1}{n_i}\sum_{j \sim i} \phi_j, \frac{\sigma^2}{n_i}\right) \]
where \(j \sim i\) denotes neighbors of region \(i\).
# 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)adjacency: Square binary adjacency matrixlevel: "group" (one effect per group) or
"obs" (one per observation)group_var: Variable name matching the grouping factor
(required if level = "group")spatial_bym2()The BYM2 reparameterization (Riebler et al., 2016) separates spatial structure from unstructured heterogeneity:
\[ \theta_i = \sigma \left(\sqrt{\rho} \cdot \phi_i^* + \sqrt{1-\rho} \cdot v_i\right) \]
where:
# 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)The \(\rho\) parameter is directly interpretable: \(\rho = 0.8\) means 80% of the random effect variance is spatially structured.
# 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")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 |
temporal_rw1(),
temporal_rw2()First-order random walk (RW1):
\[ \gamma_t | \gamma_{t-1} \sim \text{Normal}(\gamma_{t-1}, \sigma_\gamma^2) \]
Second-order random walk (RW2) for smoother trends:
\[ \gamma_t | \gamma_{t-1}, \gamma_{t-2} \sim \text{Normal}(2\gamma_{t-1} - \gamma_{t-2}, \sigma_\gamma^2) \]
# 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)
)temporal_ar1()Autoregressive order 1:
\[ \gamma_t = \rho \cdot \gamma_{t-1} + \epsilon_t, \quad \epsilon_t \sim \text{Normal}(0, \sigma_\gamma^2) \]
The correlation parameter \(\rho \in (-1, 1)\) controls persistence.
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)
)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.
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).
| Type | Description | Parameters |
|---|---|---|
| Type I | Unstructured (IID) | \(S \times T\) |
| Type II | Structured time at each location | \(S \times T\) |
| Type III | Structured space at each time | \(S \times T\) |
| Type IV | Fully structured (Kronecker) | \(S \times T\) |
| Separable | For GP spatial/temporal | \(S \times T\) |
The simplest form assumes independent random effects for each space-time combination:
\[ \delta_{st} \stackrel{iid}{\sim} \text{Normal}(0, \sigma_\delta^2) \]
Each location has its own temporal random walk:
\[ \delta_{\cdot t}^{(s)} \sim \text{RW}(\sigma^2) \]
This captures location-specific temporal trends.
Each time point has its own spatial field:
\[ \delta_{s \cdot}^{(t)} \sim \text{ICAR}(\tau) \]
This captures time-specific spatial patterns.
Precision is the Kronecker product of spatial and temporal precision:
\[ Q_\delta = Q_s \otimes Q_t \]
This is the most constrained and often most appropriate model.
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, at about 20 leapfrog steps per draw where the
centred interaction takes about 390. The spatial range of the
interaction can settle beyond the extent of the grid, where the field is
nearly constant in space and the data say little about its precision;
rhat and the divergent transitions of tau_st
and phi_st_space show whether the chains did.
Non-separability types:
"product": the separable covariance,
C_s(h) C_t(u). type = "separable" is this
arm."gneiting": Gneiting (2002) classAn 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.
# Summary format
st_effects <- spatiotemporal_effects(fit_type4, format = "summary")
head(st_effects)
# Array format (S x T x draws)
st_array <- spatiotemporal_effects(fit_type4, format = "array")
dim(st_array)
# Visualization
plot(st_effects, type = "heatmap")
plot(st_effects, type = "time_series")Compare interaction types using LOO-CV:
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 \(k > 0.67\); the ranking is a guide, and K-fold cross-validation over held-out cells gives the reliable comparison.
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.
vignette("getting-started") - Basic usagevignette("workflows") - Complete analysis examples?spatial_car, ?spatial_bym2 - Function
documentation?temporal_rw1, ?temporal_rw2,
?temporal_ar1 - Function documentation