UNIVERSITAS PADJADJARAN

Fakultas Matematika dan Ilmu Pengetahuan Alam — Statistika

Bayesian Spatial Statistics • Hands-on Workshop

Spatiotemporal Data Analysis
with R-INLA

From Spatial Patterns to Bayesian Spatiotemporal Models

3 hours Hands-on R + INLA West Java + Bandung case studies

Universitas Padjadjaran
Fakultas Matematika dan Ilmu Pengetahuan Alam
Statistika
Workshop Module — Spatiotemporal Data Analysis with R-INLA

Workshop orientation

This workshop is designed as a compact transition from spatial thinking to a practical Bayesian spatiotemporal model. The central idea is simple: observations collected over administrative regions and over time are usually not independent. Nearby regions may resemble each other, adjacent time points may be correlated, and unusual combinations of place and time may create additional interaction patterns.

Understand
spatial, temporal, and space-time dependence.

Map
Jawa Barat and Bandung using continuous and interval legends.

Fit & interpret
Bayesian models using R-INLA.

Learning outcomes

By the end of the three-hour session, participants should be able to:

distinguish spatial, temporal, and spatiotemporal sources of variation;
create and interpret continuous and interval choropleth maps for Jawa Barat and Bandung;
create a neighborhood graph from polygon data;
formulate a hierarchical Poisson model for area-level risk data;
simulate spatiotemporal count data with a known latent structure;
fit spatial, temporal, and spatiotemporal models with R-INLA;
interpret posterior summaries, relative risk, uncertainty, DIC/WAIC, and posterior exceedance probability;
communicate model results using maps and time profiles.

Three-hour flow

00:00–00:20

Spatial & temporal thinking. What dependence means and why ordinary regression is often insufficient.

00:20–00:45

Jawa Barat + Bandung mapping. Polygons, data joins, continuous legends, interval classes, neighbors, and adjacency graphs.

00:45–01:15

Bayesian hierarchy. Likelihood, latent effects, priors, and relative risk.

01:15–01:35

Simulation laboratory. Generate a known space-time signal and examine it visually.

01:35–02:15

R-INLA modeling. Baseline, spatial Besag, temporal RW1, and additive space-time model.

02:15–02:45

Posterior interpretation. RR maps, credible intervals, exceedance probability, and model comparison.

02:45–03:00

Mini challenge. Modify one assumption, refit, and defend the interpretation.

1. Why spatiotemporal models?

Suppose \(Y_{it}\) is a count observed in region \(i=1,\ldots,n\) and time \(t=1,\ldots,T\). A standard Poisson model assumes

\[ Y_{it}\mid \lambda_{it} \sim \text{Poisson}(\lambda_{it}), \]

with independent observations after conditioning on covariates. In regional data, that assumption is often too strong because:

spatial dependence: adjacent regions tend to share infrastructure, environment, climate, socioeconomic conditions, or disease transmission pathways;
temporal dependence: values in month \(t\) are usually related to month \(t-1\);
space-time interaction: a local shock can occur in a particular region during a particular time period.

A useful epidemiological or risk-mapping parameterization is

\[ Y_{it}\mid RR_{it},E_{it}\sim \text{Poisson}(E_{it}RR_{it}), \]

where \(E_{it}\) is an expected count under a reference risk and \(RR_{it}\) is the latent relative risk. If population-time is used instead, the same log-offset machinery estimates a rate rather than this standardized RR parameterization.

Interpretation rule. If \(RR_{it}=1\), observed risk is consistent with the reference level implied by \(E_{it}\). Values above 1 indicate elevated risk; values below 1 indicate lower risk. In Bayesian analysis, always interpret the posterior distribution, not only the posterior mean.

2. Software setup

Required packages

Run the installation chunk once. The INLA repository is not CRAN, so it is installed from the official R-INLA repository.

# General packages
install.packages(c(
  "sf", "spdep", "dplyr", "tidyr", "ggplot2",
  "viridis", "geodata", "terra", "knitr", "patchwork", "scales"
))

# R-INLA
install.packages(
  "INLA",
  repos = c(
    getOption("repos"),
    INLA = "https://inla.r-inla-download.org/R/stable"
  ),
  dep = TRUE
)

Load packages:

library(sf)
library(spdep)
library(dplyr)
library(tidyr)
library(ggplot2)
library(viridis)
library(geodata)
library(terra)
library(INLA)
library(patchwork)
library(scales)

Reproducibility note. This version uses the local West Java and Bandung shapefiles included in the project. No administrative boundary download is required. Keep the entire data/ folder beside the Rmd file so the workshop remains portable.

3. Local spatial data: West Java and Bandung

In the revised workshop module we explicitly use the local maps supplied in the project:

data/jabar/JABAR.shp for kabupaten/kota in West Java;
data/bandung/BANDUNG.shp for kecamatan in Bandung;
data/bandung/DataG.csv for the Bandung disease-mapping example;
data/bandung/Data.csv for the male/female aggregation example.

This makes the module fully grounded in the local spatial context. So now the workshop no longer says “West Java” in theory while quietly drawing some other map behind the curtain. The map and the math are finally on speaking terms.

Case A — Jawa Barat
27 kabupaten/kota; a province-scale spatiotemporal risk-mapping exercise. Students move from polygons and adjacency to simulated area-time counts, Besag + RW1 + interaction, posterior RR, and exceedance probability.

