---
title: "Spatial and Temporal Models"
author: "Gilles Colling"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Spatial and Temporal Models}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5,
  eval = identical(Sys.getenv("NOT_CRAN"), "true") ||
    identical(Sys.getenv("IN_PKGDOWN"), "true")
)
```

## Installation

```{r, eval = FALSE}
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.

```{r}
library(tulpaRatio)
```


# 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)

```{r}
# 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

```{r}
library(sf)

# 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:

$$
\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$.

```{r}
# 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)
```

### Parameters

- `adjacency`: Square binary adjacency matrix
- `level`: `"group"` (one effect per group) or `"obs"` (one per observation)
- `group_var`: Variable name matching the grouping factor (required if `level = "group"`)

## BYM2 Model: `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:

- $\phi_i^*$ is the scaled ICAR component
- $v_i \sim \text{Normal}(0, 1)$ is unstructured noise
- $\rho \in [0, 1]$ is the proportion of variance that is spatial

```{r}
# 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.

## Extracting Spatial Effects

```{r}
# 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):

$$
\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)
$$

```{r}
# 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):

```{r}
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:

$$
\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.

```{r}
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:

```{r}
# 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)
)
```

## Extracting Temporal Effects

```{r}
# Year-level effects
temporal_effects <- temporal(fit_rw1, summary = TRUE)
temporal_effects

# Plot trend
plot(temporal(fit_rw1))
```


# Combining Spatial and Temporal

For spatiotemporal data, include both:

```{r}
# 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) | $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$ |

### Type I: Unstructured Interaction

The simplest form assumes independent random effects for each space-time combination:

$$
\delta_{st} \stackrel{iid}{\sim} \text{Normal}(0, \sigma_\delta^2)
$$

```{r}
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:

$$
\delta_{\cdot t}^{(s)} \sim \text{RW}(\sigma^2)
$$

This captures location-specific temporal trends.

```{r}
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:

$$
\delta_{s \cdot}^{(t)} \sim \text{ICAR}(\tau)
$$

This captures time-specific spatial patterns.

```{r}
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:

$$
Q_\delta = Q_s \otimes Q_t
$$

This is the most constrained and often most appropriate model.

```{r}
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:

```{r, eval = FALSE}
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) 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

```{r}
# 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")
```

## Model Selection

Compare interaction types using LOO-CV:

```{r}
ratiod_compare(fit_type1, fit_type4)
```

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.


## 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.

```{r}
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

```{r}
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