---
title: "3. Spatio-temporal disease mapping"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{3. Spatio-temporal disease mapping}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE, purl = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", eval = FALSE)
```

When the same regions are observed at several times, add `time =` and `SDALGCP2`
fits a separable space–time model. This tutorial is self-contained.

## The model

For region \(i\) at time \(t\), counts \(Y_{it}\) with offset \(m_{it}\):
\[
Y_{it}\mid S \sim \mathrm{Poisson}\!\big(m_{it}\,e^{\eta_{it}}\big),
\qquad \eta_{it}=d_{it}^\top\beta + S_{it},
\]
with a **separable** space–time covariance
\[
\mathrm{Cov}(S_{it},S_{js}) \;=\; \sigma^2\,R_s(\lVert\cdot\rVert;\phi)_{ij}\;R_t(|t-s|;\nu)_{ts},
\]
i.e. \(\mathrm{Cov}(\mathrm{vec}\,S)=\sigma^2\,(R_t(\nu)\otimes R_s(\phi))\), where
\(R_s\) is the aggregated spatial correlation (range \(\phi\)) and \(R_t\) a temporal
Matérn correlation (range \(\nu\)). `SDALGCP2` never forms the
\((N\!\cdot\!T)\times(N\!\cdot\!T)\) covariance — it uses the Kronecker identities
\(x^\top(R_t\otimes R_s)^{-1}x=\mathrm{tr}(R_s^{-1}M R_t^{-1}M^\top)\) and
\(\log|R_t\otimes R_s|=N\log|R_t|+T\log|R_s|\) — so it scales to many regions and
times. The temporal range \(\nu\) is estimated alongside \(\beta,\sigma^2,\phi\).

## The data

One row per region **and** time, with the geometry repeated. Here four years on a
6×6 lattice:

```{r, purl = FALSE}
library(SDALGCP2)
library(sf)

set.seed(7)
shp <- st_sf(geometry = st_make_grid(
  st_as_sfc(st_bbox(c(xmin = 0, ymin = 0, xmax = 18, ymax = 18))), n = c(6, 6)))
N <- nrow(shp); times <- 2019:2022; T <- length(times)

# simulate separable space-time counts (true beta=0.5)
pts <- sda_points(shp, delta = 1.4, method = 3)
Rs  <- precompute_corr(pts, 3)$R[, , 1]
Rt  <- SDALGCP2:::.temporal_corr(seq_len(T), 1.5, 0.5)
L   <- t(chol(0.4 * kronecker(Rt, Rs)))
x1  <- rnorm(N * T); pop <- round(runif(N * T, 1000, 5000))
y   <- rpois(N * T, pop * exp(-6 + 0.5 * x1 + as.numeric(L %*% rnorm(N * T))))

dat <- st_sf(
  data.frame(cases = y, x1 = x1, pop = pop, year = rep(times, each = N)),
  geometry = st_geometry(shp)[rep(seq_len(N), T)])
```

![](t3_data.png)

## Fit

```{r, purl = FALSE}
fit <- sdalgcp(cases ~ x1 + offset(log(pop)), data = dat, time = "year",
               control = sdalgcp_control(reanchor = 2))
summary(fit)
```

```
#> Coefficients:
#>             Estimate Std.Err z value Pr(>|z|)
#> (Intercept)  -6.167   0.107  -57.61   <2e-16 ***
#> x1            0.592   0.038   15.73   <2e-16 ***
#> sigma^2       0.491   0.205    2.39    0.017 *
#> nu            0.853   0.360    2.37    0.018 *
#>
#> Spatial scale phi: 2.02
```

`nu` is the temporal range (correlation length in years); `phi` the spatial range.

## Predict and map — pick a year and a quantity

`predict()` returns a long `sf` — one row per region **and** time, with columns
`relative_risk`, `adjusted_rr` and their standard errors (the same column names as
the spatial predictor, so you can `st_write()` or map it directly). The `plot()`
method selects a **time** and a **quantity** (`"relative_risk"`, `"adjusted_rr"`,
`"relative_risk_se"`, `"adjusted_rr_se"`, `"exceedance"`):

```{r, purl = FALSE}
pr <- predict(fit)

plot(pr, time = NULL, what = "adjusted_rr")          # facet all years
plot(pr, time = 2021, what = "relative_risk")        # relative risk, 2021 only
plot(pr, time = 2021, what = "exceedance", threshold = 1.3, which = "adjusted_rr")
```

Covariate-adjusted relative risk, all years:

![](t3_arr_facet.png)

| Relative risk, 2021 | P(adjusted RR > 1.3), 2021 |
|:---:|:---:|
| ![](t3_rr_2021.png) | ![](t3_exc_2021.png) |

You can equally call `plot(fit, time = 2021, what = "relative_risk")` directly on
the fit, and `pr` (the returned `sf`) is a long region × time table of every
quantity for further analysis or mapping.

## Covariates and confounding

The space–time model supports the same covariate extensions as the spatial one.
The covariate surface is taken to be **time-invariant** (a covariate that genuinely
changes over time can always go in as an ordinary `data` column); the intensity-
scale tilting is then computed once per region and reused across the times, fitted
by a Gauss–Newton loop around the space–time likelihood.

**Raster (intensity-scale) covariates** — a spatially continuous covariate supplied
as a `terra` raster, aggregated correctly under the log link (see the
[raster tutorial](raster-covariates.html)):

```{r, purl = FALSE}
library(terra)
fit_r <- sdalgcp(cases ~ elevation + offset(log(pop)), data = dat, time = "year",
                 rasters = elevation_raster)
```

**Misaligned covariates** — a covariate observed on a different (time-invariant)
support, e.g. pollution monitors, kriged to the candidate points with a Berkson
correction (see [misaligned covariates](misaligned-covariates.html)):

```{r, purl = FALSE}
fit_m <- sdalgcp(cases ~ pm25 + offset(log(pop)), data = dat, time = "year",
                 covariates = list(pm25 = monitors_sf))
```

**Spatial confounding** — when a covariate is spatially smooth it can be collinear
with the random effect; restricted spatial regression de-confounds the coefficient
(see [spatial confounding](spatial-confounding.html)). It generalises to space–time
by constraining the random effect orthogonal to the design over all region–times:

```{r, purl = FALSE}
fit_c <- sdalgcp(cases ~ x1 + offset(log(pop)), data = dat, time = "year",
                 control = sdalgcp_control(confounding = "restricted"))
```

The restricted space–time fit reduces exactly to the spatial restricted fit when
there is a single time, and forms the full \((N\!\cdot\!T)\times(N\!\cdot\!T)\)
covariance, so it is best suited to modest \(N\!\cdot\!T\).

## Tips

Spatio-temporal fits profile the spatial range `phi` on a grid (set it with
`sdalgcp_control(phi = ...)`); the temporal range `nu` is estimated continuously.
Increase `reanchor` if the variance parameters look unstable.
```