Case B — Kota Bandung
30 kecamatan; a finer-scale urban disease-mapping exercise using the uploaded diarrhea data. Students compare observed counts, expected counts, raw RR, continuous choropleths, and interval-based risk classes.

Workshop mapping principle. Every substantive risk map is shown in two forms: (1) a continuous legend for numerical detail and (2) an interval/class legend for easier interpretation and communication.

Project structure for the maps

Spatiotemporal_Workshop_RINLA/
├── Spatiotemporal_Workshop_RINLA.Rmd
├── figures/
│   ├── logo-unpad.png
│   ├── jabar_admin_map.png
│   ├── jabar_adjacency_map.png
│   ├── jabar_spatiotemporal_rr.png
│   ├── jabar_mapping_workflow.png
│   ├── jabar_rr_continuous_interval.png
│   ├── jabar_exceedance_continuous_interval.png
│   ├── bandung_admin_map.png
│   ├── bandung_rr_map.png
│   ├── bandung_rr_continuous_interval.png
│   └── bandung_adjacency_map.png
└── data/
    ├── jabar/
    │   ├── JABAR.shp
    │   ├── JABAR.shx
    │   ├── JABAR.dbf
    │   └── data.csv
    └── bandung/
        ├── BANDUNG.shp
        ├── BANDUNG.shx
        ├── BANDUNG.dbf
        ├── BANDUNG.prj
        ├── Data.csv
        └── DataG.csv

3.1 Read the West Java polygon data

library(sf)
library(dplyr)
library(ggplot2)
library(spdep)
library(INLA)

jabar <- st_read("data/jabar/JABAR.shp", quiet = TRUE)

# the supplied file does not store CRS information explicitly,
# so we assign geographic coordinates (longitude/latitude)
if (is.na(st_crs(jabar))) {
  st_crs(jabar) <- 4326
}

jabar <- jabar |>
  st_make_valid() |>
  mutate(
    area_id   = row_number(),
    area_name = KABKOT
  )

nrow(jabar)
head(jabar[, c("area_id", "area_name")])

Important. The workshop now uses the actual West Java polygons from your uploaded file. This is preferable to downloading external boundaries because the project becomes more stable and fully portable.

3.2 Administrative map of West Java

Figure 3.1. Administrative map of West Java at kabupaten/kota level, used as the main spatial support for the workshop.

The first map is not just decoration. It determines the spatial unit of analysis. Once the unit is fixed, the data table must contain one record per area (or one record per area-time combination for a spatiotemporal analysis).

Let:

\[ i=1,2,\ldots,n \]

index the district/city in West Java. In this file, the number of areas is

\[ n = 27. \]

For each area \(i\), we will later define:

observed count \(Y_i\) or \(Y_{it}\),
expected count \(E_i\) or \(E_{it}\),
relative risk \(RR_i\) or \(RR_{it}\),
centroid coordinates for labeling and some exploratory analysis.

3.3 West Java neighborhood structure

A CAR/Besag model requires a graph describing which regions are neighbors. We use queen contiguity: two polygons are neighbors if they share at least one boundary point.

nb_jabar <- poly2nb(jabar, queen = TRUE)
summary(card(nb_jabar))

# create the graph file for INLA
nb2INLA("data/jabar/jabar.graph", nb_jabar)
g_jabar <- inla.read.graph("data/jabar/jabar.graph")

Figure 3.2. West Java neighbor graph. Blue points are centroids; orange segments connect neighboring areas under queen contiguity.

For a neighbor set \(\partial i\), the intrinsic CAR/Besag prior is written conditionally as

\[ u_i \mid u_{-i},\tau_u \sim \mathcal N\left( \frac{1}{n_i}\sum_{j\in \partial i} u_j,\; \frac{1}{\tau_u n_i} \right), \]

where \(n_i\) is the number of neighbors of area \(i\). In matrix form,

\[ \mathbf u \sim \mathcal N\left(\mathbf 0, \; \tau_u^{-1} Q^{-}\right), \qquad Q = D-W, \]

with:

\(W=(w_{ij})\) the adjacency matrix,
\(w_{ij}=1\) if areas \(i\) and \(j\) are neighbors, and \(0\) otherwise,
\(D=\operatorname{diag}(n_1,\ldots,n_n)\).

The matrix \(Q\) is singular in the intrinsic formulation, which is why the Besag effect is usually fitted with a constraint such as

\[ \sum_{i=1}^n u_i = 0. \]

3.4 West Java disease-mapping quantities on the map

The most basic disease-mapping workflow goes from population to expected count, then to observed count, and finally to risk measures.

For area \(i\):

\[ SMR_i = \frac{Y_i}{E_i}, \]

where \(Y_i\) is the observed number of cases and \(E_i\) is the expected count under a reference rate.

If the reference rate is

\[ r = \frac{\sum_i Y_i}{\sum_i N_i}, \]

then the expected count is

\[ E_i = N_i r. \]

In practice we frequently use the Poisson model

\[ Y_i \mid RR_i \sim \text{Poisson}(E_i RR_i). \]

Thus

\[ \mathbb E(Y_i\mid RR_i)=E_i RR_i, \qquad \operatorname{Var}(Y_i\mid RR_i)=E_i RR_i. \]

Figure 3.3. Illustration of the mapping workflow on West Java: observed count \(Y\), raw \(SMR=Y/E\), and a smoothed risk map for teaching the idea of spatial borrowing of strength.

Example of manual calculation for one area

Suppose one West Java district has:

population at risk \(N_i=250{,}000\),
reference rate \(r=0.0032\),
observed cases \(Y_i=980\).

Then the expected count is

\[ E_i = N_i r = 250{,}000\times 0.0032 = 800. \]

The standardized morbidity ratio is

\[ SMR_i = \frac{980}{800}=1.225. \]

Interpretation: the observed number of cases is about 22.5% higher than expected under the reference rate. In log-risk form,

\[ \log(RR_i)\approx \log(1.225)=0.203. \]

This manual step is extremely useful in class because students can see that the map color is a visual translation of a numerical quantity, not a magical heatmap summoned by software.

3.5 West Java spatiotemporal map illustration

To demonstrate the spatiotemporal idea, we attach time to the same map. Let

\[ Y_{it} \mid RR_{it} \sim \text{Poisson}(E_{it}RR_{it}), \qquad i=1,\ldots,n,\quad t=1,\ldots,T. \]

A common decomposition is

\[ \log(RR_{it}) = \beta_0 + \mathbf x_{it}^{\top}\boldsymbol\beta + u_i + v_i + \gamma_t + \phi_t + \delta_{it}, \]

where:

\(u_i\) = structured spatial effect (Besag/ICAR),
\(v_i\) = unstructured spatial effect (IID),
\(\gamma_t\) = structured temporal effect (RW1/RW2/AR1),
\(\phi_t\) = unstructured temporal shock,
\(\delta_{it}\) = space-time interaction.

West Java spatiotemporal RR illustration

Figure 3.4. Spatiotemporal illustration on the West Java map across multiple periods. The purpose is pedagogical: to help students distinguish spatial patterns from temporal evolution.

A useful classroom discussion is: which changes across time appear to be global temporal shifts, and which appear to be localized interactions? That question naturally motivates the term \(\delta_{it}\).

Continuous versus interval risk map — West Java

For teaching, a continuous color bar is useful for seeing subtle numerical differences, while an interval legend is often easier for discussion with decision makers. The workshop therefore displays both.

West Java continuous and interval RR maps

Figure 3.5. West Java risk mapping shown with a continuous scale and with explicit risk intervals. The interval map is designed for workshop interpretation rather than hiding the thresholds inside a color gradient.

The default workshop intervals are:

RR interval	Label	Interpretation
\(RR<0.80\)	Lower	clearly below the reference level
\(0.80\le RR<1.00\)	Below reference	below 1 but relatively close to reference
\(1.00\le RR<1.20\)	Moderate	modest elevation above reference
\(1.20\le RR<1.50\)	High	substantively elevated risk
\(RR\ge1.50\)	Very high	strong elevation requiring closer attention

These thresholds are teaching thresholds, not universal scientific cut-points. In an applied study, the intervals should be chosen from substantive knowledge, surveillance standards, quantiles, or another defensible rule. The code below is a preview of the mapping grammar; the fully executable version appears after the spatiotemporal data are simulated in Section 5.

map_jabar_demo <- map_sim |>
  mutate(risk_class = classify_rr(RR_true))

p_jabar_cont <- ggplot(map_jabar_demo) +
  geom_sf(aes(fill = RR_true), color = "white", linewidth = 0.25) +
  scale_fill_viridis_c(
    option = "C",
    name = "Relative risk",
    labels = scales::label_number(accuracy = 0.01)
  ) +
  labs(
    title = "Continuous scale",
    subtitle = paste("Jawa Barat — period", month_to_show)
  ) +
  workshop_map_theme()

p_jabar_class <- ggplot(map_jabar_demo) +
  geom_sf(aes(fill = risk_class), color = "white", linewidth = 0.25) +
  scale_fill_manual(
    values = risk_palette,
    drop = FALSE,
    name = "Risk interval",
    guide = guide_legend(ncol = 1, byrow = TRUE)
  ) +
  labs(
    title = "Interval / class scale",
    subtitle = "Explicit thresholds support interpretation"
  ) +
  workshop_map_theme()

p_jabar_cont + p_jabar_class +
  patchwork::plot_annotation(
    title = "Case A — Jawa Barat: Two Complementary Risk Maps",
    subtitle = "Use the continuous map for detail; use the interval map for communication."
  )

3.6 Read the Bandung map and data

Bandung is included as a finer-scale urban example. It is excellent for showing what happens when the spatial unit becomes smaller and the patterns become more heterogeneous.

bandung <- st_read("data/bandung/BANDUNG.shp", quiet = TRUE) |>
  st_make_valid() |>
  mutate(area_id = row_number(), area_name = NAMA)

bandung_dat <- read.csv("data/bandung/DataG.csv")

bandung_dat <- bandung_dat |>
  mutate(name_join = trimws(toupper(Kecamatan)))

bandung <- bandung |>
  mutate(name_join = trimws(toupper(NAMA))) |>
  left_join(bandung_dat, by = "name_join")

Why Bandung? West Java is ideal for the province-level story, while Bandung is ideal for a city-level story. Together they let students see the difference between macro-scale spatial variation and micro-area urban heterogeneity.

3.7 Administrative map of Bandung

Figure 3.6. Administrative map of Bandung at kecamatan level.

For the Bandung example, let \(i=1,\ldots,30\) index the kecamatan. The uploaded data set DataG.csv contains the following key variables:

Variable	Meaning
`Penduduk`	population at risk
`Diare`	observed number of diarrhea cases
`E`	expected count
`RR`	estimated/raw relative risk
`Group`	grouping indicator used in the earlier script

3.8 Bandung mapping illustration: observed cases and relative risk

Figure 3.7. Bandung maps for observed cases and relative risk, built from the uploaded DataG.csv and BANDUNG.shp files.

Continuous versus interval risk map — Bandung

Figure 3.8. Bandung relative risk shown both as a continuous gradient and as interpretable risk intervals. The categorical legend makes the risk thresholds explicit.

bandung_map <- bandung |>
  mutate(risk_class = classify_rr(RR))

p_bdg_cont <- ggplot(bandung_map) +
  geom_sf(aes(fill = RR), color = "white", linewidth = 0.25) +
  scale_fill_viridis_c(
    option = "C",
    name = "Relative risk",
    labels = scales::label_number(accuracy = 0.01)
  ) +
  labs(
    title = "Continuous scale",
    subtitle = "Bandung diarrhea relative risk by kecamatan"
  ) +
  workshop_map_theme()

p_bdg_class <- ggplot(bandung_map) +
  geom_sf(aes(fill = risk_class), color = "white", linewidth = 0.25) +
  scale_fill_manual(
    values = risk_palette,
    drop = FALSE,
    name = "Risk interval",
    guide = guide_legend(ncol = 1, byrow = TRUE)
  ) +
  labs(
    title = "Interval / class scale",
    subtitle = "Same RR values, easier categorical interpretation"
  ) +
  workshop_map_theme()

p_bdg_cont + p_bdg_class +
  patchwork::plot_annotation(
    title = "Case B — Kota Bandung: Diarrhea Risk Mapping",
    subtitle = "Observed local data at kecamatan level"
  )

The interval table can also be summarized numerically:

bandung_map |>
  sf::st_drop_geometry() |>
  count(risk_class, .drop = FALSE, name = "Number of kecamatan")

These maps are ideal for explaining the difference between counts and risk.

For example, a kecamatan may have a large number of cases simply because it has a large population. That does not automatically imply unusually high risk. The risk measure adjusts for the expected count:

\[ RR_i = \frac{\mu_i}{E_i}, \]

where \(\mu_i=\mathbb E(Y_i\mid \text{model})\). Under a simple plug-in interpretation, the raw risk estimate is often approximated by

\[ \widehat{RR}_i \approx \frac{Y_i}{E_i}. \]

Example from the Bandung data

Suppose a kecamatan has:

observed diarrhea cases \(Y_i=460\),
expected cases \(E_i=293\).

Then

\[ SMR_i = \frac{460}{293} \approx 1.57. \]

This means the observed count is approximately 57% higher than expected. The map color for that kecamatan should therefore fall into a higher-risk class.

3.9 Bandung neighborhood map

nb_bandung <- poly2nb(bandung, queen = TRUE)
summary(card(nb_bandung))
nb2INLA("data/bandung/bandung.graph", nb_bandung)
g_bandung <- inla.read.graph("data/bandung/bandung.graph")

Figure 3.9. Bandung neighbor structure. This graph can be used directly for Besag, BYM, or BYM2 components in INLA.

The Bandung graph is especially useful when students want to compare a province-level map and a city-level map. The smaller spatial unit usually produces:

more local heterogeneity,
stronger visual contrast between neighboring units,
higher sensitivity to sparse counts,
greater need for smoothing.

3.10 Recommended mapping workflow for the workshop

Read and validate the shapefile

Join the attribute data to the polygons

Create administrative and choropleth maps

Build the neighborhood graph for INLA

A compact but complete mapping workflow in R is:

# 1. read polygon data
jabar   <- st_read("data/jabar/JABAR.shp", quiet = TRUE)
bandung <- st_read("data/bandung/BANDUNG.shp", quiet = TRUE)
if (is.na(st_crs(jabar))) st_crs(jabar) <- 4326

# 2. read attribute data
bdg <- read.csv("data/bandung/DataG.csv")

# 3. harmonize names and join
bandung <- bandung |>
  mutate(name_join = trimws(toupper(NAMA))) |>
  left_join(
    bdg |>
      mutate(name_join = trimws(toupper(Kecamatan))),
    by = "name_join"
  )

# 4. create exploratory maps
ggplot(jabar) + geom_sf(fill = "grey95", color = "grey40")

ggplot(bandung) +
  geom_sf(aes(fill = RR), color = "white", linewidth = 0.2) +
  scale_fill_viridis_c(option = "C") +
  labs(title = "Relative Risk across Bandung Kecamatan")

# 5. create neighborhood lists
nb_jabar   <- poly2nb(jabar, queen = TRUE)
nb_bandung <- poly2nb(bandung, queen = TRUE)

# 6. export graph objects for INLA
nb2INLA("data/jabar/jabar.graph", nb_jabar)
nb2INLA("data/bandung/bandung.graph", nb_bandung)

3.11 How the maps connect to the model

The maps are not just descriptive figures; they define the algebra of the model.

If the data are stored in long format, then each row corresponds to one \((i,t)\) combination:

\[ \bigl(i,t,Y_{it},E_{it},x_{it1},\ldots,x_{itp}\bigr). \]

The polygon map gives the geometry of \(i\), while the adjacency graph gives the structure used by the Besag or BYM prior. The temporal index \(t\) allows us to build RW1, RW2, or AR(1) effects.

Thus, the mapping pipeline and the modeling pipeline are connected as follows:

Mapping step	Modeling consequence
choose spatial units	defines the index \(i\)
build adjacency	defines \(W\) and \(Q=D-W\)
compute expected counts	produces \(E_{it}\) for the Poisson offset
map raw \(SMR\)	provides exploratory risk visualization
fit INLA model	produces smoothed posterior risk
map posterior mean and exceedance probability	produces inferential disease maps

Checkpoint 2. Ask students to compare West Java and Bandung: which map is more appropriate for province-level planning, and which is more appropriate for identifying local hotspots? Then ask how the answer changes when the target is policy, surveillance, or intervention design.

4. Bayesian hierarchical formulation

A practical additive spatiotemporal model is

\[ Y_{it}\mid RR_{it},E_{it}\sim \text{Poisson}(E_{it}RR_{it}), \]

\[ \log(RR_{it}) = \beta_0 + \beta_1 x_{it} + u_i + \gamma_t + \delta_{it}. \]

The components are:

Term	Meaning	Typical INLA model
\(\beta_0\)	global log-risk	fixed intercept
\(\beta_1x_{it}\)	observed covariate effect	fixed effect
\(u_i\)	structured spatial effect	Besag / CAR
\(\gamma_t\)	structured temporal trend	RW1 or RW2
\(\delta_{it}\)	residual space-time interaction	IID in the introductory model

For the structured spatial effect, a Besag prior informally says that \(u_i\) should be similar to the average of neighboring effects. A first-order temporal random walk assumes

\[ \gamma_t - \gamma_{t-1}\sim N(0,\tau_\gamma^{-1}). \]

Why INLA? Integrated Nested Laplace Approximation provides fast deterministic approximations to posterior marginals for latent Gaussian models. For many spatial and spatiotemporal models, this is substantially faster than generic MCMC while retaining a fully Bayesian interpretation.

5. Simulate a teachable spatiotemporal dataset

We now generate data where the truth is known. This is valuable because we can ask whether the model recovers the simulated structure.

Construct spatial and temporal effects

n_area <- nrow(jabar)
n_time <- 12

# Row-standardized neighborhood matrix
W <- spdep::nb2mat(nb_jabar, style = "W", zero.policy = TRUE)

# Smooth a random signal through the neighborhood structure
z <- rnorm(n_area)
spatial_raw <- as.numeric(solve(diag(n_area) - 0.55 * W, z))
spatial_eff <- as.numeric(scale(spatial_raw)) * 0.35

# Smooth temporal signal
rw <- cumsum(rnorm(n_time, 0, 0.18))
temporal_eff <- as.numeric(scale(rw)) * 0.28

sim <- tidyr::expand_grid(
  time_id = 1:n_time,
  area_id = 1:n_area
) |>
  arrange(time_id, area_id) |>
  mutate(
    st_id = row_number(),
    x = rnorm(n()),
    E = round(runif(n(), 25, 120)),
    interaction = rnorm(n(), 0, 0.14),
    eta_true = -0.05 +
      0.30 * x +
      spatial_eff[area_id] +
      temporal_eff[time_id] +
      interaction,
    RR_true = exp(eta_true),
    lambda = E * RR_true,
    y = rpois(n(), lambda)
  )

head(sim)

Here E is a simulated expected count under the reference risk. The latent linear predictor is therefore directly on the log-relative-risk scale, so RR_true = exp(eta_true) and a value of 1 has its usual risk-mapping interpretation.

Visualize one month

month_to_show <- 12

map_sim <- jabar |>
  left_join(
    sim |> filter(time_id == month_to_show),
    by = "area_id"
  )

ggplot(map_sim) +
  geom_sf(aes(fill = RR_true), color = "white", linewidth = 0.25) +
  scale_fill_viridis_c(option = "C", name = "True RR",
                       labels = scales::label_number(accuracy = 0.01)) +
  labs(
    title = paste("Simulated relative risk — period", month_to_show),
    subtitle = "Known latent signal before fitting the Bayesian model",
    caption = "Case A — Jawa Barat | continuous legend"
  ) +
  workshop_map_theme()

Workshop map: continuous and interval scales

Now that map_sim exists, the following chunk is executable and generates the two complementary Jawa Barat maps used during the hands-on session.

map_jabar_demo <- map_sim |>
  mutate(risk_class = classify_rr(RR_true))

p_jabar_cont <- ggplot(map_jabar_demo) +
  geom_sf(aes(fill = RR_true), color = "white", linewidth = 0.25) +
  scale_fill_viridis_c(
    option = "C",
    name = "Relative risk",
    labels = scales::label_number(accuracy = 0.01)
  ) +
  labs(
    title = "Continuous scale",
    subtitle = paste("Jawa Barat — period", month_to_show),
    caption = "Continuous legends preserve small numerical differences."
  ) +
  workshop_map_theme()

p_jabar_class <- ggplot(map_jabar_demo) +
  geom_sf(aes(fill = risk_class), color = "white", linewidth = 0.25) +
  scale_fill_manual(
    values = risk_palette,
    drop = FALSE,
    name = "Risk interval",
    guide = guide_legend(ncol = 1, byrow = TRUE)
  ) +
  labs(
    title = "Interval / class scale",
    subtitle = "Explicit RR thresholds",
    caption = "Class maps are easier to discuss in a workshop or policy briefing."
  ) +
  workshop_map_theme()

p_jabar_cont + p_jabar_class +
  patchwork::plot_annotation(
    title = "Case A — Jawa Barat: Simulated Spatiotemporal Risk",
    subtitle = "The same period shown with continuous and interval legends"
  )

Temporal profile

sim |>
  group_by(time_id) |>
  summarise(mean_RR = mean(RR_true), .groups = "drop") |>
  ggplot(aes(time_id, mean_RR)) +
  geom_line(linewidth = 1) +
  geom_point(size = 2) +
  labs(
    x = "Time",
    y = "Mean simulated RR",
    title = "Average risk profile over time"
  ) +
  theme_minimal(base_size = 12)

6. Fit models with R-INLA

We fit models of increasing complexity so that participants can see what each latent component contributes.

Model A — baseline Poisson regression

m_a <- INLA::inla(
  y ~ 1 + x,
  family = "poisson",
  data = sim,
  offset = log(E),
  control.predictor = list(compute = TRUE),
  control.compute = list(dic = TRUE, waic = TRUE, cpo = TRUE)
)

m_a$summary.fixed

Model B — add structured spatial effect

m_b <- INLA::inla(
  y ~ 1 + x +
    f(
      area_id,
      model = "besag",
      graph = g_jabar,
      scale.model = TRUE,
      constr = TRUE
    ),
  family = "poisson",
  data = sim,
  offset = log(E),
  control.predictor = list(compute = TRUE),
  control.compute = list(dic = TRUE, waic = TRUE, cpo = TRUE)
)

Model C — spatial + temporal + interaction

m_c <- INLA::inla(
  y ~ 1 + x +
    f(
      area_id,
      model = "besag",
      graph = g_jabar,
      scale.model = TRUE,
      constr = TRUE
    ) +
    f(
      time_id,
      model = "rw1",
      scale.model = TRUE,
      constr = TRUE
    ) +
    f(st_id, model = "iid"),
  family = "poisson",
  data = sim,
  offset = log(E),
  control.predictor = list(compute = TRUE),
  control.compute = list(
    dic = TRUE,
    waic = TRUE,
    cpo = TRUE,
    config = TRUE,
    return.marginals.predictor = TRUE
  )
)

m_c$summary.fixed

m_c$summary.hyperpar

Interpret the fixed effect. The coefficient for x is on the log-relative-risk scale. Therefore exp(beta_x) is a multiplicative relative-risk ratio associated with a one-unit increase in x, conditional on the latent effects.

7. Compare the models

comparison <- data.frame(
  Model = c("A: fixed only", "B: + spatial", "C: + spatial + temporal + interaction"),
  DIC = c(m_a$dic$dic, m_b$dic$dic, m_c$dic$dic),
  WAIC = c(m_a$waic$waic, m_b$waic$waic, m_c$waic$waic)
)

comparison

Lower DIC/WAIC is generally preferred, but the difference should be interpreted alongside scientific plausibility, residual diagnostics, and posterior uncertainty.

CPO diagnostic

cpo_problem_rate <- mean(m_c$cpo$failure > 0)
cpo_problem_rate

## [1] 0

A non-negligible failure rate can indicate numerical or predictive issues that deserve attention.

8. Posterior relative risk

For a Poisson model with expected-count offset \(\log(E_{it})\), the fitted mean is \(E_{it}\times RR_{it}\). Therefore, posterior fitted counts can be converted to relative risk by dividing by the expected count.

sim$post_rr_mean <- m_c$summary.fitted.values$mean / sim$E
sim$post_rr_lwr  <- m_c$summary.fitted.values$`0.025quant` / sim$E
sim$post_rr_upr  <- m_c$summary.fitted.values$`0.975quant` / sim$E

sim |>
  select(area_id, time_id, RR_true, post_rr_mean, post_rr_lwr, post_rr_upr) |>
  head()

Posterior map

map_post <- jabar |>
  left_join(
    sim |>
      filter(time_id == month_to_show) |>
      select(area_id, post_rr_mean, post_rr_lwr, post_rr_upr),
    by = "area_id"
  )

map_post <- map_post |>
  mutate(post_risk_class = classify_rr(post_rr_mean))

p_post_cont <- ggplot(map_post) +
  geom_sf(aes(fill = post_rr_mean), color = "white", linewidth = 0.25) +
  scale_fill_viridis_c(
    option = "C", name = "Posterior RR",
    labels = scales::label_number(accuracy = 0.01)
  ) +
  labs(
    title = "Posterior RR — continuous",
    subtitle = paste("Jawa Barat, period", month_to_show)
  ) +
  workshop_map_theme()

p_post_class <- ggplot(map_post) +
  geom_sf(aes(fill = post_risk_class), color = "white", linewidth = 0.25) +
  scale_fill_manual(
    values = risk_palette,
    drop = FALSE,
    name = "Posterior RR interval"
  ) +
  labs(
    title = "Posterior RR — intervals",
    subtitle = "Same posterior means grouped into risk classes"
  ) +
  workshop_map_theme()

p_post_cont + p_post_class +
  patchwork::plot_annotation(
    title = "Posterior Relative Risk Mapping",
    subtitle = "Spatial smoothing + temporal structure + interaction"
  )

9. Posterior exceedance probability

A useful decision-oriented quantity is

\[ P(RR_{it} > 1 \mid \text{data}). \]

Because the model is fitted with offset = log(E), the INLA linear predictor is

\[ \eta_{it}=\log(E_{it})+\log(RR_{it}). \]

Therefore,

\[ RR_{it}>1 \quad\Longleftrightarrow\quad \eta_{it}>\log(E_{it}). \]

The safest calculation is therefore to integrate the linear-predictor marginal above the threshold \(\log(E_{it})\). This also avoids the seq.default(): 'from' must be a finite number failure that can occur when a fitted-value marginal is absent or numerically non-finite.

# Robust posterior exceedance probability:
# P(RR_it > 1 | data) = P(eta_it > log(E_it) | data)
# where eta_it is the linear predictor including the offset log(E_it).

safe_pmarginal <- function(q, marginal) {
  if (is.null(marginal) || !is.matrix(marginal) || nrow(marginal) < 2) {
    return(NA_real_)
  }

  ok <- is.finite(marginal[, 1]) & is.finite(marginal[, 2])
  marginal <- marginal[ok, , drop = FALSE]

  if (nrow(marginal) < 2 || any(diff(marginal[, 1]) <= 0)) {
    return(NA_real_)
  }

  ans <- tryCatch(
    INLA::inla.pmarginal(q = q, marginal = marginal),
    error = function(e) NA_real_
  )

  if (!is.finite(ans)) NA_real_ else ans
}

# Default fallback: Gaussian approximation based on INLA posterior
# summary for the linear predictor. This is stable for knitting.
lp_sum <- m_c$summary.linear.predictor

if (is.null(lp_sum) || nrow(lp_sum) != nrow(sim)) {
  stop("INLA did not return summary.linear.predictor. Refit m_c with control.predictor = list(compute = TRUE).")
}

log_rr_mean <- lp_sum$mean - log(sim$E)
log_rr_sd   <- lp_sum$sd

sim$prob_rr_gt_1 <- ifelse(
  is.finite(log_rr_mean) & is.finite(log_rr_sd) & log_rr_sd > 0,
  stats::pnorm(log_rr_mean / log_rr_sd),
  NA_real_
)

# If full linear-predictor marginals are available, replace the Gaussian
# approximation with direct numerical integration of the INLA marginal.
lp_marginals <- m_c$marginals.linear.predictor

if (!is.null(lp_marginals) && length(lp_marginals) == nrow(sim)) {
  prob_exact <- vapply(
    seq_len(nrow(sim)),
    function(k) {
      p_below <- safe_pmarginal(
        q = log(sim$E[k]),
        marginal = lp_marginals[[k]]
      )
      if (is.finite(p_below)) 1 - p_below else NA_real_
    },
    numeric(1)
  )

  use_exact <- is.finite(prob_exact)
  sim$prob_rr_gt_1[use_exact] <- prob_exact[use_exact]
}

# Final protection against tiny numerical excursions outside [0,1]
sim$prob_rr_gt_1 <- pmin(1, pmax(0, sim$prob_rr_gt_1))

summary(sim$prob_rr_gt_1)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.000000 0.000275 0.371185 0.476365 0.999775 1.000000

sum(is.na(sim$prob_rr_gt_1))

## [1] 0

Why this fixes the knit error. The earlier code called inla.pmarginal() directly on fitted-value marginals and could fail when a marginal had a non-finite support. The revised version first computes a stable Gaussian approximation from summary.linear.predictor; when valid full linear-predictor marginals are available, it replaces that approximation with direct INLA marginal integration. Invalid marginals are skipped instead of terminating the knit.

Map high-probability areas:

map_prob <- jabar |>
  left_join(
    sim |>
      filter(time_id == month_to_show) |>
      select(area_id, prob_rr_gt_1),
    by = "area_id"
  )

map_prob <- map_prob |>
  mutate(prob_class = classify_prob(prob_rr_gt_1))

p_prob_cont <- ggplot(map_prob) +
  geom_sf(aes(fill = prob_rr_gt_1), color = "white", linewidth = 0.25) +
  scale_fill_viridis_c(
    option = "B",
    limits = c(0, 1),
    name = "P(RR > 1)",
    labels = scales::label_percent(accuracy = 1)
  ) +
  labs(
    title = "Continuous posterior probability",
    subtitle = paste("Jawa Barat, period", month_to_show)
  ) +
  workshop_map_theme()

p_prob_class <- ggplot(map_prob) +
  geom_sf(aes(fill = prob_class), color = "white", linewidth = 0.25) +
  scale_fill_manual(
    values = prob_palette,
    drop = FALSE,
    name = "Evidence interval"
  ) +
  labs(
    title = "Evidence classes",
    subtitle = "Intervals make posterior evidence easier to communicate"
  ) +
  workshop_map_theme()

p_prob_cont + p_prob_class +
  patchwork::plot_annotation(
    title = "Posterior Exceedance Probability",
    subtitle = "Continuous probability and explicit evidence intervals"
  )

Jawa Barat continuous and interval exceedance probability maps

Workshop preview. Exceedance probability can also be presented as a continuous probability or grouped into evidence intervals such as weak, moderate, strong, and very strong evidence.

Do not label an area “high risk” only because the posterior mean exceeds 1. A more defensible statement combines the posterior mean, credible interval, and exceedance probability. For example: “The posterior mean RR is 1.18 and the posterior probability that RR exceeds 1 is 0.94.”

10. Optional advanced model: structured space-time interaction

The introductory interaction above is IID. A more structured formulation can let the spatial field evolve through time by using an INLA group model. One practical extension is:

m_adv <- INLA::inla(
  y ~ 1 + x +
    f(
      area_id,
      model = "besag",
      graph = g_jabar,
      group = time_id,
      control.group = list(model = "rw1"),
      scale.model = TRUE,
      constr = TRUE
    ),
  family = "poisson",
  data = sim,
  offset = log(E),
  control.predictor = list(compute = TRUE),
  control.compute = list(dic = TRUE, waic = TRUE, cpo = TRUE)
)

This should be introduced after participants understand the additive model because its interpretation and computational structure are more demanding.

11. Interpretation checklist

When reporting results, answer these questions in order:

What is the observational unit? District/city × time.
What is the likelihood? Poisson for counts, with an expected-count offset.
What is the target parameter? Relative risk, \(RR_{it}\).
What is the fixed-effect interpretation? Multiplicative change in RR through \(\exp(\beta)\).
What latent dependence is modeled? Spatial Besag, temporal RW1, and interaction.
How uncertain are the estimates? Posterior SD and 95% credible intervals.
Is elevated risk supported? Use \(P(RR>1\mid y)\), not only a point estimate.
Does the model improve fit/prediction? Compare WAIC/DIC/CPO and inspect residual behavior.
Is the conclusion causal? No, unless the study design and identification strategy justify causality.

12. Mini challenge — final 15 minutes

Team challenge. Choose one modification, rerun the model, and prepare a 60-second explanation.

Replace RW1 by RW2 and compare WAIC.
Remove the interaction and explain what changes in the posterior map.
Increase the simulated spatial signal from 0.35 to 0.60 and assess recovery.
Increase the time series from 12 to 24 periods.
Add a second covariate and explain its posterior multiplicative effect.

A good answer should include: model change → posterior change → substantive interpretation → uncertainty.

13. Common mistakes

Watch for these errors:

treating neighboring regions as independent;
interpreting a random effect as a directly observed risk factor;
ignoring or misinterpreting the expected-count offset;
comparing maps that use inconsistent color scales;
reporting only posterior means without credible intervals;
assuming a lower DIC/WAIC automatically makes a model scientifically correct;
calling a posterior association “causal” without a causal identification design.

14. Suggested workshop deliverable

Each participant submits a short reproducible report containing:

one map of the spatial units and neighborhood structure;
the final model formula;
a table of fixed effects and hyperparameters;
one posterior RR map;
one posterior exceedance-probability map;
a 150–200 word interpretation covering uncertainty and limitations.

15. Take-home summary

Spatial structure
Besag/CAR effects borrow information from neighboring regions and stabilize noisy area estimates.

Temporal structure
RW1/RW2 priors encode smooth evolution while allowing data-driven departures.

Spatiotemporal structure
Interaction terms capture localized departures that cannot be explained by additive spatial and temporal components alone.

Bayesian interpretation
Use posterior distributions, credible intervals, and exceedance probabilities—not point estimates alone.

Spatiotemporal Data Analysiswith R-INLA

Workshop orientation

Learning outcomes

Three-hour flow

1. Why spatiotemporal models?

2. Software setup

Required packages

3. Local spatial data: West Java and Bandung

Project structure for the maps

3.1 Read the West Java polygon data

3.2 Administrative map of West Java

3.3 West Java neighborhood structure

3.4 West Java disease-mapping quantities on the map

Example of manual calculation for one area

3.5 West Java spatiotemporal map illustration

Continuous versus interval risk map — West Java

3.6 Read the Bandung map and data

3.7 Administrative map of Bandung

3.8 Bandung mapping illustration: observed cases and relative risk

Continuous versus interval risk map — Bandung

Example from the Bandung data

3.9 Bandung neighborhood map

3.10 Recommended mapping workflow for the workshop

3.11 How the maps connect to the model

4. Bayesian hierarchical formulation

5. Simulate a teachable spatiotemporal dataset

Construct spatial and temporal effects

Visualize one month

Workshop map: continuous and interval scales

Temporal profile

6. Fit models with R-INLA

Model A — baseline Poisson regression

Model B — add structured spatial effect

Model C — spatial + temporal + interaction

7. Compare the models

CPO diagnostic

8. Posterior relative risk

Posterior map

9. Posterior exceedance probability

10. Optional advanced model: structured space-time interaction

11. Interpretation checklist

12. Mini challenge — final 15 minutes

13. Common mistakes

14. Suggested workshop deliverable

15. Take-home summary

Spatiotemporal Data Analysis
with R-INLA