Bayesian Spatiotemporal Analysis, Mapping and Forecasting

Four-Day International Workshop Course Book: Spatial Concepts, Bayesian Inference and R-INLA

INLA check: params$run_inla = TRUE and the INLA package is available. INLA examples will be attempted during render.

Workshop at a glance. Four days (16 hours), led by Henk Folmer and I Gede Nyoman Mindra Jaya. Day 1 develops spatial concepts, Bayesian reasoning and R-INLA. Days 2–4 introduce the three Jaya–Folmer papers through simulated dengue-style data over the real 30-kecamatan Kota Bandung map. Each day combines a concise explanation, mathematical model, mapped illustration, reproducible R code and exercises.

DAY 1	DAY 2 — PAPER 1	DAY 3 — PAPER 2	DAY 4 — PAPER 3
Spatial & Bayesian concepts; R-INLA	Pure spatiotemporal mapping	AHC clustering & varying coefficients	Prediction & confounding

Day 1 - Spatial Concepts, Bayesian Foundations and Introduction to R-INLA

Bandung illustration. Start from the 30 real kecamatan polygons. Simulate monthly dengue counts, define rook/queen neighbors, calculate expected cases and SIR, then inspect temporal and spatial dependence. Interactive widgets below illustrate prior learning and BYM2 mixing.

PPTX alignment: spatial concepts (slides 4–46). The spatial lecture proceeds through spatial econometrics and its applications, Tobler’s law, spatial dependence versus heterogeneity, contiguity and distance weights, temporal and spatial lag operators, and global spatial autocorrelation. The short formulas and Bandung maps below are the practical counterparts of these slides.

A compact spatial reference. Spatial dependence describes statistical association among geographically linked observations; spatial heterogeneity describes variation in the distribution or relationships from place to place. Neither concept implies a causal spillover. An econometric model with \(Wy\) on the right-hand side makes a stronger structural assumption than a diagnostic Moran’s I calculation.

PPTX concept	Core formulation	Bandung interpretation
Rook / queen adjacency	\(a_{ij}=1\) for a shared edge / edge or vertex	Which kecamatan can share information?
Power and exponential weights	\(w_{ij}=d_{ij}^{-p}\) or \(e^{-d_{ij}/h}\), \(i\ne j\)	Nearby centroids receive larger weights
Spatial lag	\((Wy)_i=\sum_jw_{ij}y_j\)	Neighboring dengue risk average
Temporal lag	\(Ly_{it}=y_{i,t-1}\), \(L^ky_{it}=y_{i,t-k}\)	Risk in the previous month or season
Higher-order spatial paths	\(W^ky\)	Influence propagated across \(k\) graph steps
Moran’s I	\(I=(n/S_0)\,\frac{(y-\bar y)^TW(y-\bar y)}{\sum_i(y_i-\bar y)^2}\)	Global evidence of clustered or dispersed risk

Use a symmetric binary graph for ICAR/BYM2 model construction. A row-standardized \(W\) is useful for spatial lag interpretation but is not generally symmetric. Statistical significance for Moran’s I requires a specified null distribution (for example, permutation); \(I>0\) alone is descriptive.

Core sources for this session: spatial dependence and spatiotemporal thinking are introduced using Tobler (1970), Waller and Gotway (2004), and Cressie and Wikle (2011). The official course structure follows the uploaded Course Outline (n.d.).

Workshop motivation.

The purpose of this course is to familiarize participants with the theory and practice of estimation, testing, and forecasting of spatiotemporal panel models. The emphasis is practical but mathematically explicit: every model is introduced through its purpose, its formal structure, and a reproducible Bandung City example.

By the end of Day 1, participants should understand the language used throughout the course. Before fitting any Bayesian model, the analyst must know what is meant by space, time, spatial dependence, temporal dependence, and spatiotemporal dependence.

Inline citation policy for this course book. Claims about disease transmission, spatial dependence, Bayesian priors, INLA, model-selection criteria, forecasting, clustering, and public-health interpretation are cited in the text at the point where the claim is made. The reference list at the end is only a companion list; it is not a substitute for inline citation.

The notions of space, spatial dependence, temporal dependence, and spatiotemporal dependence

Spatial dependence is introduced here using the general principle that nearby places often contain related information, while spatiotemporal dependence extends this logic across both space and time (Tobler, 1970; Cressie and Wikle, 2011).

Spatiotemporal modeling combines two simple questions:

1. Where? Does disease risk differ from one location to another?
2. When? Does disease risk change from one time period to another?

A standard table may only show numbers, but a disease-surveillance table has hidden geography and hidden time structure. If the same district appears repeatedly over many months, and neighboring districts influence one another, then ordinary independent analysis is usually too simple.

Think of a disease map like a weather map. Rain in one district may be related to rain in a nearby district, and rain today may be related to rain yesterday. Infectious disease risk often behaves similarly: it can cluster geographically and persist temporally.

In this workshop, space is represented by administrative areas. For the case study, the areas are Bandung City districts (kecamatan). A spatial unit \(i\) is therefore not only a name in a table; it is a polygon with a location, area, boundary, centroid, and neighbors.

Spatial dependence means that nearby areas tend to have related risks:

\[ \mathrm{Corr}(Y_{it},Y_{jt}) \neq 0 \quad \text{when } i \sim j \]

where \(i \sim j\) means that areas \(i\) and \(j\) are neighbors. In infectious disease mapping, spatial dependence can arise from shared climate, mobility, vector habitat, sanitation, socioeconomic conditions, reporting systems, or cross-border transmission.

Temporal dependence means that the risk at time \(t\) is related to previous or future periods:

\[ \mathrm{Corr}(Y_{it},Y_{i,t-1}) \neq 0 \]

For dengue, temporal dependence is often linked to rainfall, humidity, mosquito life cycles, reporting delays, and seasonal behavior.

Spatiotemporal dependence means that the spatial pattern itself changes over time, or the temporal pattern differs by area:

\[ \mathrm{Corr}(Y_{it},Y_{j,t-1}) \neq 0 \]

This is the most important setting for surveillance. A district may become high risk after its neighbors were high risk in previous months, or a seasonal increase may be strong in one part of the city but weak elsewhere.

Concept	Statistical_question	Disease_mapping_interpretation
Space	Where is each observation located?	District polygons, centroids, boundaries, adjacency, and distance
Spatial dependence	Are neighboring areas similar at the same time?	Nearby districts may share vector habitat, mobility, or reporting conditions
Temporal dependence	Is the current period related to past periods?	Risk may persist because of seasonality, transmission chains, or delayed effects
Spatiotemporal dependence	Does spatial risk evolve over time?	Hotspots can move, expand, disappear, or interact with neighboring districts

Bayesian statistics vs frequentist statistics.

The comparison between Bayesian and frequentist statistics is presented as a practical modeling distinction, especially around how uncertainty is represented and updated (Gelman et al., 2013).

Both frequentist and Bayesian approaches are useful. The difference is how they treat unknown quantities and uncertainty.

Topic	Frequentist_view	Bayesian_view
Unknown parameters	Fixed but unknown constants	Random quantities described by probability distributions
Probability	Long-run frequency under repeated sampling	Degree of uncertainty conditional on available information
Prior information	Usually not part of the model	Entered explicitly through priors and hyperpriors
Uncertainty interval	Confidence interval from repeated-sampling logic	Credible interval from the posterior distribution
Disease map	Often a point estimate plus standard error	Posterior mean, credible interval, and exceedance probability for each area-time unit
Forecasting	Prediction interval from fitted sampling model	Posterior predictive distribution integrating parameter uncertainty

For spatiotemporal disease mapping, the Bayesian approach is attractive because the object of interest is not one coefficient only. We need a full map of latent risks \(\theta_{it}\), spatial effects \(u_i\), temporal effects \(\gamma_t\), interactions \(\delta_{it}\), and uncertainty for every district-month.

First disease-mapping quantities: cases, expected cases, SIR, and relative risk

The first mistake in disease mapping is to compare districts only by raw case counts. A district with many people naturally tends to have more cases. Therefore, a fair comparison needs a baseline. The baseline asks: how many cases would we expect if this district had the same average disease rate as the entire study region?

That baseline is called the expected number of cases. After we have expected cases, we compare observed cases to expected cases. This comparison is the standardized incidence ratio or SIR.

A useful way to remember the idea is:

\[ \mathrm{SIR}=\frac{\text{what actually happened}}{\text{what would be expected under average risk}}. \]

In disease mapping, raw case counts are not sufficient for comparing areas. Areas with larger populations usually have more cases, even when their underlying risk is not higher. For this reason, we use expected cases and relative risk.

The central question is not “which district has the largest number of cases?” but rather:

Given the population at risk, which district and time period has more disease than expected?

This distinction matters for decision-making. A large metropolitan district may report many cases simply because many people live there. A smaller district may report fewer cases but still have unusually high risk after adjusting for its population. Relative risk makes these comparisons fairer.

For an early warning system, Day 1 builds the foundation:

• The surveillance data must be organized as an area-time panel.
• Each observation must have cases and population at risk.
• Expected cases provide the baseline.
• SIR gives the first crude risk signal.
• Maps and time-series plots reveal whether risk varies across space, time, or both.

Workshop scenario. Imagine a provincial disease-control center monitoring monthly dengue-like fever across 30 urban districts. Decision-makers want to know where vector control teams should be deployed next month. Raw case counts are politically visible, but relative risk is scientifically more appropriate because it accounts for different population sizes.

Let:

• \(Y_{it}\): number of cases in area \(i = 1,\ldots,n\) and period \(t = 1,\ldots,T\).
• \(N_{it}\): population at risk.
• \(E_{it}\): expected number of cases.
• \(\theta_{it}\): relative risk.

The basic formulas are:

\[ p = \frac{\sum_{i=1}^{n}\sum_{t=1}^{T}Y_{it}}{\sum_{i=1}^{n}\sum_{t=1}^{T}N_{it}}, \qquad E_{it} = N_{it}p \]

\[ \widehat{\theta}_{it}^{ML} = \mathrm{SIR}_{it} = \frac{Y_{it}}{E_{it}} \]

Interpretation:

• \(\mathrm{SIR}_{it} = 1\): risk is equal to the study-region average.
• \(\mathrm{SIR}_{it} > 1\): observed cases are higher than expected.
• \(\mathrm{SIR}_{it} < 1\): observed cases are lower than expected.

SIR is an unbiased estimator of relative risk, but it can be highly unstable when the number of cases or the population at risk is small. The thesis therefore uses Bayesian spatiotemporal smoothing so that information from spatial and temporal neighbors contributes to risk estimation.

Mathematical formulation.

How to read the first model formula.

Beginners often find the Poisson formula intimidating. Read it piece by piece:

\[ Y_{it}\mid \theta_{it} \sim \mathrm{Poisson}(E_{it}\theta_{it}). \]

This means:

Symbol	Plain-language meaning
\(Y_{it}\)	observed number of cases in area \(i\) at time \(t\)
\(E_{it}\)	expected number of cases in area \(i\) at time \(t\)
\(\theta_{it}\)	relative risk in area \(i\) at time \(t\)
\(E_{it}\theta_{it}\)	model-based mean number of cases
Poisson	probability distribution for count data

So the model says: the observed cases are random counts, and their average level is determined by expected cases multiplied by relative risk.

The most basic model is:

\[ Y_{it}\mid \theta_{it} \sim \mathrm{Poisson}(E_{it}\theta_{it}) \]

with the log link:

\[ \log(\mu_{it}) = \log(E_{it}) + \log(\theta_{it}) \]

If the data are overdispersed or contain many zeros, the negative binomial model can be used:

\[ Y_{it}\mid \theta_{it}, \kappa \sim \mathrm{NegBin}(E_{it}\theta_{it}, \kappa) \]

with:

\[ \mathbb{E}(Y_{it}) = E_{it}\theta_{it}, \qquad \mathrm{Var}(Y_{it}) = E_{it}\theta_{it} + \frac{(E_{it}\theta_{it})^2}{\kappa} \]

Important assumptions at this stage:

• Cases are non-negative counts.
• The population at risk is measured for the same area and time period as the cases.
• The reference rate \(p\) represents the average disease rate over the study region and study period.
• The offset \(\log(E_{it})\) must enter the model with coefficient fixed at 1.
• SIR is descriptive; it is not yet a smoothed model-based estimate.

The first diagnostic question is whether SIR appears noisy. If many small-population areas jump from very low to very high SIR across adjacent months, the analyst should be cautious and move toward Bayesian smoothing in Day 2.

Applied example: Bandung City simulation

Beginner workflow for the simulation.

The simulation is not random decoration. It is a controlled learning environment. Participants can learn the full workflow without confidential health data.

The workflow is:

1. Create or load district polygons.
2. Create monthly population and case data.
3. Calculate the reference disease rate.
4. Calculate expected cases.
5. Calculate SIR.
6. Build neighbor relationships between districts.
7. Visualize maps and time trends.
8. Use these outputs as preparation for Bayesian modeling.

Why use simulated disease counts? Real surveillance data may contain privacy restrictions, reporting delays, or political sensitivity. Simulated counts allow every participant to run the same code, compare results, and focus on learning the method.

We now simulate all workshop examples on the actual administrative shape of Bandung City, not on an artificial square grid. The boundary used here is the GADM 4.1 level-3 administrative boundary for Indonesia, filtered to the 30 districts (kecamatan) of Kota Bandung. The disease counts remain simulated, but the spatial topology, district names, relative location, area, and adjacency structure follow Bandung’s real administrative map.

This is the teaching compromise used throughout the e-book:

• Real geography: Bandung district polygons, centroids, areas, and neighbors.
• Simulated surveillance data: monthly dengue-like cases, population at risk, weather-like covariates, and mobility proxy.
• Reproducible code: no confidential patient-level data are needed.
• Policy realism: outputs are interpreted as if a health office were planning district-level surveillance and vector-control actions in Bandung.

The simulation is intentionally designed to mimic realistic Bandung surveillance complexity:

• Bandung has heterogeneous district areas and population density.
• Northern districts are given higher elevation-like and humidity-like profiles.
• Rainy-season months increase dengue-like risk through precipitation and humidity.
• Spatially neighboring districts share some latent risk.
• Local pockets of high risk can appear even when the city-wide trend is moderate.

Important file-handling note. This edition ships with the actual Kota Bandung kecamatan boundary in data/bandung_kecamatan_gadm41.geojson (30 polygons) supplied by the course organizer. The Rmd loads that local file first. It never downloads a boundary. By default require_real_bandung_map: true makes rendering stop if the actual geometry is missing, so a schematic rectangular fallback can never be mistaken for Kota Bandung. Only the epidemiological counts, populations, weather variables and relative risks are simulated. The administrative boundary source should be identified in applied reports (GADM, 2022).

This first output is a data-structure check. Each row is one Bandung district-month. The most important columns are:

• district: Bandung district name.
• date and month_no: monthly time index.
• population: simulated population at risk.
• cases: simulated monthly dengue-like cases.
• expected: expected cases under the city-wide reference rate.
• sir: crude standardized incidence ratio.
• precip, temp, humidity, solar, mobility: simulated risk-factor proxies.
• rr_true: the unobserved risk surface used only for simulation validation.

Bandung district boundary overview.

district_id	district	area_km2	baseline_population	neighbors
1	Andir	4.09	57074	5
2	Antapani	4.11	91327	4
3	Arcamanik	7.47	105530	5
4	Astanaanyar	2.70	48719	6
5	Babakan Ciparay	7.05	81417	3
6	Bandung Kidul	5.38	67253	5
7	Bandung Kulon	6.76	79825	2
8	Bandung Wetan	3.33	95098	5
9	Batununggal	4.80	99209	2
10	Bojongloa Kaler	2.99	49055	3
11	Bojongloa Kidul	5.10	70021	4
12	Buahbatu	7.31	92921	5
13	Cibeunying Kaler	4.64	84914	3
14	Cibeunying Kidul	3.94	102055	3
15	Cibiru	6.85	93220	2
16	Cicendo	7.65	96806	4
17	Cidadap	8.22	99782	2
18	Cinambo	4.25	69286	4
19	Coblong	7.26	94443	5
20	Gedebage	9.68	126299	3
21	Kiaracondong	5.68	98213	4
22	Lengkong	5.81	95039	4
23	Mandalajati	4.74	82564	4
24	Panyileukan	5.25	84648	3
25	Rancasari	6.88	87147	2
26	Regol	4.61	72704	5
27	Sukajadi	5.10	69480	4
28	Sukasari	6.26	80487	2
29	Sumur Bandung	3.41	91667	5
30	Ujung Berung	6.02	83798	4

This table confirms the spatial units used by the simulation. The supplied GeoJSON provides the actual 30 Kota Bandung kecamatan boundaries. If that file is missing, the default configuration stops rendering. area_km2 is computed from projected polygons, baseline_population is simulated, and neighbors is derived from queen contiguity. Neighbor counts influence which areas share information in spatial smoothing.

ggplot(districts) +
  geom_sf(aes(fill = area_km2), color = "white", linewidth = 0.35) +
  geom_sf_text(aes(label = district_id), size = 3, color = "#21313c") +
  scale_fill_gradient(low = "#edf7f5", high = "#457b9d", name = "Area (km2)") +
  labs(
    title = "Bandung City District Boundary",
    subtitle = paste0(boundary_source, "; labels show district IDs")
  )

This map is the spatial frame for the entire workshop. Every later SIR map, model-based relative-risk map, cluster map, forecast map, and high-resolution map is drawn inside this boundary. If the subtitle says Teaching fallback, the polygons are simplified teaching polygons, not official administrative boundaries. The numeric labels link the map to the district inventory table.

Data summary.

What beginners should check in every disease-mapping dataset.

Before fitting any model, check the dataset like a public-health analyst, not only like a programmer. The model can only be as good as the data structure.

A minimum area-time dataset should contain:

Required column	Why it matters
Area identifier	Needed to group observations by district or region
Time identifier	Needed to order observations and build trends
Cases	Response variable for count models
Population at risk	Needed to calculate expected cases
Expected cases	Needed as offset or baseline
Geometry or area code	Needed to connect data to maps
Neighbor information	Needed for spatial smoothing models
Covariates, if used	Needed for explanatory regression models

districts	months	total_cases	total_population_time	reference_rate	mean_sir	median_sir	max_sir
30	36	18789	101439061	0.0002	1.0072	0.728	5.7587

This summary gives the scale of the simulated Bandung surveillance system. total_cases is the total simulated burden over all district-months. reference_rate is the city-wide baseline rate used to compute expected cases. mean_sir is usually close to 1 by construction, but median_sir and max_sir reveal skewness: outbreaks create a small number of district-months with much higher risk than the city average.

Minimum data validation checklist.

Before fitting any model, participants should verify the following conditions:

Check	Rationale	Result
No negative cases	Case counts must be non-negative integers	TRUE
Positive population	Population at risk must be strictly positive	TRUE
Positive expected cases	The Poisson offset must be positive	TRUE
Complete area-time panel	Every district should appear in every month	TRUE
Valid SIR values	SIR should be finite after expected cases are computed	TRUE

These checks look simple, but they prevent common workshop and field-data errors: duplicated months, missing population denominators, mismatched area codes, and invalid offsets. In real applications, most modeling failures happen before the model is fitted.

Monthly SIR maps.

selected_months <- c(1, 6, 12, 24, 30, 36)

sir_map <- districts |>
  select(district_id, district) |>
  left_join(
    panel |> filter(t %in% selected_months) |>
      mutate(month_label = format(date, "%Y-%m")),
    by = "district_id"
  )

ggplot(sir_map) +
  geom_sf(aes(fill = sir), color = "white", linewidth = 0.25) +
  facet_wrap(~month_label, ncol = 3) +
  scale_fill_gradient2(
    low = "#f1d9b1", mid = "#fffaf3", high = "#973d1c",
    midpoint = 1, limits = c(0, max(sir_map$sir, na.rm = TRUE)),
    name = "SIR"
  ) +
  labs(
    title = "Standardized Incidence Ratio (SIR) Map",
    subtitle = "Values above 1 indicate observed risk above the study-region average"
  )

Interpretation of this output:

• Each facet is one month in the Bandung simulation.
• Light-colored districts have SIR below the city average, while dark brown districts have higher SIR.
• The same district can move between low and high SIR because dengue-like risk is seasonal and stochastic.
• This is a crude map: it does not yet borrow information from neighboring districts or adjacent months.
• Participants should look for persistent spatial patterns and sudden local spikes. Those two patterns motivate the Bayesian pure model in Day 2.

Worked calculation: SIR, spatial lag and Moran’s I (Bandung, month 12). The code below computes the statistics from the simulated observations and the queen-contiguity matrix built from actual kecamatan polygons. It reports the arithmetic inputs and numerical result rather than printing the full data structure.

d12 <- panel[panel$t==12,][order(panel$district_id[panel$t==12]),]
y <- d12$sir
z <- y-mean(y)
S0 <- sum(W_row)
morans_I <- length(y)/S0 * as.numeric(t(z)%*%W_row%*%z)/sum(z^2)
lag_sir <- as.numeric(W_row%*%y)
worked <- data.frame(
  district=d12$district[1:5], cases=d12$cases[1:5],
  expected=round(d12$expected[1:5],2),
  SIR=round(y[1:5],3), neighbor_lag=round(lag_sir[1:5],3))
knitr::kable(worked,caption="Five Bandung districts: observed SIR and weighted neighbor SIR")

Five Bandung districts: observed SIR and weighted neighbor SIR

district	cases	expected	SIR	neighbor_lag
Andir	18	11.27	1.597	0.966
Antapani	14	18.03	0.776	0.725
Arcamanik	13	20.84	0.624	0.866
Astanaanyar	4	9.62	0.416	0.983
Babakan Ciparay	8	16.08	0.498	0.576

cat("Moran's I =",round(morans_I,4),"; n =",length(y),
    "; sum of weights =",round(S0,1))

## Moran's I = 0.2412 ; n = 30 ; sum of weights = 30

The manual numerator is \(z^TWz\) and the denominator is \(z^Tz\). Here \(I=(n/S_0)(z^TWz)/(z^Tz)\). The spatial lag is \((Wy)_i=\sum_jw_{ij}y_j\). For example, with three equally weighted neighbors whose SIR values are \(0.8,1.2,1.6\), the spatial lag is \((0.8+1.2+1.6)/3=1.2\).

One-sided randomization check for positive spatial association

Morans_I	upper_tail_permutation_p	permutations
0.2412	0.03	199

Temporal inspection.

panel |>
  group_by(date) |>
  summarise(
    cases = sum(cases),
    expected = sum(expected),
    sir_city = cases / expected,
    .groups = "drop"
  ) |>
  ggplot(aes(date, sir_city)) +
  geom_hline(yintercept = 1, linetype = 2, color = "#60717d") +
  geom_line(color = "#457b9d", linewidth = 1.1) +
  geom_point(color = "#2a9d8f", size = 2) +
  labs(
    title = "Aggregate SIR for Bandung City Simulation",
    x = NULL,
    y = "Aggregate SIR"
  )

This time-series plot collapses all Bandung districts into one monthly city-wide signal. It helps answer whether the simulated outbreak is mainly seasonal. If the aggregate SIR rises during rainy months, then a temporal component such as RW1/RW2 or seasonal terms is justified. If only a few districts spike while the city-wide line stays moderate, then spatial and space-time interaction components become more important.

Day 1 exercise. Keep the Bandung boundary fixed, then change the number of months, rainfall seasonality, and simulated outbreak pockets. Observe how SIR changes when population and case counts are small. Discuss why raw SIR should be smoothed before it is used as a policy signal.

Five R commands to begin with. After the first Bandung simulation chunk, students can inspect objects using basic R functions, before trying any model fitting.

head(panel[,c("district","date","cases","expected","sir")], 3)

## # A tibble: 3 × 5
##   district date       cases expected   sir
##   <chr>    <date>     <int>    <dbl> <dbl>
## 1 Andir    2024-01-01    33     10.6  3.12
## 2 Andir    2024-02-01    39     10.6  3.67
## 3 Andir    2024-03-01    32     10.7  2.99

mean(panel$cases)

## [1] 17.39722

length(unique(panel$district))

## [1] 30

Bayesian econometric concepts: a focused introduction

Learning goal. This section follows Henk Folmer’s recommended format. It introduces only those Bayesian ideas needed to understand Papers 1–3: Bayes’ rule, likelihood versus probability, prior/hyperprior/posterior, prior families, INLA in broad terms, prediction, and prediction with INLA. Conditional probability is assumed known. Detailed spatial/temporal priors, interaction types, diagnostics and computation appear with the respective papers.

Why Bayesian inference here? Paper 1 estimates and forecasts dengue risk across Bandung districts; Paper 2 estimates cluster-specific weather effects; Paper 3 assesses different forecast specifications. In each case we begin with a probability model, learn from observed cases, and report both estimates and uncertainty.

Bayes’ rule and the meaning of likelihood

\[ \boxed{\displaystyle \pi(\theta\mid y)= \frac{p(y\mid\theta)\pi(\theta)} {\int p(y\mid u)\pi(u)\,du} \propto L(\theta;y)\pi(\theta)} \]

The prior \(\pi(\theta)\) describes uncertainty before the current observations; the likelihood \(L(\theta;y)\) compares parameter values given the observed data; the posterior \(\pi(\theta\mid y)\) describes updated uncertainty. The denominator is the marginal likelihood or model evidence.

A probability distribution fixes \(\theta\) and varies the possible observation \(Y\). A likelihood fixes the observed \(y\) and varies \(\theta\). For four observed dengue cases:

\[ P(Y=4\mid\lambda)=\frac{\lambda^4e^{-\lambda}}{4!}, \qquad L(\lambda;4)\propto\lambda^4e^{-\lambda}. \]

For \(\lambda=2\), \(P(Y=4\mid\lambda)=0.0902\); for \(\lambda=4\), this probability is \(0.1954\). The data are more compatible with a mean of 4 than 2, but neither figure is the posterior probability that \(\lambda\) takes a particular value.

lambda <- seq(.1, 10, length.out = 200)
plot(lambda, dpois(4, lambda), type = "l", lwd = 2,
     col = "#AE6A3B", xlab = "Candidate count mean",
     ylab = "Likelihood", main = "Fixed observation y = 4")

round(c(mean_2 = dpois(4, 2), mean_4 = dpois(4, 4)), 4)

## mean_2 mean_4 
## 0.0902 0.1954

Prior, hyperprior, posterior: main types and examples

For a district-month Poisson model \(Y_{it}\mid\theta_{it}\sim\mathrm{Poisson}(E_{it}\theta_{it})\), the unknown \(\theta_{it}\) is relative dengue risk. A prior describes that risk. A hyperprior is a distribution on a parameter governing the prior, such as spatial effect standard deviation \(\sigma\) or precision \(\tau=1/\sigma^2\):

\[ \underbrace{p(y\mid\theta)}_{\mathrm{likelihood}}\, \underbrace{\pi(\theta\mid\tau)}_{\mathrm{prior}}\, \underbrace{\pi(\tau)}_{\mathrm{hyperprior}} \longrightarrow \underbrace{\pi(\theta,\tau\mid y)}_{\mathrm{posterior}}. \]

Type	Example	Why it is useful
Informative prior	\(\theta\sim\mathrm{Gamma}(20,20)\)	Strong prior concentration near risk 1
Weakly informative prior	\(\log\theta\sim N(0,1)\)	Regularization without a very narrow prior
Conjugate prior	\(\theta\sim\mathrm{Gamma}(a,b)\)	Transparent analytical updating with Poisson counts
Reference (Jeffreys) prior	\(\pi(\lambda)\propto\lambda^{-1/2}\) for Poisson mean	Objective-rule example; propriety must be checked
Gaussian effect prior	\(u_i\mid\sigma\sim N(0,\sigma^2)\)	Describe unobserved variability across districts
Scale hyperprior	\(\sigma\sim\mathrm{HalfNormal}(0,s)\)	Represent uncertainty about smoothing magnitude
Penalized-complexity hyperprior	\(P(\sigma>U)=\alpha\)	Calibrate preference for a simpler base model

The Gamma notation here is shape–rate, with mean \(a/b\). Other families include Beta priors for probabilities and normal priors for regression effects. CAR/BYM2, RW1/RW2 and interaction priors are explained when Paper 1 introduces their substantive roles.

A second example: a Gaussian conjugate prior. The introductory chapter of Bayesian Inference with INLA illustrates how posterior means balance a normal prior and an observed sample mean. Suppose \(Y_1,\ldots,Y_n\mid\mu\sim N(\mu,\tau^{-1})\), where the observation precision \(\tau\) is known, and the prior is \(\mu\sim N(\mu_0,\tau_0^{-1})\). Then

\[ \mu\mid\boldsymbol y\sim N\!\left(\mu_1,\tau_1^{-1}\right),\qquad \mu_1=\frac{\tau_0\mu_0+n\tau\bar y}{\tau_0+n\tau},\qquad \tau_1=\tau_0+n\tau. \]

For example, let the prior mean be \(10\), prior precision \(\tau_0=2\), sample mean \(\bar y=12\), \(n=5\), and sampling precision \(\tau=1\). The posterior mean is \((2\times10+5\times12)/7=11.429\) and the posterior standard deviation is \(1/\sqrt7=0.378\). If the sample size increases while other quantities remain fixed, the data receive more weight. Unlike a purely verbal explanation of ‘borrowing strength’, the formula makes that weighting explicit.

mu0 <- 10; prior_precision <- 2
ybar <- 12; n <- 5; sampling_precision <- 1
posterior_precision <- prior_precision + n * sampling_precision
posterior_mean <- (prior_precision * mu0 +
                     n * sampling_precision * ybar) / posterior_precision
round(c(posterior_mean = posterior_mean,
        posterior_sd = 1 / sqrt(posterior_precision)), 3)

## posterior_mean   posterior_sd 
##         11.429          0.378

The joint posterior can involve many unknown coefficients and random effects, while a marginal posterior focuses on one: \(\pi(\theta_j\mid y)=\int\pi(\boldsymbol\theta\mid y)\,d\boldsymbol\theta_{-j}\). Its mean, standard deviation, quantiles, and threshold probabilities are precisely the summaries that participants will later read from R-INLA output. This is the useful conceptual bridge from the reference e-book to our three applied papers.

Worked Bandung example. Suppose \(y=9\) dengue cases and \(E=7\) expected cases, with \(\theta\sim\mathrm{Gamma}(2,2)\). Then

\[ \underbrace{\mathrm{Gamma}(2,2)}_{\mathrm{prior}} \times \underbrace{\mathrm{Poisson}(7\theta)}_{\mathrm{likelihood}} \Longrightarrow \underbrace{\theta\mid y\sim\mathrm{Gamma}(11,9)}_{\mathrm{posterior}}. \]

The crude SIR is \(9/7=1.286\), prior mean \(1\), and posterior mean \(11/9=1.222\). A posterior distribution gives credible intervals and probabilities that relative risk exceeds a decision threshold, not only a point estimate.

y <- 9; E <- 7; a <- 2; b <- 2
rr <- seq(.01, 4, length.out = 300)
plot(rr, dgamma(rr, a, b), type = "l", lwd = 2,
     col = "#315B93", xlab = "Relative risk", ylab = "Density",
     main = "Prior, likelihood, and posterior")
lines(rr, dgamma(rr, y + 1, E), col = "#B98435", lwd = 2)
lines(rr, dgamma(rr, a + y, b + E), col = "#068E83", lwd = 3)
legend("topright", c("Prior", "Normalized likelihood", "Posterior"),
       col = c("#315B93", "#B98435", "#068E83"), lwd = 2, bty = "n")

round(c(crude_SIR = y/E, posterior_mean = (a+y)/(b+E),
        prob_RR_gt_1 = 1-pgamma(1, a+y, b+E)), 3)

##      crude_SIR posterior_mean   prob_RR_gt_1 
##          1.286          1.222          0.706

The likelihood curve is rescaled to integrate to one for visual comparison only; the underlying likelihood is not a posterior probability density. Use the interactive laboratory below to test whether prior influence is stronger with sparse observations.

LIVE BAYESIAN LAB · DAY 1

How much does the prior matter?

Prior shape, a: 2

Prior rate, b: 2

Observed cases, y: 9

Expected cases, E: 7

Prior mean1.00

Raw SIR1.29

Posterior mean1.22

P(RR > 1)—

Blue = Gamma prior; gold = normalized likelihood (Gamma(y+1, E) kernel); teal = Gamma posterior. The likelihood is divided by its area only to make its shape visually comparable with probability densities; likelihood is not intrinsically a probability density over RR. The posterior is Gamma(a+y, b+E). The graph updates immediately.

Estimation using INLA: the big picture

For a latent Gaussian spatiotemporal model, Integrated Nested Laplace Approximation (INLA) approximates marginal posterior distributions. In broad terms: choose a count likelihood; formulate a predictor with fixed and latent spatial–temporal effects; set priors and hyperpriors; run R-INLA; inspect fitted risks, uncertainty and posterior exceedance probabilities. It combines analytical approximations and numerical integration. It is not MCMC, and it is not exact for every model.

Participants should install R-INLA once before knitting using the instructions at https://www.r-inla.org/download/; software installation is deliberately not performed by any code chunk during rendering.

The specific latent fields, spatial adjacency, interaction models and numerical estimation are addressed in Day 2 (Paper 1), followed by cluster-varying effects in Day 3 and forecasting comparisons in Day 4. No Laplace-approximation derivations are needed at this foundational stage.

Definition of prediction and prediction with INLA

Estimation concerns an unknown model parameter or latent risk. Prediction concerns an unobserved outcome \(Y_*\) or latent risk \(\theta_*\). Forecasting is prediction for a future period. Bayesian predictive inference averages over parameter uncertainty:

\[ \boxed{\displaystyle p(y_*\mid y)=\int p(y_*\mid\theta)\,\pi(\theta\mid y)\,d\theta. } \]

A plug-in predicted count is not the same as the posterior predictive distribution: the latter accounts for uncertainty in both the data-generating process and inferred parameters.

To predict with INLA, add a row whose outcome is missing, provide its covariates and exposure, then compute fitted-value/posterior summaries. The following small count example uses a simple trend to explain the interface—not the full spatial model used in Paper 1:

city_month <- panel |>
  dplyr::group_by(t) |>
  dplyr::summarise(y = sum(cases), E = sum(expected), .groups = "drop")
next_month <- data.frame(t = max(city_month$t) + 1L,
                         y = NA_integer_, E = tail(city_month$E, 1))
new_data <- rbind(city_month, next_month)
fit_day1_demo <- INLA::inla(
  y ~ 1 + t, family = "poisson", E = E, data = new_data,
  control.predictor = list(compute = TRUE))
round(tail(fit_day1_demo$summary.fitted.values[
  , c("mean", "0.025quant", "0.975quant")], 1), 2)

##                      mean 0.025quant 0.975quant
## fitted.Predictor.37 -0.21      -0.24      -0.18

The last row describes the future month. In a genuine spatiotemporal model, future latent effects also need a specified temporal evolution and valid spatial indexing; merely adding an NA row does not automatically provide these. Day 2 explains posterior risk mapping, and Day 4 compares held-out forecasts.

Day 1 discussion. What is the difference between a likelihood and a posterior? Which hyperprior could control district-to-district smoothing? Why can two priors change an exceedance map? How does forecasting differ from estimating risk in an already observed month?

Day 2

Paper 1 · Pure Disease Mapping

Spatiotemporal risk, Poisson and negative binomial models, INLA estimation and Bandung risk maps.

Day 2 - Paper 1: Pure Bayesian Spatiotemporal Disease Mapping

Bandung illustration. Compare the raw district SIR against the smoothed posterior relative risk. Interpret a 95% credible interval and explain where a high-risk map has strong or weak evidence.

Source coverage. This chapter follows all 65 slides of the separately supplied Paper 1 presentation: motivation and epidemiological quantities (slides 1–19), priors and four interaction structures (20–38), INLA and model comparison (39–49), and the published Bandung application, graphics, forecasts and conclusions (50–65). Mathematical content is typeset below, while original scientific graphs are included individually—not as screenshots of the slides.

From observed cases to relative risk. Set \(E_{it}=N_{it}p_t\), with population \(N_{it}\) and a time-specific reference incidence \(p_t\). The raw SIR is \(Y_{it}/E_{it}\). The model assumes \(Y_{it}\mid R_{it}\sim\mathrm{Poisson}(E_{it}R_{it})\), \(\log R_{it}=\alpha+u_i+v_i+\gamma_t+\delta_t+\psi_{it}\). In that decomposition, \(u_i\) is spatially structured, \(v_i\) IID spatial heterogeneity, \(\gamma_t\) structured time, \(\delta_t\) IID time, and \(\psi_{it}\) the interaction.

Count-family choice. Poisson imposes conditional \(\operatorname{Var}(Y_{it}\mid R)=E_{it}R_{it}\). The negative binomial allows \(\operatorname{Var}(Y_{it}\mid\mu,\kappa)=\mu+\mu^2/\kappa\). Random effects can themselves create marginal overdispersion; the two explanations should be compared rather than conflated.

Worked variance example. With expected mean \(\mu=12\) and negative-binomial size \(\kappa=4\), the Poisson variance is \(12\) but the negative-binomial variance is \(12+12^2/4=48\).

mu <- 12; kappa <- 4
data.frame(model=c("Poisson","Negative binomial"),
           mean=c(mu,mu),variance=c(mu,mu+mu^2/kappa))

##               model mean variance
## 1           Poisson   12       12
## 2 Negative binomial   12       48

Structured effects and prediction. For ICAR, the conditional average of \(u_i\) is influenced by neighbors; RW1 penalizes first temporal differences; RW2 penalizes second differences. Knorr-Held types I–IV cross structured/unstructured space with structured/unstructured time. Identifiability constraints are necessary for intrinsic priors. The original study compares candidate specifications and interprets estimated risks, seasonal patterns, exceedance probabilities and geographic variation. Continuous isopleth visualization via inverse-distance weighting is a post-estimation display technique; it does not change the underlying district-level estimand.

Assessing fit and predictive value. DIC and WAIC summarize fit with complexity penalties; CPO and PIT can reveal poor predictive performance or calibration. Posterior exceedance \(P(RR_{it}>1.25\mid y)\) directly communicates risk-threshold evidence. Out-of-sample forecasting should reserve future observations and never select models using those outcomes.

Core sources for this session: SIR and disease mapping follow Clayton and Kaldor (1987), Waller and Gotway (2004), and Lawson (2018). CAR priors and spatiotemporal interaction models follow Besag et al. (1991), Knorr-Held (2000), Rue and Held (2005), and Blangiardo and Cameletti (2015).

The standardized incidence ratio, Poisson model, and relative risk as a pure model

The SIR and Poisson relative-risk formulation are standard tools in disease mapping, especially when observed cases are compared with expected cases through an offset (Clayton and Kaldor, 1987; Waller and Gotway, 2004; Lawson, 2018).

Beginner model-building ladder.

The pure model can be understood as a ladder:

Step	Model idea	What is added?
1	Raw counts	Only observed cases
2	Expected cases	Population adjustment
3	SIR	Crude relative risk
4	Poisson model	Probability model for counts
5	Spatial effect	Neighboring areas can share risk
6	Temporal effect	Adjacent months can share risk
7	Space-time interaction	Local changes can differ over time
8	Posterior risk map	Smoothed risk with uncertainty

This ladder is useful because beginners can diagnose where they are in the modeling process. Do not jump to Step 8 before understanding Steps 1 to 4.

The standardized incidence ratio is:

\[ \mathrm{SIR}_{it}=\frac{Y_{it}}{E_{it}} \]

It is the crude estimate of relative risk. The model-based target is the latent relative risk:

\[ \theta_{it}=\frac{\mu_{it}}{E_{it}}, \qquad \mu_{it}=\mathbb{E}(Y_{it}) \]

Under the Poisson disease model:

\[ Y_{it}\mid \theta_{it} \sim \mathrm{Poisson}(E_{it}\theta_{it}) \]

The log link gives:

\[ \log(\mu_{it})=\log(E_{it})+\log(\theta_{it}) \]

The term \(\log(E_{it})\) is an offset. Its coefficient is fixed at 1, so the model estimates log-relative risk rather than raw case counts.

The pure model decomposes log-relative risk as:

\[ \eta_{it} = \log(\theta_{it}) = \alpha + u_i + v_i + \gamma_t + \phi_t + \delta_{it} \]

where:

• \(\alpha\): global average risk on the log scale.
• \(u_i\): structured spatial effect, capturing autocorrelation among neighboring areas.
• \(v_i\): unstructured spatial effect, capturing area-level heterogeneity.
• \(\gamma_t\): structured temporal effect, such as a random walk.
• \(\phi_t\): unstructured temporal effect.
• \(\delta_{it}\): space-time interaction.

Each component has a specific teaching interpretation:

Component	Statistical_role	Public_health_interpretation
\(\alpha\)	Global intercept	Overall average log-risk across all areas and periods
\(u_i\)	Structured spatial random effect	Risk shared by neighboring areas, such as environmental similarity or local transmission
\(v_i\)	Unstructured spatial random effect	Area-specific heterogeneity not explained by adjacency
\(\gamma_t\)	Structured temporal random effect	Smooth epidemic trend or seasonality over time
\(\phi_t\)	Unstructured temporal random effect	Sudden monthly or weekly shocks
\(\delta_{it}\)	Space-time interaction	Local outbreaks that deviate from the main spatial and temporal patterns

Likelihood function of the Poisson incidence model.

For the full Bandung panel, the likelihood is:

\[ L(\Theta,\Lambda;y) = \prod_{i=1}^{n}\prod_{t=1}^{T} \frac{(E_{it}\theta_{it})^{Y_{it}}\exp(-E_{it}\theta_{it})}{Y_{it}!} \]

where:

\[ \theta_{it} = \exp(\alpha+u_i+v_i+\gamma_t+\phi_t+\delta_{it}) \]

This likelihood says that observed cases are noisy realizations around expected cases multiplied by latent relative risk. The Bayesian model estimates the latent risk surface rather than treating the crude SIR as the final answer.

Prior distributions of spatial, temporal, and interaction effects

The structured and unstructured random-effect priors used in disease mapping are grounded in CAR/GMRF and spatiotemporal interaction modeling traditions (Besag et al., 1991; Knorr-Held, 2000; Rue and Held, 2005).

Beginner intuition for random-effect priors.

Random-effect priors define how information is shared. They are not only mathematical details; they determine the type of smoothing.

Prior component	What it smooths	Beginner interpretation
Structured spatial effect	Neighboring areas	Nearby districts borrow information from each other
Unstructured spatial effect	Area-specific noise	A district may have unique risk not explained by neighbors
Structured temporal effect	Adjacent periods	Consecutive months borrow information from each other
Unstructured temporal effect	Time-specific shocks	One month may have a city-wide unusual shock
Space-time interaction	Local changes over time	A district can deviate from the general spatial and temporal patterns

A good disease-mapping model often needs both structured and unstructured components because real surveillance data contain both smooth patterns and unexpected shocks.

Intrinsic CAR prior for structured spatial effects.

\[ p(u\mid \tau_u) \propto \tau_u^{(n-1)/2} \exp\left[ -\frac{\tau_u}{2} \sum_{i\sim j}(u_i-u_j)^2 \right] \]

Leroux CAR prior.

\[ u_i\mid u_{-i}, \rho,\tau_u \sim N\left( \frac{\rho\sum_j w_{ij}u_j}{\rho\sum_jw_{ij}+1-\rho}, \frac{1}{\tau_u(\rho\sum_jw_{ij}+1-\rho)} \right) \]

If \(\rho\) is close to 1, spatial structure dominates. If \(\rho\) is close to 0, the effects become more independent.

Temporal priors.

First-order random walk:

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

Second-order random walk:

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

Unstructured temporal effects:

\[ \phi_t\sim N(0,\tau_\phi^{-1}) \]

Knorr-Held space-time interaction priors.

The interaction \(\delta_{it}\) can be formulated in four types:

1. Type I: unstructured spatial x unstructured temporal.
2. Type II: structured spatial x unstructured temporal.
3. Type III: unstructured spatial x structured temporal.
4. Type IV: structured spatial x structured temporal.

Using precision-matrix notation:

\[ Q_\delta = Q_t \otimes Q_s \]

where \(\otimes\) is the Kronecker product. The interaction type is selected using DIC, WAIC, CPO/PIT, or predictive performance.

Joint posterior distribution and marginal posterior distribution.

\[ p(\Theta,\Lambda\mid y) \propto p(y\mid \Theta,\Lambda)p(\Theta\mid \Lambda)p(\Lambda) \]

where \(\Theta\) denotes the main parameters and \(\Lambda\) denotes hyperparameters. In INLA, the main targets are marginal posteriors:

\[ p(\theta_k\mid y) = \int p(\theta_k\mid \lambda,y)p(\lambda\mid y)d\lambda \]

The marginal posterior \(p(\theta_{it}\mid y)\) is the main object for risk mapping. From this distribution, we can report posterior mean relative risk, credible intervals, and posterior exceedance probabilities such as \(\Pr(\theta_{it}>1.25\mid y)\).

Applied example: quick baseline model for Bandung.

The following GLM is not a full replacement for the pure Bayesian model, but it is useful for classroom demonstration: it shows the offset, district effects, and seasonal effects.

The applied workflow in this session is:

1. Fit a simple offset model to reinforce the role of expected cases.
2. Compare observed SIR with fitted relative risk.
3. Introduce the INLA template for the full Bayesian pure model.
4. Extend the data with future rows for forecasting.
5. Convert forecasted risk into a hotspot probability.

glm_pure <- glm(
  cases ~ offset(log(expected)) + factor(district_id) + factor(month_no),
  family = poisson(),
  data = panel
)

panel <- panel |>
  mutate(
    mu_glm = predict(glm_pure, type = "response"),
    rr_glm = mu_glm / expected
  )

overdispersion <- sum(residuals(glm_pure, type = "pearson")^2) / glm_pure$df.residual
overdispersion

## [1] 1.895369

The overdispersion statistic is a quick diagnostic. A value much larger than 1 indicates that the Poisson variance assumption may be too restrictive, motivating a negative binomial model or richer random effects.

Model	AIC	Interpretation
Poisson GLM	6734.65	Baseline count model with mean equal to variance
Negative binomial GLM	6452.37	Allows variance to exceed the mean

This comparison is included for teaching purposes. In the thesis framework, overdispersion can also be handled through random effects and Bayesian hierarchical modeling, not only by changing the likelihood family.

ggplot(panel, aes(sir, rr_glm)) +
  geom_point(alpha = 0.45, color = "#457b9d") +
  geom_abline(slope = 1, intercept = 0, linetype = 2, color = "#e76f51") +
  coord_equal() +
  labs(
    title = "Observed SIR vs GLM-Based Relative Risk",
    x = "Observed SIR",
    y = "Model-based relative risk"
  )

The scatterplot compares crude SIR with model-smoothed relative risk. Points far from the diagonal show district-months where the model moderates noisy crude risk. In practice, this moderation is useful when a small denominator produces an extreme SIR that is not supported by the broader spatial-temporal pattern.

Minimal R-INLA pattern. The full spatial and temporal model is computed by the following background chunk. This compact example shows the important formula arguments without overwhelming new R learners.

library(INLA)
fit <- inla(cases ~ 1, family="poisson", E=expected,
            data=panel, control.predictor=list(compute=TRUE))
fit$summary.fixed

##                       mean          sd  0.025quant       0.5quant 0.975quant
## (Intercept) -0.00002655231 0.007295392 -0.01432526 -0.00002655231 0.01427215
##                       mode kld
## (Intercept) -0.00002655231   0

Estimation using INLA and marginal posterior extraction

Beginner interpretation of INLA output.

INLA output can look large and technical. Beginners should focus first on four outputs:

1. Fixed effects: baseline level and covariate effects if covariates are included.
2. Hyperparameters: smoothing strength and variance/precision components.
3. Fitted values: estimated mean cases or risk for each observation.
4. Marginal posteriors: uncertainty distributions for parameters and latent effects.

For disease mapping, the most important interpretation is usually not a single coefficient. It is the map of posterior mean relative risk and the map of posterior probability that risk exceeds an action threshold.

This workshop fits a compact BYM2 + RW1 model. The more elaborate spatial, temporal and interaction alternatives are explained above, but only the compact specification is fitted in the classroom to keep computation manageable and the model interpretable.

##        DIC     WAIC
## 1 6701.863 6769.591

After fitting the INLA model, the fitted values are on the expected case scale. Relative risk is obtained by dividing fitted means by expected cases:

In a production analysis, compute exceedance probabilities directly from the marginal posterior of \(\theta_{it}\) whenever possible. The main point is that uncertainty is not an optional add-on; it is part of the disease mapping result.

Forecasting procedure and mapping

In INLA, forecasting is implemented by adding future rows and setting cases = NA. The posterior predictive distribution is then computed for the missing observations. For the quick exercise, we use the offset GLM already fitted above.

date	mean_rr	max_rr	hotspot_count
2027-01-01	1.587	3.451	13
2027-02-01	2.064	4.490	20
2027-03-01	2.009	4.369	20
2027-04-01	1.571	3.417	13
2027-05-01	1.042	2.267	5
2027-06-01	0.624	1.357	0

This forecast table summarizes Bandung-wide risk for each future month. mean_rr is the average predicted relative risk across the 30 districts, max_rr identifies the strongest predicted local risk, and hotspot_count counts districts whose exceedance probability is above the operational threshold. A rising hotspot_count is more actionable than a small change in the city-wide mean because it tells the health office how many districts may need field response.

future_map <- districts |>
  select(district_id, district) |>
  left_join(
    future |>
      filter(h %in% c(1, 3, 6)) |>
      mutate(month_label = format(date, "%Y-%m")),
    by = "district_id"
  )

ggplot(future_map) +
  geom_sf(aes(fill = rr_forecast), color = "white", linewidth = 0.25) +
  facet_wrap(~month_label, nrow = 1) +
  scale_fill_gradient2(
    low = "#f1d9b1", mid = "#fffaf3", high = "#973d1c",
    midpoint = 1, name = "Forecast RR"
  ) +
  labs(
    title = "Relative Risk Forecast",
    subtitle = "Example forecasts 1, 3, and 6 months ahead"
  )

The forecast map shows where future risk is expected to concentrate inside Bandung. In a real early warning system, this map would be used to pre-position staff, larvicide, fogging capacity, public communication, and environmental inspection teams before the peak-risk month occurs.

Model selection and predictive diagnostics.

Criteria used in the thesis:

\[ \mathrm{DIC} = \bar{D} + p_D \]

\[ \mathrm{WAIC} = -2(\mathrm{lppd} - p_{\mathrm{WAIC}}) \]

\[ \mathrm{MPL} = \sum_{i,t}\log(\mathrm{CPO}_{it}) \]

\[ \mathrm{MAE} = \frac{1}{nT}\sum_{i,t}|Y_{it}-\widehat{Y}_{it}|, \qquad \mathrm{RMSE} = \sqrt{\frac{1}{nT}\sum_{i,t}(Y_{it}-\widehat{Y}_{it})^2} \]

Practical principles:

• Smaller DIC/WAIC values are generally preferred.
• Larger MPL values are generally preferred.
• A PIT distribution close to uniform indicates better predictive calibration.
• Model selection should not be purely mechanical; epidemiological interpretation must remain part of the decision.

Day 2 exercise. Run the INLA template if the package and computation time are available. Compare Poisson, negative binomial, RW1, RW2, and interaction-type specifications. Create a table of DIC, WAIC, MAE, and RMSE.

Paper 1 in depth: likelihood, priors, interaction, and INLA

Scientific motivation (PPTX slides 1–4). Disease surveillance produces dependent observations in both space and time. Population sizes differ, neighboring districts share environmental and social conditions, and dengue occurrence changes with season and outbreaks. A pure model explains this variability using latent spatial, temporal and space–time effects, without requiring observed weather regressors. This differs from Day 3, where explanatory variables and changing weather coefficients are explicitly modeled. Both a district-level choropleth and a continuous isopleth map can display risk, but smoothing map colors does not turn areal observations into individual-level measurements.

Expected cases and the offset (slides 5–14). The publication explicitly contrasts a period-specific standard with a pooled standard:

\[ p_t=\frac{\sum_{i=1}^nY_{it}}{\sum_{i=1}^nN_{it}},\qquad \bar p=\frac{\sum_{i=1}^n\sum_{t=1}^T Y_{it}}{\sum_{i=1}^n\sum_{t=1}^T N_{it}}, \qquad E_{it}^{(t)}=N_{it}p_t,\quad E_{it}^{(\mathrm{pool})}=N_{it}\bar p. \]

The published study chooses the pooled rate so that a high month is not automatically normalized to a relative risk near one. The choice must be declared before comparisons are made. For a classroom calculation, suppose three Bandung districts have populations 50,000, 80,000 and 70,000 and 20, 48 and 42 monthly cases. The pooled one-month incidence is \(110/200{,}000=0.00055\). Therefore the first district has \(E_1=27.5\) cases and \(\mathrm{SIR}_1=20/27.5=0.7273\). The second district has \(E_2=44\) and \(\mathrm{SIR}_2=48/44=1.0909\).

pop <- c(50000, 80000, 70000)
cases_example <- c(20, 48, 42)
p_reference <- sum(cases_example) / sum(pop)
expected_example <- pop * p_reference
sir_example <- cases_example / expected_example
round(c(reference_rate = p_reference, expected_1 = expected_example[1],
        SIR_1 = sir_example[1], expected_2 = expected_example[2],
        SIR_2 = sir_example[2]), 5)

## reference_rate     expected_1          SIR_1     expected_2          SIR_2 
##        0.00055       27.50000        0.72727       44.00000        1.09091

The offset enters the sampling model as \(\log\mu_{it}=\log E_{it}+\eta_{it}\), where \(\theta_{it}=\exp(\eta_{it})\). Consequently, \(\mu_{it}=E_{it}\theta_{it}\) and the estimated effect is interpretable as relative risk, not simply a difference in raw case counts.

Poisson likelihood and gamma–Poisson mixing (slides 15–20). For counts \(Y_{it}\), let \(\mu_{it}=E_{it}\theta_{it}\). The simple model is

\[ P(Y_{it}=y\mid\mu_{it})=\frac{\mu_{it}^{y}e^{-\mu_{it}}}{y!}, \qquad \mathbb E(Y_{it})=\operatorname{Var}(Y_{it})=\mu_{it}. \]

When the variance exceeds the mean even after accounting for main effects, use \(\varepsilon_{it}\sim\mathrm{Gamma}(\kappa,\kappa)\) in shape–rate parameterization and \(Y_{it}\mid\varepsilon_{it}\sim\mathrm{Poisson}(\mu_{it}\varepsilon_{it})\). Integrating out \(\varepsilon_{it}\) gives

\[ P(Y_{it}=y\mid\mu_{it},\kappa)= \frac{\Gamma(y+\kappa)}{\Gamma(\kappa)y!} \left(\frac{\kappa}{\kappa+\mu_{it}}\right)^\kappa \left(\frac{\mu_{it}}{\kappa+\mu_{it}}\right)^y, \qquad \operatorname{Var}(Y_{it})=\mu_{it}+\frac{\mu_{it}^2}{\kappa}. \]

This is NB2. As \(\kappa\to\infty\), the negative binomial approaches the Poisson. Note that an alternative overdispersion parameter \(a=1/\kappa\) gives \(\operatorname{Var}(Y)=\mu+a\mu^2\). The paper reports the latter type of extra-Poisson variability; do not equate an estimated \(a\) with the size \(\kappa\).

mean_cases <- 12
nb_size <- 4
y_grid <- 0:38
plot(y_grid, dpois(y_grid, lambda = mean_cases), type = "h", lwd = 3,
     col = "#4978A8", xlab = "Cases", ylab = "Probability",
     main = "Same mean, different count variability")
points(y_grid + .18, dnbinom(y_grid, size = nb_size, mu = mean_cases),
       type = "h", lwd = 3, col = "#CA6A38")
legend("topright", legend = c("Poisson", "Negative binomial"),
       col = c("#4978A8", "#CA6A38"), lwd = 3, bty = "n")

round(c(poisson_variance = mean_cases,
        nb_variance = mean_cases + mean_cases^2 / nb_size), 2)

## poisson_variance      nb_variance 
##               12               48

Full latent hierarchy (slides 21–24, 30–35 and 37–38). Using the paper’s symbols, the pure log-relative risk is

\[ \eta_{it}=\log\theta_{it} =\alpha+\omega_i+\nu_i+\phi_t+\gamma_t+\delta_{it}, \]

where \(\omega_i\) is structured spatial heterogeneity; \(\nu_i\) is independent spatial heterogeneity; \(\phi_t\) is structured temporal evolution; \(\gamma_t\) is independent temporal noise; and \(\delta_{it}\) is their interaction. A seasonal \(\kappa_{m(t)}\) may be added when a recurring month-of-year effect is modeled separately from the multi-year trend; it must not be confused with the negative-binomial size parameter above.

In hierarchical form, the three levels are:

\[ \underbrace{Y_{it}\mid\eta_{it},\kappa}_{\text{data model}} \quad\longleftarrow\quad \underbrace{\boldsymbol\eta\mid\boldsymbol\tau,\rho}_{\text{latent effects}} \quad\longleftarrow\quad \underbrace{\pi(\boldsymbol\tau,\rho,\kappa)}_{\text{hyperpriors}}. \]

The spatial ICAR conditional specification is

\[ \omega_i\mid\boldsymbol\omega_{-i},\tau_\omega \sim N\!\left(\frac{1}{n_i}\sum_{j\sim i}\omega_j, \frac{1}{\tau_\omega n_i}\right),\qquad Q_\omega=\tau_\omega(D-W), \]

with \(n_i\) neighbors, adjacency matrix \(W\), and degree matrix \(D\). A sum-to-zero constraint is needed because \(D-W\) is singular for a connected graph. The independent alternative is \(\nu_i\overset{\mathrm{iid}}\sim N(0,\tau_\nu^{-1})\). A Leroux variant replaces the intrinsic precision by \(Q_\rho=\tau[(1-\rho)I+\rho(D-W)]\) for \(0\leq\rho<1\). These priors say what kind of risk patterns are plausible before seeing the observed cases.

Random-walk trends, seasonal patterns and interaction (slides 25–36). The temporal priors penalize differences rather than penalizing absolute levels:

\[ \pi(\boldsymbol\phi\mid\tau_\phi,\mathrm{RW1})\propto \exp\!\left[-\frac{\tau_\phi}{2}\sum_{t=2}^T(\phi_t-\phi_{t-1})^2\right], \]

\[ \pi(\boldsymbol\phi\mid\tau_\phi,\mathrm{RW2})\propto \exp\!\left[-\frac{\tau_\phi}{2}\sum_{t=3}^T (\phi_t-2\phi_{t-1}+\phi_{t-2})^2\right]. \]

RW1 encourages adjacent months to have similar risks. RW2 encourages adjacent slopes to be similar, allowing a smoother changing trend. A cyclic annual seasonal effect, if fitted, links December back to January and is not interchangeable with a noncyclic RW1 on all months. The figures below simulate the two priors directly: no packages beyond base R are required.

set.seed(208)
T_demo <- 36
rw1 <- cumsum(rnorm(T_demo, sd = .12))
rw2 <- cumsum(cumsum(rnorm(T_demo, sd = .024)))
rw1 <- rw1 - mean(rw1)
rw2 <- rw2 - mean(rw2)
plot(seq_len(T_demo), rw1, type = "l", lwd = 2, col = "#3575B8",
     xlab = "Month", ylab = "Centered log-risk",
     ylim = range(rw1, rw2), main = "How different temporal priors shape trends")
lines(seq_len(T_demo), rw2, col = "#C46B35", lwd = 2)
legend("topleft", c("RW1: smooth levels", "RW2: smooth slopes"),
       col = c("#3575B8", "#C46B35"), lwd = 2, bty = "n")

The four Knorr-Held interactions differ in which dimension is correlated. If \(I_n\) and \(I_T\) are identity matrices, \(Q_s\) is a spatial ICAR precision, and \(Q_t\) is a temporal RW precision, then (up to the indexing convention and identifiability constraints)

\[ \begin{aligned} \mathrm{I:}\;&Q_{\delta}=\tau_\delta(I_T\otimes I_n),\\ \mathrm{II:}\;&Q_{\delta}=\tau_\delta(I_T\otimes Q_s),\\ \mathrm{III:}\;&Q_{\delta}=\tau_\delta(Q_t\otimes I_n),\\ \mathrm{IV:}\;&Q_{\delta}=\tau_\delta(Q_t\otimes Q_s). \end{aligned} \]

Type IV simultaneously shares information among neighboring districts and adjacent months. It does not assert that every neighbor evolves identically. An interaction is a local departure from the sum of the main spatial and temporal effects; if a structured interaction is added, appropriate constraints are needed to separate it from those main effects.

# Illustrate persistence of a spatially localized outbreak over three months.
illustration_map <- districts |>
  dplyr::mutate(base = as.numeric(scale(sf::st_coordinates(sf::st_centroid(sf::st_geometry(districts)))[, 1])))
illustration_map <- tidyr::crossing(dplyr::select(illustration_map, district_id, district, base), month = 1:3) |>
  dplyr::mutate(effect = .25 * base + .35 * sin(month / 3 + base))
ggplot2::ggplot(dplyr::left_join(districts, illustration_map,
                                 by = c("district_id", "district"))) +
  ggplot2::geom_sf(ggplot2::aes(fill = effect), colour = "white", linewidth = .15) +
  ggplot2::facet_wrap(~month, nrow = 1) +
  ggplot2::scale_fill_gradient2(low = "#3D7B9E", mid = "white", high = "#BD5D31") +
  ggplot2::labs(title = "Simulated correlated local effects on Bandung boundaries",
                subtitle = "Conceptual illustration—not fitted Type IV posterior effects")

INLA estimation (slides 39–46). A latent Gaussian model has \(\boldsymbol x\mid\boldsymbol\vartheta\sim N(0,Q(\boldsymbol\vartheta)^{-1})\) (with an appropriate mean or fixed effects). Bayes’ rule gives

\[ \pi(\boldsymbol x,\boldsymbol\vartheta\mid\boldsymbol y) \propto\pi(\boldsymbol y\mid\boldsymbol x,\boldsymbol\vartheta) \pi(\boldsymbol x\mid\boldsymbol\vartheta)\pi(\boldsymbol\vartheta), \]

and INLA approximates hyperparameter posteriors and integrates them numerically:

\[ \widetilde\pi(x_j\mid\boldsymbol y) =\int\widetilde\pi(x_j\mid\boldsymbol\vartheta,\boldsymbol y) \widetilde\pi(\boldsymbol\vartheta\mid\boldsymbol y)\,d\boldsymbol\vartheta. \]

This is why the workflow yields a whole marginal posterior for log-risk rather than just a maximum-likelihood point estimate. A normal approximation to a log-risk marginal with posterior mean \(m\) and standard deviation \(s\) yields an instructive numerical exceedance calculation:

\[ \Pr(\theta>1.25\mid y) =\Pr(\eta>\log 1.25\mid y) \approx 1-\Phi\!\left(\frac{\log(1.25)-m}{s}\right). \]

post_mean_logrisk <- log(1.40)
post_sd_logrisk <- .20
probability <- 1 - pnorm(log(1.25), post_mean_logrisk, post_sd_logrisk)
round(c(posterior_mean_RR_approx = exp(post_mean_logrisk + post_sd_logrisk^2 / 2),
        probability_RR_above_1_25 = probability), 3)

##  posterior_mean_RR_approx probability_RR_above_1_25 
##                     1.428                     0.715

For actual analysis, compute exceedance from the numerical marginal density, not the normal approximation, where available. INLA also returns posterior summaries of fixed effects, hyperparameters, fitted means, CPO and PIT. The original publication’s forecasters append missing future responses and obtain region-by-month predictive posteriors; predictions cannot be validated using future outcomes that were also used for fitting.

Paper 1 model comparison, historical results, and original scientific figures

Model assessment (slides 47–49). The paper considers Poisson versus negative-binomial sampling, RW1 versus RW2 temporal structure, and Type I–IV interactions. Its selection uses DIC, WAIC, conditional predictive ordinate (CPO), marginal predictive likelihood (MPL), pseudo-Bayes factors and prediction errors. For observation \(k\), \(\mathrm{CPO}_k=p(y_k\mid\boldsymbol y_{-k})\) and \(\mathrm{MPL}=\sum_k\log\mathrm{CPO}_k\). Predictive integral transform values are \(\mathrm{PIT}_k=F(y_k\mid\boldsymbol y_{-k})\); because counts are discrete, PIT interpretation should account for discreteness (for example via randomized PIT). A plot of PIT values and residual autocorrelations tests different weaknesses than DIC.

Published data and likelihood decision (slides 50–52). The original study is based on 30 Bandung districts, 2009–2016, with 35,496 registered dengue cases. It reports an extra-Poisson dispersion estimate of \(0.0238\) (standard error \(0.0038\), 95% interval \(0.0177\)–\(0.0327\)), supporting the negative-binomial likelihood. The pooled reference rate used in the publication differs from the more pedagogical month-specific reference introduced earlier. None of these historical numbers were generated from the synthetic dataset in this e-book.

Study-area geography reproduced from the Paper 1 presentation

Published interaction comparison (slides 52–54). With the stated inverse-gamma hyperprior configuration, Type IV has the lowest reported DIC and WAIC and highest MPL among the four interaction specifications. It is important to compare the same likelihood, response and dataset across these values; the simulation’s own DIC values are not directly comparable to published estimates.

Interaction	DIC	WAIC	MPL	MAE	RMSE
Type I	17,490.21	17,500.95	−8,753.33	4.37	6.27
Type II	16,552.63	16,498.89	−8,264.53	3.15	4.44
Type III	16,247.67	16,217.70	−8,299.80	2.19	2.93
Type IV	15,866.68	15,910.80	−8,004.75	2.58	3.53

Type III gives lower MAE and RMSE than Type IV in the article’s table, so no model wins every criterion. This is an important workshop discussion: in-sample absolute error, posterior predictive scoring and effective model complexity answer different questions. The original authors give primary weight to their Bayesian selection criteria and continue with Type IV.

The paper’s standard-deviation summaries show posterior mean SDs for unstructured space \(0.147\), spatial ICAR \(0.148\), unstructured time \(0.097\), RW1 trend \(0.560\), annual seasonal effect \(0.059\) and interaction \(0.174\). Their respective reported critical ratios (estimate / posterior standard error) are approximately \(1.56\), \(1.64\), \(3.29\), \(6.54\), \(4.38\) and \(20.08\). These are descriptive ratios rather than posterior probabilities or formal Bayesian credible-interval tests. The relatively weak separate spatial main effects motivate the reduced temporal specification retained in the paper, while the structured interaction still allows spatially varying evolution.

Published temporal and district-level graphics (slides 55–58). The original article analyzes multi-year monthly structured effects, recurring within-year seasonality, and district deviations. The figures below are original scientific figures extracted individually from the PowerPoint—no full-slide screenshots. The e-book’s runnable Bandung illustrations above remain independently editable R plots.

Paper 1: structured temporal and uncertainty effects, original scientific figure

The temporal plot distinguishes estimated mean trend from a credible interval. Strong seasonality should be interpreted on the scale and baseline indicated by the figure; it is not equivalent to a trend in raw city-wide counts.

Paper 1: estimated annual seasonal dengue pattern, original scientific figure

The annual cycle suggests elevated risk around the first part of the year and reduced risk later in the year. The model additionally allows each district to depart from this shared pattern.

Paper 1: spatiotemporal interaction in selected districts, original scientific figure

Historical risk maps (slide 59). The paper’s 2009 and 2016 maps show that the geography of elevated dengue relative risk can move over time: higher-risk patterns were concentrated in northern Bandung in 2009 and shifted toward southern districts by 2016. Both maps represent historical posterior summaries, not simulated case counts.

Paper 1: published dengue relative risk in 2009

Paper 1: published dengue relative risk in 2016

Time-series forecasts (slide 60). The publication compares observed and estimated relative risk for 2009–2016 and projects 2017 for selected districts. These curves show why a single Bandung-wide trend is insufficient: Sumur Bandung, Cidadap and other districts have different baseline risks and prospective paths.

Paper 1: historical observation, fitted risk and 2017 forecast by district

Forecast risk and exceedance maps (slides 61–62). A forecasted posterior mean and the probability of exceeding the operational threshold are different maps. The first answers “how high might risk be?”; the second answers “how certain are we that risk exceeds 1.25?” Both must be interpreted alongside the lead time, exposure baseline and uncertainty.

Paper 1: monthly posterior mean dengue relative risk forecasts for 2017

Paper 1: monthly posterior exceedance probability of RR above 1.25 for 2017

Recompute an exceedance map with our simulated Bandung observations. The model fitted earlier yields a posterior marginal for every district-month. The code below creates a new simulated-data probability map from the approximated marginal log-risk. The original 2017 map above is preserved only as the publication’s research evidence.

last_simulated_month <- max(dat_inla$t)
exceedance_map <- districts |>
  dplyr::select(district_id, district) |>
  dplyr::left_join(dat_inla |>
                     dplyr::filter(t == last_simulated_month) |>
                     dplyr::select(district_id, hotspot_prob_125),
                   by = "district_id")
ggplot2::ggplot(exceedance_map) +
  ggplot2::geom_sf(ggplot2::aes(fill = hotspot_prob_125),
                   colour = "white", linewidth = .2) +
  ggplot2::scale_fill_gradient(low = "#FFF7E8", high = "#A74225",
                               limits = c(0, 1), name = "P(RR > 1.25)") +
  ggplot2::labs(title = "Simulated Bandung: posterior risk exceedance",
                subtitle = "Most recent simulated month; Gaussian marginal approximation")

Why this matters for early warning (slides 63–65). The authors’ conclusion is not merely that one map is attractive. It is that a hierarchical model can stabilize small-area risks, describe spatially heterogeneous epidemic trajectories, evaluate out-of-sample performance and deliver probabilistic warnings. The historical paper’s original limitations also remain relevant: surveillance quality, population denominator uncertainty, prior sensitivity, spatial unit choice and unmeasured ecological drivers must be examined before operational decisions. The publication and the references in its final slides provide the scientific provenance for the formulas, interaction taxonomy and original figures in this chapter.

Day 3

Paper 2 · Clusters and Varying Effects

Agglomerative clustering, posterior risk boundaries and Bayesian regression with heterogeneous effects.

Day 3 - Paper 2: Clustering and Spatiotemporally Varying Coefficients

Bandung illustration. Use the district risk time series for agglomerative clustering. Display the candidate clusters on the same Bandung map and compare cluster-specific weather effects.

PPTX alignment: Paper 2 (slides 158–198). Why smoothing can erase discontinuities; two-stage spatiotemporal AHC; generation of candidate partitions; cluster fixed-effect, cluster random-effect and single-cluster models; priors and hyperpriors; identifiability and confounding; configuration selection; effects of weather variables; maps of high/low-risk clusters and the original study’s conclusions.

Two-stage logic. Stage 1 calculates a dissimilarity among space–time risk profiles and uses agglomerative clustering to produce nested candidate configurations. Stage 2 evaluates the candidate configuration using a Bayesian regression with covariate effects allowed to differ across risk regimes. A simplified expression is \(\log R_{it}=\alpha+\beta_{c(i,t)}^Tx_{it}+u_i+\gamma_t+\psi_{it}\). The clustering stage and the Bayesian stage answer different questions: the former defines candidate partitions; the latter assesses risk-factor relationships and estimates uncertainty conditional on the candidate structure.

Three model variants discussed in the slides. A cluster fixed-effect model assigns coefficients to selected clusters; a cluster random-effect model places a hierarchical distribution on those coefficients; the one-cluster model reduces to global coefficients. Differences between these specifications reflect assumptions about between-cluster borrowing of information. Intercepts and spatial effects must be constrained to avoid unidentifiable shifts, and correlated covariates can destabilize local coefficient estimation.

Bandung interpretation. Weather variables (including precipitation, humidity, sunshine and evaporation in the published application) can have different associations across risk clusters. Compare the district cluster map with the smoothed risk map and ask whether a boundary looks stable across months. The classroom clustering code is an accessible analogue of the paper’s procedure, not a numerical reproduction of the full spatiotemporal optimization.

Core sources for this session: agglomerative hierarchical clustering follows Ward (1963), while the applied Bayesian spatiotemporal clustering and varying-coefficient framework is connected to Jaya (2022).

Big idea.

Bayesian smoothing can make maps more stable, but it can also cause over-smoothing. This means that local boundaries between high-risk and low-risk areas may be blurred. The thesis addresses this using two stages:

1. Generate candidate clusters using agglomerative hierarchical clustering (AHC) based on space-time adjacency and risk similarity.
2. Select the optimal cluster configuration using Bayesian spatiotemporal regression with coefficients varying by cluster.

The big idea is that disease risk is not always smoothly changing across neighboring areas. Real outbreaks can have discontinuities: a district may be high-risk because of local housing, drainage, mobility, vegetation, water storage, or reporting practices, while a neighboring district remains lower-risk. If the model forces too much smoothness, it may miss exactly the local boundaries that public-health teams need.

Clustering adds a second lens:

• smoothing stabilizes noisy estimates,
• clustering preserves discontinuities,
• varying coefficients reveal whether risk factors have different effects in different local contexts.

Workshop scenario. Suppose rainfall increases risk in dense urban clusters but has a weaker effect in well-drained districts. A global coefficient for rainfall may average these patterns and produce an interpretation that is technically correct but operationally weak. A varying-coefficient model allows the rainfall effect to differ across clusters.

Rationale.

The rationale for hierarchical spatiotemporal clustering is that the risk surface may contain local discontinuities. A pure CAR/RW model is designed to smooth neighboring areas and adjacent periods. That is valuable when SIR is noisy, but it may blur real differences between local epidemiological regimes.

Clusters answer a different question:

\[ \text{Which area-time units behave similarly enough to be modeled together?} \]

Varying coefficients then answer:

\[ \text{Do risk factors have the same effect in every cluster, or do they differ by local context?} \]

This is especially important for dengue because rainfall, humidity, solar radiation, population density, mobility, and sanitation may not operate uniformly across the city.

Identifying cluster configurations using agglomerative clustering

Agglomerative hierarchical clustering is used as a structured way to build cluster configurations, but the resulting groups must still be interpreted in epidemiological context (Ward, 1963; Jaya, 2022).

Beginner intuition for agglomerative clustering.

Agglomerative clustering starts with many small groups and gradually merges them into larger groups. Imagine every district-month begins as its own cluster. The algorithm repeatedly joins the two most similar clusters until only a few groups remain.

The analyst then decides which number of clusters is useful. Too few clusters may hide important heterogeneity. Too many clusters may be difficult to interpret and may overfit noise.

A practical cluster should be:

• statistically supported;
• epidemiologically interpretable;
• geographically or temporally meaningful;
• useful for communication with decision-makers;
• stable enough to guide policy.

Space-time units are written as:

\[ \Omega = \{\omega_{11},\ldots,\omega_{it},\ldots,\omega_{nT}\} \]

We cluster based on:

\[ z_{it} = \log(\mathrm{SIR}_{it}) \]

Single-linkage distance:

\[ d_{ab} = \min\{ \lVert z_{it}-z_{js}\rVert: \omega_{it}\in C_a,\omega_{js}\in C_b \} \]

Centroid distance:

\[ d_{ab} = \lVert \bar{z}_{C_a} - \bar{z}_{C_b}\rVert \]

Ward distance:

\[ d_{ab} = \mathrm{ESS}(C_a\cup C_b) - \{\mathrm{ESS}(C_a)+\mathrm{ESS}(C_b)\} \]

Contiguity constraint:

\[ C_a \text{ may be merged with } C_b \quad \text{if they are spatial neighbors at the same time, or the same area in adjacent periods.} \]

The single cluster model with global fixed effects and varying coefficients

Beginner intuition for varying coefficients.

In ordinary regression, one coefficient applies everywhere:

\[ \log(\theta_{it}) = \alpha + \beta x_{it}. \]

This says a one-unit change in \(x\) has the same association with risk in every district and month. A varying-coefficient model relaxes that assumption:

\[ \log(\theta_{it}) = \alpha + \beta_{g(i,t)} x_{it}, \]

where \(g(i,t)\) is the cluster for area \(i\) and time \(t\). Now the covariate effect can differ by cluster. This is useful when the same environmental or social factor has different implications in different local contexts.

Start with the single-cluster model with global fixed effects:

\[ \eta_{it} = \alpha + u_i + v_i + \gamma_t + \phi_t + \delta_{it} + x_{it}^{\top}\beta \]

This model assumes one global relationship between covariates and log-relative risk. For example, the precipitation coefficient \(\beta_{\mathrm{rain}}\) is the same for all Bandung districts and all months.

The cluster-specific extension is:

\[ \eta_{it} = \alpha + u_i + v_i + \gamma_t + \phi_t + \delta_{it} + x_{it}^{\top}\beta + \sum_{c=1}^{K-1} I(\omega_{it}\in C_c) \left(a_c + x_{it}^{\top}b_c\right) \]

Interpretation:

• \(\beta\): global risk-factor effects.
• \(a_c\): intercept difference for cluster \(c\).
• \(b_c\): deviation in risk-factor effects for cluster \(c\).

If \(b_{c,\mathrm{rain}}>0\), then rainfall has a stronger effect in cluster \(c\) than in the global baseline. If \(b_{c,\mathrm{rain}}<0\), rainfall has a weaker effect in that cluster.

Prior and joint prior distributions of the single-cluster model.

For the single-cluster model with global fixed effects:

\[ \eta_{it} = \alpha + u_i+v_i+\gamma_t+\phi_t+\delta_{it}+x_{it}^{\top}\beta \]

a typical prior structure is:

\[ \alpha \sim N(0,\sigma_{\alpha}^{2}), \qquad \beta_k \sim N(0,\sigma_{\beta}^{2}) \]

\[ u\mid \tau_u \sim \mathrm{CAR}(\tau_u), \qquad v_i\mid \tau_v \sim N(0,\tau_v^{-1}) \]

\[ \gamma\mid\tau_\gamma \sim \mathrm{RW1}(\tau_\gamma) \quad \text{or} \quad \mathrm{RW2}(\tau_\gamma), \qquad \phi_t\mid\tau_\phi \sim N(0,\tau_\phi^{-1}) \]

\[ \delta\mid\tau_\delta \sim N(0,\tau_\delta^{-1}I) \quad \text{or} \quad \delta\sim \text{Knorr-Held interaction type I-IV} \]

The joint prior can be written compactly as:

\[ p(\Theta\mid\Lambda) = p(\alpha) p(\beta) p(u\mid\tau_u) p(v\mid\tau_v) p(\gamma\mid\tau_\gamma) p(\phi\mid\tau_\phi) p(\delta\mid\tau_\delta) \]

and the hyperprior is:

\[ p(\Lambda) = p(\tau_u)p(\tau_v)p(\tau_\gamma)p(\tau_\phi)p(\tau_\delta) \]

For the cluster-specific model, priors are added for the cluster intercept deviations and slope deviations:

\[ a_c\sim N(0,\tau_a^{-1}), \qquad b_{ck}\sim N(0,\tau_{b_k}^{-1}) \]

The complete posterior is then proportional to likelihood times all priors:

\[ p(\Theta,\Lambda\mid y) \propto \left[ \prod_{i,t} \frac{(E_{it}\theta_{it})^{Y_{it}}\exp(-E_{it}\theta_{it})}{Y_{it}!} \right] p(\Theta\mid\Lambda)p(\Lambda) \]

Selection of the optimal cluster configuration.

Agglomerative clustering can generate many possible cluster configurations. The optimal configuration should not be selected only because it looks attractive on a map. It should be evaluated by:

• DIC and WAIC for fit-complexity trade-off;
• CPO/PIT or MPL for predictive calibration;
• out-of-sample MAE and RMSE for forecasting;
• interpretability of cluster-specific coefficients;
• spatial and temporal coherence of the cluster map;
• feasibility of translating clusters into public-health action.

Candidate_K	Statistical_check	Substantive_check
2 clusters	Simple, low variance, may underfit local regimes	Useful for broad high/low risk policy separation
3-4 clusters	Balanced complexity for many city-level applications	Often interpretable as core, transition, and peripheral regimes
5-8 clusters	May improve fit but risks overfitting	Needs strong local explanation and clear intervention meaning
>8 clusters	Usually high complexity for a four-day workshop dataset	Use only if data volume and policy need justify detail

Measurement error, lags, and multicollinearity.

In disease surveillance, weather covariates often have the following problems:

\[ X_{it}^{obs} = X_{it}^{true} + \epsilon_{it} \]

Lagged effects can be written as:

\[ \eta_{it} = \cdots + \sum_{\ell=0}^{L}\beta_{\ell}X_{i,t-\ell} \]

Multicollinearity can be checked using VIF:

\[ \mathrm{VIF}_k = \frac{1}{1-R_k^2} \]

High VIF values indicate that covariate \(X_k\) is strongly explained by other covariates.

covariate	VIF
precip_lag1	5.109
temp	3.967
humidity_lag1	4.796
solar	3.333

As a practical rule, VIF values above 5 suggest moderate collinearity and values above 10 suggest serious collinearity. In disease modeling, this should trigger a substantive discussion: should the variables be combined, lagged differently, standardized, or selected based on biological plausibility?

A short clustering example. The next practical section shows how candidate patterns are constructed and mapped; it also contains a more complete optional analysis with covariates.

avg <- aggregate(cases ~ district_id, data=panel, FUN=mean)
groups <- cutree(hclust(dist(avg$cases), method="ward.D2"), k=3)
districts$demo_cluster <- groups[match(districts$district_id,avg$district_id)]
plot(sf::st_geometry(districts), col=c("#ebd3b1","#c98f5d","#8b4a2e")[districts$demo_cluster])

Case study: dengue disease in Bandung - classroom clustering

The quick example below clusters districts using summary risk and covariate features. For a thesis-level application, space-time contiguity constraints should be activated.

The classroom version deliberately uses a compact feature matrix so that participants can focus on interpretation before moving to constrained AHC. Each district is summarized by:

• average log-SIR,
• volatility of log-SIR,
• average lagged precipitation,
• average temperature,
• average lagged humidity,
• average solar radiation.

The resulting clusters are not final epidemiological truth. They are a teaching device for understanding how local risk profiles can be grouped before a more rigorous constrained Bayesian model is fitted.

This cluster map groups Bandung districts with similar simulated risk and weather profiles. It should be interpreted as an exploratory classification, not a final causal conclusion. If neighboring districts fall into different clusters, participants should ask whether the difference is epidemiologically plausible or whether a contiguity-constrained clustering method is needed.

Bayesian-regression teaching analogue with cluster effects.

term	Estimate	Std. Error	z value	Pr(>|z|)
cluster2	-0.5941	0.0274	-21.6480	0.0000
cluster3	-0.3514	0.0295	-11.9236	0.0000
cluster4	-0.8465	0.0381	-22.2231	0.0000
scale(precip_lag1)	0.2004	0.0357	5.6189	0.0000
scale(temp)	0.2180	0.0299	7.2940	0.0000
scale(humidity_lag1)	0.1955	0.0322	6.0715	0.0000
scale(solar)	-0.0777	0.0276	-2.8178	0.0048
cluster2:scale(precip_lag1)	0.0160	0.0441	0.3631	0.7165
cluster3:scale(precip_lag1)	0.0207	0.0440	0.4695	0.6387
cluster4:scale(precip_lag1)	-0.0127	0.0671	-0.1892	0.8499
cluster2:scale(temp)	-0.0379	0.0375	-1.0085	0.3132
cluster3:scale(temp)	-0.1129	0.0410	-2.7536	0.0059
cluster4:scale(temp)	-0.0580	0.0569	-1.0196	0.3079
cluster2:scale(humidity_lag1)	-0.0482	0.0438	-1.1012	0.2708
cluster3:scale(humidity_lag1)	-0.0656	0.0428	-1.5353	0.1247
cluster4:scale(humidity_lag1)	-0.0785	0.0598	-1.3118	0.1896
cluster2:scale(solar)	-0.0015	0.0355	-0.0431	0.9656
cluster3:scale(solar)	0.0278	0.0371	0.7481	0.4544
cluster4:scale(solar)	-0.0478	0.0514	-0.9302	0.3523

This coefficient table shows whether weather effects differ across Bandung risk clusters. The key teaching point is not only statistical significance; it is whether the direction and magnitude of the effect changes by local context. For example, precipitation may have a stronger simulated association in a humid northern cluster than in a dense central cluster.

How to read the interaction terms:

• A main covariate coefficient is the effect for the reference cluster.
• A cluster:covariate term is the difference between a given cluster and the reference cluster.
• The cluster-specific effect is obtained by adding the main effect and the interaction effect.
• If cluster-specific effects differ meaningfully, a single global model may hide important local mechanisms.

This map is the cluster-specific model’s fitted relative risk. Compared with raw SIR, it is smoother; compared with the pure model, it can preserve sharper differences among local risk regimes. For public health, this helps distinguish places that need the same intervention from places that need different intervention logic.

AHC template with space-time contiguity.

The full paper develops spatiotemporally constrained clustering. In this four-hour teaching module, we implement Ward’s AHC directly with hclust() and visualize the selected groups on the Bandung map. The additional constrained search algorithm is omitted to keep all practical R code short and executable.

INLA template for spatiotemporally varying coefficients

##                               mean           sd    0.025quant      0.5quant
## (Intercept)          -3.4998404425 0.5713455829 -4.6200710623 -3.4996189280
## cluster2             -0.4724133114 0.1914006087 -0.8472008780 -0.4735713510
## cluster3             -0.3110401608 0.2030276472 -0.7107894042 -0.3117012617
## cluster4             -0.6607261822 0.2359651711 -1.1222061936 -0.6623712500
## precip_lag1           0.0047421760 0.0004047070  0.0039444885  0.0047437018
## temp                  0.0865842489 0.0206046881  0.0463768208  0.0865459043
## humidity_lag1         0.0211076806 0.0019989191  0.0171780440  0.0211113212
## solar                -0.0050516561 0.0006459070 -0.0063164679 -0.0050520644
## cluster2:precip_lag1 -0.0002665459 0.0003181318 -0.0008903635 -0.0002665484
## cluster3:precip_lag1  0.0003385663 0.0003260789 -0.0003008306  0.0003385622
## cluster4:precip_lag1 -0.0001890195 0.0004772467 -0.0011248364 -0.0001890253
##                         0.975quant          mode                 kld
## (Intercept)          -2.3811736783 -3.4996919967 0.00000000196696289
## cluster2             -0.0912471216 -0.4736023166 0.00000006867123823
## cluster3              0.0924863248 -0.3116865381 0.00000008641718131
## cluster4             -0.1899645925 -0.6623645475 0.00000007535762404
## precip_lag1           0.0055313045  0.0047437457 0.00000000382777724
## temp                  0.1270308970  0.0865502993 0.00000000875312129
## humidity_lag1         0.0250167717  0.0211113913 0.00000000095806774
## solar                -0.0037844460 -0.0050520491 0.00000000084464811
## cluster2:precip_lag1  0.0003572864 -0.0002665484 0.00000000005516838
## cluster3:precip_lag1  0.0009779864  0.0003385622 0.00000000005519622
## cluster4:precip_lag1  0.0007468305 -0.0001890253 0.00000000005523837

Day 3 exercise. Change the number of clusters from 2 to 8. Compare the interpretation of precipitation, temperature, humidity, and solar radiation effects. Discuss when a global model becomes misleading because local effects are very different.

Paper 2 in depth: the complete two-stage methodology

The accompanying presentation (Identifying Spatiotemporal Clusters by Means of Agglomerative Hierarchical Clustering and Bayesian Regression Analysis) follows Jaya and Folmer (2021, Geographical Analysis, DOI: 10.1111/gean.12264). It presents a two-stage procedure, rather than a standard one-step clustering analysis. This section consolidates the equations, assumptions, prior specification, algorithm and empirical discussion from the 41 presentation slides. Scientific figures are reproduced individually where they clarify the results; the PowerPoint slide layouts are not reproduced.

Study units and standardized risk (slides 1–5). Let \(\mathcal A_{it}\) be the \(i\)-th district in month \(t\), with \(n\) districts and \(T\) time periods. The observed number of cases is \(Y_{it}\), the expected number is \(E_{it}\), and the unknown relative risk is \(\theta_{it}\). The starting point is

\[Y_{it}\mid\theta_{it}\sim\operatorname{Poisson}(E_{it}\theta_{it}),\qquad \widehat\theta_{it}^{\rm SIR}=\frac{Y_{it}}{E_{it}},\qquad S_{it}=\log(\widehat\theta_{it}^{\rm SIR}).\]

For example, if a Bandung district has 24 dengue cases but only 12 expected cases, its SIR equals \(24/12=2.0\) and its log-SIR is \(\log 2\approx0.693\). If the neighboring district has 6 cases and 12 expected, its SIR is \(0.5\), with log-SIR \(\approx-0.693\). The absolute log-SIR difference is therefore \(1.386\). For observed zero counts, directly taking the logarithm is undefined; the classroom examples use a small positive adjustment, whereas the likelihood model uses the raw counts and expected values without that adjustment.

Stage 1: adjacency-constrained agglomeration (slides 6–11). Start with \(nT\) singleton clusters. At each iteration, merge the most similar eligible pair. Eligibility is determined by actual space–time contiguity: two units may be linked when they are neighboring districts in the same month, or the same district in consecutive months. Consequently, two nearby log-SIR values cannot merge directly if their units are not spatiotemporally connected. The distinction is crucial: ordinary hclust() on a risk vector does not enforce this rule.

\[ a\bigl((i,t),(j,s)\bigr)= \begin{cases} 1,& i\sim j\ \text{and}\ t=s,\\ 1,& i=j\ \text{and}\ |t-s|=1,\\ 0,&\text{otherwise}. \end{cases} \]

For two clusters \(C_k\) and \(C_l\), define \(d(C_k,C_l)\) using a chosen linkage criterion. The eligible dissimilarity is

\[ D_{kl}=\begin{cases} d(C_k,C_l),&\exists\,u\in C_k,\ v\in C_l:\ a(u,v)=1,\\ \infty,&\text{otherwise}. \end{cases} \qquad (k^*,l^*)=\arg\min_{k\ne l}D_{kl}. \]

The paper describes single linkage (minimum pairwise distance), centroid linkage (difference in cluster means), and Ward linkage (minimum increase in within-cluster error sum of squares). For scalar log-SIR values, Ward’s merge cost has a particularly transparent expression:

\[ \Delta_{kl}^{\rm Ward} =\frac{n_kn_l}{n_k+n_l}(\bar S_k-\bar S_l)^2. \]

Calculate one Ward merge by hand. Suppose \(C_k\) has two values \((0.1,0.3)\) and \(C_l\) has one value \((0.5)\). Their means are \(0.2\) and \(0.5\). If adjacent, the Ward cost is \((2\times1)/(2+1)\,(0.2-0.5)^2=0.06\). A nonadjacent pair with a smaller numeric cost would still be ineligible. This is why the method preserves geographically meaningful local boundaries.

cluster_A <- c(0.1, 0.3)
cluster_B <- 0.5
ward_cost <- length(cluster_A) * length(cluster_B) /
  (length(cluster_A) + length(cluster_B)) *
  (mean(cluster_A) - mean(cluster_B))^2
round(ward_cost, 4)

## [1] 0.06

Complete stage-1 algorithm in words. (1) Choose the dissimilarity and linkage. (2) Put every area–month into its own cluster. (3) Construct the dissimilarity matrix, assigning infinity to pairs that have no space–time connection. (4) Merge the eligible pair with minimum dissimilarity. (5) Update the cluster memberships, spatial–temporal connectivity, and merge costs. (6) Save the resulting candidate partition and repeat until one cluster remains. For a panel of 30 Bandung districts over 12 months there are \(30\times12=360\) singleton units and, in principle, candidate partitions with one to 360 clusters.

Stage 2: cluster-specific Bayesian regression (slides 12–15). The second stage evaluates the candidate partitions using a Poisson log-linear model, with structured and unstructured effects and weather covariates. Write the pure spatiotemporal component as

\[ \eta_{it}^{\rm pure}=\alpha+\omega_i+\nu_i+\phi_t+\gamma_t+\delta_{it}, \]

where \(\omega_i\) and \(\phi_t\) are structured spatial and temporal terms, \(\nu_i\) and \(\gamma_t\) are unstructured terms, and \(\delta_{it}\) is the interaction. Let \(x_{it}\) be the vector of weather covariates and let \(c(i,t)\) identify the cluster. The single-cluster baseline (Eq. 7 in the paper) is

\[M_0:\quad\eta_{it}=\eta_{it}^{\rm pure}+x_{it}^{\mathsf T}\boldsymbol\beta.\]

For a partition with \(g\) clusters, the fixed cluster effects (Eq. 5) become

\[ M_{\rm fixed}:\quad\eta_{it}=\eta_{it}^{\rm pure}+x_{it}^{\mathsf T}\boldsymbol\beta +\sum_{l=1}^{g-1}\mathbf1\{c(i,t)=l\}\left(o_l+x_{it}^{\mathsf T}\boldsymbol\zeta_l\right). \]

The last cluster is the reference to avoid redundant dummy parameters. Its coefficients are \(\boldsymbol\beta\); for cluster \(l\) they are \(\boldsymbol\beta_l=\boldsymbol\beta+\boldsymbol\zeta_l\). Under the random cluster effects (Eq. 6), the cluster deviations \(o_{c(i,t)}\) and \(\boldsymbol\zeta_{c(i,t)}\) receive hierarchical distributions and borrow information across clusters:

\[ M_{\rm random}:\quad\eta_{it}=\eta_{it}^{\rm pure} +x_{it}^{\mathsf T}\boldsymbol\beta +o_{c(i,t)}+x_{it}^{\mathsf T}\boldsymbol\zeta_{c(i,t)}. \]

For illustration, suppose the global precipitation coefficient is \(\beta_{\rm rain}=0.15\) and the deviation for one high-risk cluster is \(\zeta_{\rm rain}=0.20\). The within-cluster coefficient is \(0.35\). A one-unit standardized precipitation increase corresponds to a relative-risk multiplier \(\exp(0.35)\approx1.42\) in that cluster, compared with \(\exp(0.15)\approx1.16\) in the reference cluster. These are illustrative numbers, not estimated coefficients from the published Bandung case study.

global_rain <- 0.15
cluster_deviation <- 0.20
c(reference_RR = exp(global_rain),
  cluster_RR = exp(global_rain + cluster_deviation))

## reference_RR   cluster_RR 
##     1.161834     1.419068

Priors and hyperpriors (slides 16–21). The model is hierarchical: weakly informative Gaussian distributions can describe the intercept and global coefficients, CAR structures can represent spatial autocorrelation, and RW1 can describe the temporal trend. The spatial CAR precision has diagonal entries proportional to the number of neighbors and off-diagonal entries \(-1\) for adjacent units, with zero otherwise. This can be written as \(Q_{\omega}=\tau_{\omega}(D-W)\), subject to identification constraints. Unstructured space and time terms receive independent Gaussian priors; interaction types I–IV depend on whether the spatial and temporal parts are structured. The presentation discusses inverse-Gamma hyperpriors for temporal variances, half-Cauchy distributions for spatial scale parameters, and penalized-complexity (PC) priors for cluster-effect scales. Different prior choices have different regularization implications and must be checked in a sensitivity analysis.

Prior shrinkage example. A PC prior can express \(P(\sigma_{\rm cluster}>U)=\alpha_0\). Under an exponential prior on the standard-deviation scale, \(p(\sigma)=\lambda_0\exp(-\lambda_0\sigma)\), the corresponding rate is \(\lambda_0=-\log(\alpha_0)/U\). For \(U=1\) and \(\alpha_0=.05\), \(\lambda_0\approx2.996\). This favors small cluster deviations unless the observed data support larger effects.

U <- 1
alpha0 <- 0.05
lambda0 <- -log(alpha0) / U
c(pc_prior_rate = lambda0,
  check_tail_probability = exp(-lambda0 * U))

##          pc_prior_rate check_tail_probability 
##               2.995732               0.050000

Confounding and identifiability (slides 22–24). Cluster intercepts, global intercepts and spatial effects can become redundant. A cluster covariate deviation can also compete with a smooth latent effect that explains the same pattern. Suitable reference coding, constraints, centering and prior regularization are needed to obtain interpretable components. However, additional hard constraints may change estimation and fit, so they should be justified and assessed rather than imposed mechanically. The article reports a DIC of 2476.26 for its preferred model without an additional linear constraint versus 2486.28 with the constraint in that comparison; the lower figure favors the former within that specific candidate comparison.

The actual identifying constraints. In intrinsic CAR, RW and interaction constructions, a constant vector may lie in the null space of the precision matrix. The article therefore discusses sum-to-zero restrictions on structured effects and, where appropriate, cluster deviations:

\[ \sum_{i=1}^{n}\omega_i=0,\qquad \sum_{t=1}^{T}\phi_t=0,\qquad \sum_{i=1}^{n}\delta_{it}=0\ \text{for each }t,\qquad \sum_{t=1}^{T}\delta_{it}=0\ \text{for each }i, \]

\[\sum_{l=1}^{g}\zeta_{lk}=0\quad\text{for covariate }k.\]

These restrictions address identifiability with an intercept; they do not automatically solve confounding between a covariate and an unobserved spatial effect. To illustrate the latter, write \(\log\boldsymbol\theta=\alpha\boldsymbol1+X\boldsymbol\beta+\boldsymbol\varrho\). If \(X\) has full column rank, decompose the latent effect using \(P_X=X(X^{\mathsf T}X)^{-1}X^{\mathsf T}\) and its complement \((I-P_X)\). Restricting \(X^{\mathsf T}\boldsymbol\varrho=0\) moves shared variation into the fixed-effect part, but changes the model’s meaning and potentially its fit. The paper examines such additional constraints empirically rather than treating them as automatically preferable.

Multicollinearity and ridge shrinkage. When weather measurements and lags strongly correlate, their regression coefficients can be unstable. The article relates ridge regularization to a Gaussian-prior Bayesian MAP estimate. A local weighted approximation to a Poisson regression gives

\[ \widehat{\boldsymbol\beta}_{\rm ridge}(\lambda) =\left(X^{\mathsf T}WX+\lambda I\right)^{-1}X^{\mathsf T}W\boldsymbol z, \]

where \(W\) is the diagonal working-weight matrix and \(\boldsymbol z\) the working response; in ordinary least squares set \(W=I\) and \(\boldsymbol z=y\). A normal prior with precision proportional to \(\lambda\) shrinks slopes toward zero. In an intercept model, one ordinarily leaves the intercept unpenalized. The following tiny numerical illustration demonstrates the stability–shrinkage trade-off without requiring an optimization package.

set.seed(20261126)
x1 <- seq(-1, 1, length.out = 20)
x2 <- x1 + rnorm(20, sd = 0.02)  # Highly correlated predictors.
y <- 1 + 0.8 * x1 + rnorm(20, sd = 0.2)
X <- cbind(x1, x2)
for (penalty in c(0.01, 0.5, 2)) {
  estimate <- solve(crossprod(X) + diag(penalty, 2),
                    crossprod(X, y - mean(y)))
  print(round(c(lambda = penalty, beta1 = estimate[1],
                beta2 = estimate[2]), 3))
}

## lambda  beta1  beta2 
##  0.010  0.374  0.487 
## lambda  beta1  beta2 
##  0.500  0.412  0.422 
## lambda  beta1  beta2 
##  2.000  0.377  0.383

Model-selection criteria. The presentation compares candidate configurations using DIC, WAIC, CPO, PIT, MPL and posterior predictive fit. DIC is often written as

\[\operatorname{DIC}=\overline{D(\theta)}+p_D,\qquad D(\theta)=-2\log p(y\mid\theta)+C,\]

with \(p_D\) the effective complexity penalty. The paper’s cluster-selection scree plot uses a stabilization/elbow heuristic, not a universal hard threshold. In the published application, 14 clusters were selected from the 360 candidate configurations using the DIC scree-plot approach on a subset of 150 configurations. Changing \(K\) can alter coefficient estimates and the inferred boundaries between high-risk and low-risk regimes.

How CPO, PIT and MPL are calculated. For observation \(y_j\), let \(y_{-j}\) denote the data excluding that observation. The leave-one-out conditional predictive ordinate and marginal predictive likelihood are

\[\operatorname{CPO}_j=p(y_j\mid y_{-j}),\qquad \operatorname{MPL}=\sum_{j=1}^{nT}\log(\operatorname{CPO}_j).\]

The probability integral transform (PIT) evaluates the predictive cumulative distribution at the observed outcome, with randomization needed for a formally uniform PIT under discrete counts. A nearly flat PIT histogram is compatible with good probabilistic calibration; a U-shaped pattern may indicate underdispersion, while a peaked pattern may indicate overdispersion. The article also uses a pseudo-Bayes-factor approximation between candidate models,

\[\operatorname{PBF}_{1,2}=\exp\{\operatorname{MPL}(M_1)-\operatorname{MPL}(M_2)\}.\]

For example, if the MPL values are \(-120\) and \(-117\), then \(\operatorname{PBF}_{1,2}=e^{-3}\approx0.0498\), indicating a relative preference for \(M_2\) under this predictive comparison.

MPL_M1 <- -120
MPL_M2 <- -117
exp(MPL_M1 - MPL_M2)

## [1] 0.04978707

Classroom illustration of a model-selection curve. The following short R code compares candidate district-level groupings by AIC, which is fast and available in base R. AIC is only a pedagogical analogue of the article’s Bayesian DIC procedure. It must not be described as the published 14-cluster selection.

candidate_K <- 2:8
candidate_AIC <- vapply(candidate_K, function(k) {
  district_group <- cutree(hc, k = k)
  dat <- as.data.frame(panel)
  dat$group <- factor(district_group[dat$district_id])
  fit <- glm(cases ~ offset(log(expected)) + group +
               precip_lag1 + temp + humidity_lag1 + solar +
               factor(month_no), family = poisson(), data = dat)
  AIC(fit)
}, numeric(1))
plot(candidate_K, candidate_AIC, type = "b", pch = 19,
     col = "#A95F36", lwd = 2, xlab = "Number of district groups",
     ylab = "AIC (teaching analogue, not published DIC)",
     main = "Candidate clustering: fit versus complexity")

A faithful classroom implementation of the spatiotemporal adjacency rule. The earlier hclust() illustration is useful for beginners but clusters district summaries without constraining geography. The next executable example implements the paper’s core space–time eligibility rule explicitly on all 30 actual Bandung kecamatan over a simulated twelve-month cycle. It uses the queen-adjacency matrix already calculated from the actual supplied Bandung polygons, plus consecutive observations of the same district. We keep the Ward criterion and save the candidate configurations with 4, 8, 14 and 20 clusters. This is an instructional implementation of Stage 1, not the authors’ original optimized program.

# District-by-month observations: 30 areas x 12 months = 360 units.
ahc_data <- as.data.frame(panel)
ahc_data <- ahc_data[ahc_data$t <= 12, ]
ahc_data <- ahc_data[order(ahc_data$t, ahc_data$district_id), ]
rownames(ahc_data) <- NULL
n_units <- nrow(ahc_data)
stopifnot(n_units == 360)

# Eligible links: geographic neighbors in the same month OR
# identical district in consecutive months.
same_month <- outer(ahc_data$t, ahc_data$t, "==")
same_district <- outer(ahc_data$district_id,
                       ahc_data$district_id, "==")
next_month <- abs(outer(ahc_data$t, ahc_data$t, "-")) == 1
eligible <- (W_binary[ahc_data$district_id,
                      ahc_data$district_id] > 0 & same_month) |
            (same_district & next_month)
diag(eligible) <- FALSE

# Every unit starts in its own cluster. A Ward merge updates
# cluster sizes, cluster means and the eligible-neighbor graph.
members <- as.list(seq_len(n_units))
cluster_size <- rep(1, n_units)
cluster_mean <- log((ahc_data$cases + 0.5) /
                    (ahc_data$expected + 0.5))
save_at <- c(20, 14, 8, 4)
saved_partitions <- list()
while (length(members) > min(save_at)) {
  pair_index <- which(eligible & upper.tri(eligible), arr.ind = TRUE)
  if (nrow(pair_index) == 0) {
    stop("The Bandung space-time adjacency graph is disconnected.")
  }
  a <- pair_index[, 1]
  b <- pair_index[, 2]
  merge_cost <- cluster_size[a] * cluster_size[b] /
    (cluster_size[a] + cluster_size[b]) *
    (cluster_mean[a] - cluster_mean[b])^2
  winner <- which.min(merge_cost)
  a <- a[winner]
  b <- b[winner]  # Always b > a from upper.tri().

  new_n <- cluster_size[a] + cluster_size[b]
  cluster_mean[a] <- (cluster_size[a] * cluster_mean[a] +
                      cluster_size[b] * cluster_mean[b]) / new_n
  cluster_size[a] <- new_n
  members[[a]] <- c(members[[a]], members[[b]])
  eligible[a, ] <- eligible[a, ] | eligible[b, ]
  eligible[, a] <- eligible[a, ]

  members <- members[-b]
  cluster_size <- cluster_size[-b]
  cluster_mean <- cluster_mean[-b]
  eligible <- eligible[-b, -b, drop = FALSE]
  diag(eligible) <- FALSE
  if (length(members) %in% save_at) {
    labels <- integer(n_units)
    for (j in seq_along(members)) labels[members[[j]]] <- j
    saved_partitions[[as.character(length(members))]] <- labels
  }
}
stopifnot(all(as.character(save_at) %in% names(saved_partitions)))
ahc_data$cluster14 <- factor(saved_partitions[["14"]])
table(ahc_data$cluster14)[1:min(14, nlevels(ahc_data$cluster14))]

## 
##  1  2  3  4  5  6  7  8  9 10 11 12 13 14 
## 17 87 35  3 27  2 69  1 50 42  7  2 11  7

The mapping below makes the space–time continuity visible. The same district can change cluster membership between months because it is the district-month that is clustered. A cluster may extend through several districts and months. The colour identifies an algorithmic group, not automatically a high-risk or low-risk label.

display_months <- ahc_data[ahc_data$t %in% c(1, 6, 12),
                           c("district_id", "t", "cluster14")]
ahc_map <- dplyr::left_join(
  districts |> dplyr::select(district_id),
  display_months, by = "district_id")
ggplot2::ggplot(ahc_map) +
  ggplot2::geom_sf(ggplot2::aes(fill = cluster14),
                   colour = "white", linewidth = .2) +
  ggplot2::facet_wrap(~t, nrow = 1,
                      labeller = ggplot2::label_both) +
  ggplot2::scale_fill_viridis_d(option = "C", name = "Cluster") +
  ggplot2::labs(title = "Space-time constrained AHC: 14 candidate clusters",
                subtitle = "Actual Bandung polygons; simulated monthly SIR") +
  ggplot2::theme_minimal(base_size = 11)

A Bayesian DIC exercise for Stage 2. We now evaluate the same constrained-AHC partitions with a small Poisson Bayesian model in R-INLA. It uses an expected-count offset, global weather covariates, a candidate cluster factor, and spatial/temporal random effects. Keeping the specification the same across candidate \(K\) allows a meaningful comparison of the DIC values within this teaching example. The full published model has richer structured and unstructured components, cluster-specific slopes, priors and identification constraints; do not equate its reported 14-cluster solution with whatever minimizes the simplified classroom DIC.

fit_cluster_dic <- function(k) {
  dat <- ahc_data
  dat$candidate_group <- factor(saved_partitions[[as.character(k)]])
  fit <- INLA::inla(
    cases ~ 1 + candidate_group + precip_lag1 + temp +
      humidity_lag1 + solar +
      f(district_id, model = "iid") + f(t, model = "rw1"),
    family = "poisson", data = dat, E = expected,
    control.predictor = list(compute = TRUE),
    control.compute = list(dic = TRUE, waic = TRUE))
  c(DIC = fit$dic$dic, WAIC = fit$waic$waic)
}
dic_results <- t(vapply(c(4, 8, 14, 20), fit_cluster_dic,
                        numeric(2)))
dic_results <- data.frame(K = c(4, 8, 14, 20), dic_results)
print(round(dic_results, 2), row.names = FALSE)

##   K          DIC     WAIC
##   4      1943.15  1951.07
##   8      1898.55  1982.48
##  14   -230554.76  5840.85
##  20 -15184746.04 22073.11

plot(dic_results$K, dic_results$DIC, type = "b", pch = 19,
     col = "#A95F36", lwd = 2, xlab = "Number of space-time clusters",
     ylab = "DIC", main = "R-INLA: teaching example of Stage 2")

Why does this implementation matter? It makes the difference between the two stages observable in R: Stage 1 merges only eligible space–time neighbors according to log-SIR, while Stage 2 compares the statistical fit of their resulting partitions. Changing only the number of clusters does not recover the original paper’s posterior results, and DIC should not be used as the sole criterion for deciding outbreak-control boundaries. Nevertheless, participants can now trace each important mathematical object in the paper to a numerical calculation and a Bandung map.

Paper 2: published Bandung findings and maps

Actual application (slides 25–28). The study used the real administrative structure of 30 Bandung districts and dengue surveillance data from January 2012 to December 2016. The presentation reports 21,483 cases over the five-year period, with observed monthly district incidences ranging from 2 to 230. The weather risk factors were temperature, precipitation, sunshine duration, evaporation and humidity; delayed effects at zero, one and two months were considered. The paper uses the mean annual monthly pattern for its clustering illustration to reduce the computational burden and emphasize recurring seasonality. The e-book’s runnable example, by contrast, uses simulated cases and weather on the supplied Bandung boundaries, so its numerical outcomes should not be confused with these published estimates.

Temporal weather patterns (slide 29). Averages are not a sufficient summary when weather mechanisms have distinct seasonal cycles. The participants can plot the simulated Bandung precipitation, temperature, humidity and solar proxies by month, compare their peaks and discuss why dengue incidence may respond with a lag. The original slide also includes evaporation, which is discussed conceptually; the teaching dataset does not contain a separately validated evaporation measurement.

seasonal_weather <- panel |>
  dplyr::group_by(month_no) |>
  dplyr::summarise(Precipitation = mean(precip), Temperature = mean(temp),
                   Humidity = mean(humidity), Solar = mean(solar),
                   .groups = "drop")
weather_long <- tidyr::pivot_longer(seasonal_weather, -month_no,
                                   names_to = "Variable", values_to = "Value")
ggplot2::ggplot(weather_long, ggplot2::aes(month_no, Value)) +
  ggplot2::geom_line(colour = "#9A532B", linewidth = .8) +
  ggplot2::geom_point(colour = "#9A532B") +
  ggplot2::facet_wrap(~Variable, scales = "free_y") +
  ggplot2::scale_x_continuous(breaks = c(1, 3, 5, 7, 9, 11)) +
  ggplot2::labs(title = "Simulated Bandung weather: average annual cycle",
                x = "Month", y = "Simulated value") +
  ggplot2::theme_minimal(base_size = 11)

Published model assessment (slides 30–33). The article shows a DIC scree plot across candidate cluster configurations, an observed-versus-predicted plot, a probability integral transform (PIT) assessment and posterior coefficient estimates. A small PIT or predictive-check exercise is useful because a lower DIC does not alone demonstrate calibrated predictions. The two extracted illustrations below are original scientific figures only, not screenshots of full PowerPoint slides:

Paper 2 original observed-versus-predicted dengue cases: scientific figure

Paper 2 original PIT histogram: scientific figure

Spatiotemporally varying weather coefficients (slides 34–37). The paper examines how weather relationships change across the selected spatiotemporal clusters. The sign and magnitude of a coefficient must be interpreted relative to its covariate definition, lag and reference cluster. The publication’s multi-panel coefficient map gives an overview for temperature, precipitation, sunshine, evaporation and humidity across the annual cycle. These are published estimates and should be distinguished from the simulated cluster-specific model plots earlier in Day 3.

Paper 2 original varying-coefficient maps across weather factors and months

Risk and cluster maps (slides 38–39). A map of posterior relative risk shows the strength of predicted risk, while a map of high- and low-risk clusters displays the partition. These are related but distinct outputs. In the paper, high-risk clusters concentrate in northern, central and southern Bandung, become more prevalent from January through April, peak around May, and then decline in the dry-season months. The full annual progression matters: a district that is high-risk in May need not be high-risk in September.

Paper 2 original twelve-month posterior relative-risk maps

Posterior classification is different from a simple mean-risk threshold. The second map uses posterior exceedance probabilities to label districts as high-risk, low-risk or statistically inconclusive. For a threshold \(r_0\), compute \(P(\theta_{it}>r_0\mid y)\); an exceedance probability near 0.5 expresses substantial uncertainty even if the posterior mean exceeds \(r_0\). In the following original figure, grey areas are not confidently classified, while red and cyan indicate high- and low-risk districts, respectively.

Paper 2 original monthly posterior high-low-insignificant classification

Paper 2 original monthly high-risk and low-risk cluster maps

Conclusions (slides 40–41). The study proposes a practical compromise between smoothing noisy small-area rates and identifying local discontinuities. Stage 1 produces spatiotemporally contiguous candidate clusters based on log-SIR; stage 2 selects a useful configuration through Bayesian regression with global and cluster-specific effects. Risk is associated with weather, but the reported patterns are not automatically causal effects. The key scientific contribution is the joint treatment of local cluster boundaries, changing weather coefficients and posterior uncertainty. For operational surveillance, future applications should also assess the stability of selected clusters to prior choice, lag selection, spatial resolution and time aggregation.

Day 4 - Paper 3: Spatiotemporal Forecasting and Confounding

Bandung illustration. Reserve the last six months for forecasting. Compare covariate-only, pure and full models on unseen monthly cases, then inspect their spatial forecast patterns and MAE/RMSE.

PPTX alignment: Paper 3 (slides 199–227). Explanation versus prediction; pure versus covariate forecasting; omitted factors, coefficient instability and future covariate availability; Bayesian GMRF assumptions, priors and forecast design; Monte Carlo experimental design; continuous, count and binary outcomes; M1/M2/M3 comparisons and interpretation of the reported results.

Prediction is different from coefficient interpretation. A covariate can explain historical variation and still reduce future prediction accuracy because it is confounded with a latent spatial effect, changes over time, or is not observed for the forecast horizon. Forecasting must condition on the information actually available at the prediction date. In the instructional exercise, future covariates are scenario inputs known because they were simulated; in practice their forecasting uncertainty needs separate treatment.

The model comparison in the PPTX. M1 uses a potentially confounded covariate without the full latent structure; M2 uses pure spatiotemporal random effects; M3 includes both covariates and latent effects. The original study compares these structures under Gaussian, Poisson and Bernoulli response families, various levels of confounding and repeated simulations. Use RMSE/MSPE for continuous outcomes, count-appropriate scores for Poisson counts, and Brier or log score for binary prediction. One small teaching simulation cannot establish the paper’s general ranking.

Reported methodological takeaway. In the study’s experiments the covariate-only M1 model is generally outperformed; for the examined count and binary settings the pure model commonly performs better than the full model, while the Gaussian result depends on the degree of confounding. The slides also discuss how temporal dimension and model complexity affect predictive accuracy. These conclusions are conditional on the paper’s design, not universal forecasting laws.

Core sources for this session: the distinction between explanation and prediction follows Shmueli (2010), while forecast evaluation follows Hyndman and Athanasopoulos (2021). Applied disease-risk forecasting is linked to Jaya (2022).

Explanation vs prediction

This section follows the distinction between explanatory modeling and predictive modeling: a model that explains associations is not automatically the best forecasting model (Shmueli, 2010).

A model can be useful in different ways. Sometimes we want to understand causes or associations. Sometimes we want the most accurate warning signal. These goals overlap, but they are not identical.

A weather analogy helps:

• An explanatory model asks why rain occurs: humidity, pressure, wind, geography.
• A predictive model asks whether it will rain tomorrow and how accurate the warning is.

In disease surveillance, an explanatory model helps understand risk factors, while a predictive model helps allocate resources before cases increase.

Explanation and prediction are related but not identical.

An explanatory model asks:

Which risk factors are associated with disease risk, and how should we interpret their effects?

A predictive model asks:

How accurately can we forecast future cases or future relative risk?

In applied disease surveillance, both objectives matter. A health department wants to know why risk is high, but it also needs a forecast that helps allocate resources before cases rise.

Dimension	Explanation	Prediction
Main objective	Interpret mechanisms and associations	Forecast future observations or risk
Preferred output	Coefficients, credible intervals, effect signs	Forecast error, calibration, forecast maps
Model choice	Theoretically meaningful covariates and priors	Variables and structures that improve out-of-sample performance
Risk	A highly predictive model may be hard to interpret	A highly interpretable model may forecast poorly
Workshop balance	Use covariates and random effects for interpretation	Evaluate forecasts using held-out Bandung months

Simulation evidence using held-out Bandung months

The Bandung data are simulated, so we know that weather-like variables and spatial structure were used to generate risk. We can use this simulated setting to show the difference between explanatory and predictive performance.

The exercise below trains models on the first 30 months and evaluates them on the last 6 months.

Model	MAE	RMSE	Mean_predicted_cases	Mean_observed_cases
Explanation model with weather-like covariates	2.955	4.173	9.631	10.806
Prediction model with district, seasonality, and smooth time	2.996	4.282	9.443	10.806

The table is simulation evidence, not a universal theorem. It shows how model performance can differ when the goal changes. The model with weather-like covariates is easier to interpret epidemiologically. The prediction model focuses more directly on district baseline risk, seasonality, and temporal shape.

test_eval |>
  group_by(date) |>
  summarise(
    observed = sum(cases),
    explanation = sum(pred_explanation),
    prediction = sum(pred_prediction),
    .groups = "drop"
  ) |>
  tidyr::pivot_longer(-date, names_to = "series", values_to = "cases") |>
  ggplot(aes(date, cases, color = series)) +
  geom_line(linewidth = 1) +
  geom_point(size = 2) +
  scale_color_manual(values = c(
    observed = "#21313c",
    explanation = "#2a9d8f",
    prediction = "#e76f51"
  )) +
  labs(
    title = "Held-Out Forecast Evaluation",
    subtitle = "Last six simulated Bandung months are treated as future observations",
    x = NULL,
    y = "Monthly cases",
    color = NULL
  )

This figure compares observed future cases with two model forecasts. If a model follows the held-out trajectory more closely, it has better short-run predictive performance. In a real surveillance system, this evaluation must be repeated over many rolling forecast origins.

Forecasting procedure.

A practical disease-risk forecasting workflow should answer the following questions:

1. What is the forecast target: cases, SIR, relative risk, or exceedance probability?
2. What is the forecast horizon: one month ahead, three months ahead, or one year ahead?
3. What data are available before the forecast date?
4. Which covariates are known in advance and which must themselves be forecast?
5. How will forecast accuracy be evaluated?
6. How will forecast uncertainty be communicated to decision-makers?
7. What action threshold will trigger public-health response?

Important. A model should not use future information when producing a forecast. For example, if rainfall for next month is not known at decision time, it must be forecast or excluded from an operational early-warning model.

A practical forecasting procedure for spatiotemporal panel data is:

1. Choose the forecast horizon, for example 1, 3, or 6 months ahead.
2. Hold out recent periods to evaluate predictive performance.
3. Fit the candidate models using only the training period.
4. Generate posterior predictive distributions for future periods.
5. Compare predictive accuracy using MAE, RMSE, calibration, and spatial hotspot classification.
6. Map future relative risk and posterior exceedance probabilities.

For Bayesian INLA models, future cases are set to NA, and the model returns posterior predictions for those missing outcomes. The quick example below uses the fitted prediction GLM to show the mechanics.

date	mean_rr	max_rr	districts_above_125	probable_hotspots
2027-01-01	1.501	3.270	19	12
2027-02-01	1.940	4.227	23	20
2027-03-01	1.876	4.087	22	20
2027-04-01	1.458	3.177	17	12
2027-05-01	0.962	2.095	7	4
2027-06-01	0.572	1.246	0	0

future_course_map <- districts |>
  select(district_id, district) |>
  left_join(
    future_course |>
      filter(h %in% c(1, 3, 6)) |>
      mutate(month_label = format(date, "%Y-%m")),
    by = "district_id"
  )

ggplot(future_course_map) +
  geom_sf(aes(fill = rr_forecast), color = "white", linewidth = 0.25) +
  facet_wrap(~month_label, nrow = 1) +
  scale_fill_gradient2(
    low = "#f1d9b1", mid = "#fffaf3", high = "#973d1c",
    midpoint = 1, name = "Forecast RR"
  ) +
  labs(
    title = "Forecast Relative Risk Map",
    subtitle = "Forecasting is evaluated by future performance, not only in-sample fit"
  )

Day 4 exercise. Repeat the train/test split with a rolling origin: train through month 24 and forecast 25-30, then train through month 30 and forecast 31-36. Compare MAE, RMSE, and hotspot classification stability.

Day 4.4 The third paper: competing forecast specifications

Paper 3. Jaya and Folmer (2025), Spatiotemporal forecasting models with and without a confounded covariate, Journal of Geographical Systems 27, 113–146, DOI: 10.1007/s10109-024-00454-z. In this paper the methodological comparison is between (M1) a multivariate model with the potentially confounded covariate, (M2) a univariate pure latent spatiotemporal model, and (M3) a full model incorporating both. The paper’s Monte Carlo study encompasses continuous, count and binary responses. The simplified Poisson laboratory below illustrates this comparison using the same simulated Bandung kecamatan panel developed in Day 1.

Reading guide for the paper and PPTX Part 2. Slides 43–48 introduce pure versus causal forecasting and the research problem; slides 49–54 explain the GMRF model and Bayesian forecast; slides 55–57 describe the Monte Carlo experiment; slides 58–66 discuss the forecasting results and conclusions. This is a focused methodological case study, not a claim that one model wins for every dataset.

The forecasting question. Suppose a public-health analyst wants to predict dengue-related outcomes in each Bandung kecamatan several months ahead. Rainfall can help explain an observed outbreak, but future rainfall may itself need to be forecast. A spatially patterned covariate may also partially reproduce an unobserved, spatially patterned disease process. Consequently, a regression coefficient that appears useful for explaining historical observations does not guarantee a better out-of-sample forecast. The paper asks whether incorporating such a confounded covariate improves prediction over a pure spatiotemporal model that learns directly from the outcome’s spatial and temporal dependence.

Three distinct sources of difficulty. First, if \(x_{it}\) is highly correlated with an unobserved risk component \(\delta_{it}\), the regression and latent effect compete to explain similar variation. Second, the future value \(x_{i,T+\ell}\) may be unknown; replacing it with a forecast \(\widehat{x}_{i,T+\ell}\) introduces additional prediction error. Third, extrapolating spatial and temporal random effects requires a coherent dependence model: including an unrelated random intercept for every future period does not automatically generate a good forecast. The paper investigates these problems jointly in a Bayesian latent-Gaussian framework.

The central model and its meaning. For area \(i=1,\ldots,n\) and month \(t=1,\ldots,T\), the original study considers the linear predictor

\[ \eta_{it}=\beta_0+\beta_1x_{it}+\omega_i+\phi_t+\delta_{it}, \]

where \(\omega_i\) captures persistent spatial variation, \(\phi_t\) captures common temporal change, and \(\delta_{it}\) captures the space–time interaction. The observation model depends on the outcome:

\[ \begin{aligned} \text{Continuous: }&Y_{it}\mid\eta_{it}\sim N(\eta_{it},\sigma_y^2),\\ \text{Count: }&Y_{it}\mid\eta_{it}\sim\operatorname{Poisson}\{\exp(\eta_{it})\},\\ \text{Binary: }&Y_{it}\mid\eta_{it}\sim\operatorname{Bernoulli}\{\operatorname{logit}^{-1}(\eta_{it})\}. \end{aligned} \]

The same underlying idea therefore applies to a continuous risk index, number of cases, or a yes/no event, but the link function, forecast target and uncertainty are different. The workshop’s earlier dengue example uses an additional \(\log E_{it}\) offset to forecast relative risk; that offset is an adaptation for disease mapping, not the exact count-generation equation used in Paper 3’s Monte Carlo experiment.

What exactly are M1, M2 and M3? The article uses the following nested formulations. In particular, M1 omits all three structured effects; M2 retains latent spatiotemporal structure but sets the covariate coefficient to zero; M3 combines the two.

\[ \begin{aligned} M_1:&\quad\eta_{it}=\beta_0+\beta_1 x_{it},\\ M_2:&\quad\eta_{it}=\beta_0+\omega_i+\phi_t+\delta_{it},\\ M_3:&\quad\eta_{it}=\beta_0+\beta_1 x_{it}+\omega_i+\phi_t+\delta_{it}. \end{aligned} \]

This distinction is conceptually important: M2 is not simply an intercept-only model, and M3 is not guaranteed to win merely because it has more terms. We later add a common seasonality/trend to the small Bandung classroom implementation to make future-month extrapolation easy to understand; those teaching formulas should not be mistaken for the article’s exact M1–M3 definitions.

What does confounding mean here? In the paper, confounding arises when a covariate shares variation with one or more spatial, temporal or spatiotemporal latent effects. The focused Monte Carlo experiment correlates the covariate with the standardized interaction component:

\[ x_{it}=\rho_\delta\widetilde{\delta}_{it} +\sqrt{1-\rho_\delta^2}\,z_{it}, \qquad z_{it}\sim N(0,1), \]

where \(\widetilde{\delta}\) is standardized and \(\rho_\delta\in\{0.3,0.6,0.9\}\) describes low, moderate and high confounding. When \(\rho_\delta\) increases, \(x\) resembles the hidden interaction more closely; it becomes more difficult to distinguish their predictive contributions. This is statistical confounding of spatial patterns, not by itself proof of a causal relation.

A calculation participants can do without R. Let a standardized hidden risk at one district-month equal \(\widetilde{\delta}=1\), the independent shock equal \(z=0.5\), and take \(\rho_\delta=0.6\). Then

\[ x=0.6(1)+\sqrt{1-0.6^2}(0.5) =0.6+0.8(0.5)=\boxed{1.0}. \]

At \(\rho_\delta=0.9\), the same pair gives \(x=0.9+\sqrt{0.19}(0.5)\approx1.118\). In the second scenario, a larger share of the covariate’s systematic component is tied to the hidden spatial–temporal pattern. The numerical comparison uses the same latent effect and independent shock so the effect of changing \(\rho_\delta\) is visible.

set.seed(20261124)
hidden_effect <- rnorm(180)
independent_shock <- rnorm(180)
rho_values <- c(0.3, 0.6, 0.9)
par(mfrow = c(1, 3), mar = c(4, 4, 2.5, 1))
for (rho in rho_values) {
  x <- rho * hidden_effect + sqrt(1 - rho^2) * independent_shock
  plot(hidden_effect, x, pch = 19, cex = .45,
       col = grDevices::adjustcolor("#B16B36", alpha.f = .45),
       xlab = "Latent space-time effect", ylab = "Covariate x",
       main = paste0("Confounding = ", rho))
  abline(lm(x ~ hidden_effect), col = "#693C28", lwd = 2)
}

par(mfrow = c(1, 1))

The three panels are a simulation illustration. They are not figures reproduced from the paper. Participants should observe that a high empirical correlation can arise by construction even when the independent shock remains present.

Bayesian forecasting rather than fitted-value substitution. Let \(\boldsymbol\theta\) collect the coefficients, latent effects and hyperparameters, and let \(\mathcal D_T\) contain observations through month \(T\). For a future outcome,

\[ p(Y_{i,T+\ell}\mid\mathcal D_T) =\int p(Y_{i,T+\ell}\mid\boldsymbol\theta,\mathcal D_T) p(\boldsymbol\theta\mid\mathcal D_T)\,d\boldsymbol\theta. \]

This integral averages future outcomes over posterior uncertainty. For a Poisson outcome, for example, \(E(Y^*\mid\mathcal D_T)=E\{\exp(\eta^*)\mid\mathcal D_T\}\), which generally differs from \(\exp\{E(\eta^*\mid\mathcal D_T)\}\). In R-INLA we append rows with an NA response for future months, define their spatial and temporal indices, and request fitted and latent marginals. Future covariates require their own forecast or explicitly stated scenario. A point forecast alone does not summarize posterior predictive uncertainty.

Why use a GMRF? A Gaussian Markov random field has a sparse precision matrix \(Q=\Sigma^{-1}\) reflecting conditional dependence. Adjacent districts and nearby months can be coupled while most distant pairs are conditionally independent given neighbors. The paper develops spatial and temporal precision structures and a space–time interaction with a Kronecker product, schematically \(Q_\delta\propto Q_\phi\otimes Q_\omega\). INLA exploits this sparsity to approximate Bayesian posterior marginals without requiring a full MCMC simulation for every parameter. Model identification, the constraints on structured effects and the time-process specification remain substantive modelling choices.

What was varied in the published Monte Carlo? The paper does not base its conclusion on one held-out Bandung map. Its simulated study varies (i) spatial units \(n=16,49,100\) for continuous and count outcomes, or \(n=16,25,49\) for binary outcomes; (ii) time points \(T=10,20,30\); (iii) forecast horizons \(\ell=2\) and \(\ell=8\); (iv) spatiotemporal interaction parameter \(\lambda=0.3,0.6,0.9\); and (v) confounding strength \(\rho_\delta=0.3,0.6,0.9\). The coefficient is \(\beta_1=0.5\). The intercept is \(10\) for Gaussian, \(1\) for Poisson, and \(0.1\) for Bernoulli data. For the Gaussian response, its simulated conditional observation variance is \(0.01\). Paper 3 uses weak Gaussian coefficient priors and hyperpriors for the latent-precision parameters; the precise choice and sensitivity should be reported when applying the approach to real surveillance data.

A worked forecast-loss example. Suppose the four observed counts in a short holdout are \((2,4,3,5)\); M2 predicts \((2.5,3.5,3.5,4.5)\); and M3 predicts \((3,5,3,6)\). Then

\[ \begin{aligned} \operatorname{MSPE}(M_2)&=\frac{(.5)^2+(-.5)^2+(.5)^2+(-.5)^2}{4}=\boxed{0.25},\\ \operatorname{MSPE}(M_3)&=\frac{1^2+1^2+0^2+1^2}{4}=\boxed{0.75}. \end{aligned} \]

This toy calculation shows how to compare forecasts; the invented four values are not results from the article. We cannot decide which model is best from training residuals or in-sample DIC alone; we must score future or genuinely held-out observations.

observed <- c(2, 4, 3, 5)
forecast_M2 <- c(2.5, 3.5, 3.5, 4.5)
forecast_M3 <- c(3, 5, 3, 6)
c(M2_MSPE = mean((observed - forecast_M2)^2),
  M3_MSPE = mean((observed - forecast_M3)^2))

## M2_MSPE M3_MSPE 
##    0.25    0.75

What does the paper actually find? Read the following comparisons as results conditional on its simulated data-generating mechanism, not as a universal rule about regression covariates.

Response	Reported comparison and interpretation
Continuous (Gaussian)	M3 can forecast better under weak confounding because the covariate contributes independent information. Its advantage shrinks as the interaction and confounding become stronger; with both \(\lambda=0.9\) and \(\rho_\delta=0.9\), M2 outperforms M3 in the reported combinations. Longer horizons generally raise MSPE.
Count (Poisson)	M1 has the highest MSPE in the examined settings. M2 consistently improves over M3 in the article’s count simulations; errors in forecasting the confounded covariate can degrade the full model. M2 is comparatively stable across the investigated horizons.
Binary (Bernoulli)	M1 again has the highest MSPE. M2 is usually at least as accurate as M3, but the difference is often small or zero; one should not claim a large universal binary-data advantage. The paper first predicts event probability and then classifies at a 0.5 threshold for its binary MSPE.

Interpret the Gaussian exception. When \(x_{it}\) has a component that cannot be reconstructed from latent spatial and temporal dependence, retaining it in M3 can help. If the same covariate mostly reproduces the latent process, the additional fitted coefficient and the need to forecast \(x\) can increase uncertainty without contributing much new forecast information. Thus, the research question is incremental predictive value after spatial–temporal structure, not merely whether \(\beta_1\) is statistically different from zero.

Interpret count and binary losses. Poisson outcomes are sensitive to the exponential link: a moderate error on the log scale can become a disproportionately large error in the expected count. For binary data, probability forecasts should ideally be assessed with Brier score, log score and calibration plots, in addition to a classification loss. Paper 3’s reported binary MSPE is calculated after thresholding predicted probabilities at 0.5; our short exercise below provides both a probability-based Brier score and an RMSE-style count/continuous score so that the two tasks are not confused.

An accessible three-family simulation (ordinary R, no long INLA fit). We now create 30 simulated Bandung-area effects observed over 36 months, correlate the explanatory variable with that hidden effect, and reserve the last six months. For readability, the simple classroom versions of M2 and M3 use district indicators and a time trend as a proxy for the paper’s full GMRF. The three examples below are executable and demonstrate training/testing; they do not attempt to reproduce the article’s Monte Carlo results or its exact rankings.

set.seed(20261125)
n_area <- 30
n_month <- 36
demo <- expand.grid(area = factor(seq_len(n_area)), month = seq_len(n_month))
demo <- demo[order(demo$area, demo$month), ]
area_effect <- rnorm(n_area, sd = .35)
demo$hidden <- area_effect[as.integer(demo$area)] +
  .20 * sin(2 * pi * demo$month / 12) + rnorm(nrow(demo), sd = .15)
hidden_z <- as.numeric(scale(demo$hidden))
rho_demo <- .6
demo$x <- rho_demo * hidden_z +
  sqrt(1 - rho_demo^2) * rnorm(nrow(demo))
demo$trend <- demo$month / 12
demo$season <- sin(2 * pi * demo$month / 12)
demo$continuous <- 10 + .5 * demo$x + .5 * demo$hidden +
  rnorm(nrow(demo), sd = .1)
demo$count <- rpois(nrow(demo), exp(1 + .5 * demo$x + .4 * demo$hidden))
demo$binary <- rbinom(nrow(demo), size = 1,
  prob = plogis(.1 + .5 * demo$x + .4 * demo$hidden))
training <- demo$month <= 30

compare_family <- function(response, type) {
  family_used <- switch(type, gaussian = gaussian(),
                        count = poisson(), binary = binomial())
  formulas <- list(
    M1 = reformulate(c("x", "trend", "season"), response),
    M2 = reformulate(c("area", "trend", "season"), response),
    M3 = reformulate(c("x", "area", "trend", "season"), response))
  do.call(rbind, lapply(names(formulas), function(name) {
    model <- glm(formulas[[name]], data = demo[training, ],
                 family = family_used)
    predicted <- predict(model, newdata = demo[!training, ],
                         type = "response")
    actual <- demo[[response]][!training]
    data.frame(outcome = type, model = name,
               score = mean((actual - predicted)^2))
  }))
}
simulation_scores <- rbind(compare_family("continuous", "gaussian"),
                           compare_family("count", "count"),
                           compare_family("binary", "binary"))
simulation_scores$score <- round(simulation_scores$score, 4)
print(simulation_scores, row.names = FALSE)

##   outcome model  score
##  gaussian    M1 0.0317
##  gaussian    M2 0.2249
##  gaussian    M3 0.0131
##     count    M1 2.8012
##     count    M2 5.7160
##     count    M3 2.8201
##    binary    M1 0.2435
##    binary    M2 0.2724
##    binary    M3 0.2552

All three score values are mean squared prediction errors on different outcome scales. The continuous and count outcomes use squared errors of posterior/conditional-mean proxies; the binary outcome uses squared errors of probabilities, which is the Brier score. Compare M1, M2 and M3 within an outcome, not scores across different outcomes. The original paper thresholds binary probabilities and uses its own Monte Carlo design, so numerical scores are not directly comparable to its reported binary MSPE.

Map a model’s forecast for Bandung. The exercise’s first panel below shows the average simulated held-out count in each of the 30 kecamatan, and the second shows the M2 teaching-model forecast. Their overall geographic pattern is easier to compare than 180 printed rows. Polygon boundaries come from the supplied Bandung GIS data; outcome values are simulated, not real dengue counts.

fit_demo_M2 <- glm(count ~ area + trend + season,
                   data = demo[training, ], family = poisson())
heldout_demo <- demo[!training, ]
heldout_demo$forecast <- predict(fit_demo_M2, heldout_demo, type = "response")
area_observed <- aggregate(count ~ area, heldout_demo, mean)
area_predicted <- aggregate(forecast ~ area, heldout_demo, mean)
map_values <- data.frame(district_id = seq_len(n_area),
                         observed = area_observed$count,
                         predicted = area_predicted$forecast)
bandung_forecasts <- dplyr::left_join(districts, map_values, by = "district_id")
stopifnot(nrow(bandung_forecasts) == n_area,
          all(is.finite(bandung_forecasts$predicted)))
map_long <- tidyr::pivot_longer(bandung_forecasts,
  cols = c("observed", "predicted"),
  names_to = "Series", values_to = "Average_cases")
ggplot2::ggplot(map_long) +
  ggplot2::geom_sf(ggplot2::aes(fill = Average_cases),
                   colour = "#FFFFFF", linewidth = .25) +
  ggplot2::facet_wrap(~Series) +
  ggplot2::scale_fill_gradient(low = "#FFF2DE", high = "#A9502B") +
  ggplot2::labs(title = "Simulated Bandung: six-month count forecast",
                subtitle = "Observed holdout mean vs M2 teaching approximation",
                fill = "Mean count") +
  ggplot2::theme_minimal(base_size = 11)

A reproducible sensitivity exercise. Change rho_demo <- .6 to .3 and .9 in the simulation, then compare the three MSEs. Repeat with a different random seed and a longer forecast horizon. A change in ranking is a useful lesson: the observed winner depends on the strength and structure of spatial dependence, the covariate’s independent information, the observation family, forecasting the covariate correctly, and the available history. It does not invalidate a systematically repeated published simulation.

Discussion with participants. (1) Why can M3 explain the training data well yet forecast poorly? (2) Would M3’s advantage change if all future values of \(x\) were known perfectly? (3) Why is a 0.5 classification threshold a different decision from forecasting event probability? (4) Which GMRF components allow an unobserved future time index to be extrapolated? (5) What would you report to a health department: an expected-case map only, or predictive intervals, exceedance probabilities and calibration as well?

Forecast information and experimental design.

Let the training period be months 1–\((T-6)\) and the test period be the last six months. The participants should predict cases without conditioning on future observed cases. Here the future weather proxies are supplied as known simulation scenarios. In a practical forecasting operation, their own forecasts and uncertainty must be provided separately. The latent persistent district effect can be learned from training data, while the linear time trend is extrapolated to future months. This teaching implementation deliberately differs from the paper’s full simulation design and should not be presented as an exact replication.

\[\begin{aligned} M_1: &\quad \log(\mu_{it}/E_{it})=\alpha+\beta x_{it}+s(t),\\ M_2: &\quad \log(\mu_{it}/E_{it})=\alpha+u_i+s(t),\\ M_3: &\quad \log(\mu_{it}/E_{it})=\alpha+\beta x_{it}+u_i+s(t). \end{aligned}\]

For the simple experiment \(s(t)\) consists of seasonal terms and a linear extrapolating trend. The paper’s full model can employ richer latent GMRF constructions. This distinction matters because an INLA model that fits history well need not extrapolate a time-indexed random interaction correctly.

##   model      MAE     RMSE
## 1    M1 10.28237 13.19566
## 2    M2 10.27133 13.11614
## 3    M3 10.27182 13.11661

Interpret the comparison.

The held-out MAE and RMSE evaluate the mean count forecasts. They do not assess the calibration of full posterior predictive intervals. Lower is better for both scores, but results are conditional on the simulation, the information set and the fitted model class. The original article reports general findings from a Monte Carlo design that should be discussed separately from these educational numbers.

Exercise. Refit using a different six-month split and compare the winner. Discuss how prior strength and district heterogeneity affect M2/M3. Add a predictive-interval coverage study using marginal or posterior predictive approximations, and distinguish covariate confounding from information leakage.

Paper 3 figures: original results and executable illustrations

The four high-resolution figures below are reproduced directly from the figures supplied with the workshop slides. They are published/presentation results, not outputs of the short Bandung R exercise. Keeping them next to the discussion is essential: the participant must be able to inspect the factorial Monte Carlo outcomes and the forecast trajectories without opening PowerPoint. The full slide sequence follows immediately after the R example.

Original continuous-outcome MSE results. Columns vary interaction strength \(\lambda\) and forecasting horizon \(\ell\); rows vary numbers of regions \(n\) and time periods \(T\). Within a panel, confounding \(\rho_\delta\) varies along the horizontal axis. The comparison is between M2 (light) and M3 (orange). The pattern is conditional on the interaction, confounding and horizon; it is not enough to judge one isolated bar.

Full-size original continuous-data simulation MSE, M2 versus M3

Original continuous-outcome forecast trajectories. Black represents observed outcomes, while red, green and blue show predictions from M1, M2 and M3. The grey interval indicates the forecast horizon. Rows correspond to selected regions and confounding levels; columns correspond to interaction strengths. Examine forecast-period disagreement, rather than which colored line is closest during training.

Full-size original continuous-data observed and forecast trajectories

Original count-outcome MSE results. The same factorial comparison shows that M2 has lower MSPE than M3 in the paper’s investigated Poisson settings. The additional uncertainty in predicting a confounded covariate can affect the Poisson mean through its exponential link.

Full-size original count-data simulation MSE, M2 versus M3

Original count-outcome forecast trajectories. Spatially structured effects can recover recurring area-specific count variation. Under the paper’s count simulation, the extra confounded covariate is often unnecessary for forecasting, although it may still have a substantive explanatory role.

Full-size original count-data observed and forecast trajectories

Now generate a comparable teaching figure in R. Use the simulated Bandung data already created above. The next chunk refits three familiar GLMs for the continuous and count outcomes, makes genuine out-of-sample forecasts for months 31–36, and draws the observed series together with M1, M2 and M3. Its line colors intentionally match the workshop’s original figure (black, red, green, blue). District numbers 1 and 5 are two of the 30 actual Bandung polygons, but all response measurements are simulated. Unlike the paper, the M2/M3 implementations use district fixed effects as an easily understood approximation to latent spatial effects.

training <- demo$month <= 30
predict_three <- function(outcome, family_used) {
  model_M1 <- glm(reformulate(c("x", "trend", "season"), outcome),
                  data = demo[training, ], family = family_used)
  model_M2 <- glm(reformulate(c("area", "trend", "season"), outcome),
                  data = demo[training, ], family = family_used)
  model_M3 <- glm(reformulate(c("x", "area", "trend", "season"), outcome),
                  data = demo[training, ], family = family_used)
  subset_data <- demo[demo$area %in% c("1", "5"), ]
  subset_data$Observed <- subset_data[[outcome]]
  subset_data$M1 <- predict(model_M1, subset_data, type = "response")
  subset_data$M2 <- predict(model_M2, subset_data, type = "response")
  subset_data$M3 <- predict(model_M3, subset_data, type = "response")
  long <- tidyr::pivot_longer(subset_data,
    cols = c("Observed", "M1", "M2", "M3"),
    names_to = "Model", values_to = "Value")
  long$Outcome <- outcome
  long$District <- paste("Bandung district", as.character(long$area))
  long
}
demo_lines <- rbind(
  predict_three("continuous", gaussian()),
  predict_three("count", poisson()))
ggplot2::ggplot(demo_lines, ggplot2::aes(month, Value, colour = Model)) +
  ggplot2::annotate("rect", xmin = 30.5, xmax = 36.5,
                    ymin = -Inf, ymax = Inf, fill = "grey65", alpha = .2) +
  ggplot2::geom_line(linewidth = .65) +
  ggplot2::facet_grid(Outcome ~ District, scales = "free_y") +
  ggplot2::scale_colour_manual(values = c(Observed = "black",
    M1 = "red3", M2 = "forestgreen", M3 = "blue3")) +
  ggplot2::labs(title = "Bandung classroom simulation: observed vs M1/M2/M3",
    subtitle = "Shaded months 31–36 are withheld; not a numerical replication of Paper 3",
    x = "Month", y = "Response / forecast") +
  ggplot2::theme_minimal(base_size = 11)

For a numerical comparison of those same models, the following computation produces actual holdout MSPE values, rather than displaying code without a calculated outcome. The results depend on the set.seed() value in the preceding simulation. MSE on the Gaussian and Poisson scales must be interpreted separately.

print(simulation_scores[simulation_scores$outcome %in% c("gaussian", "count"), ],
      row.names = FALSE)

##   outcome model  score
##  gaussian    M1 0.0317
##  gaussian    M2 0.2249
##  gaussian    M3 0.0131
##     count    M1 2.8012
##     count    M2 5.7160
##     count    M3 2.8201

Important computational distinction. The four original slide figures summarize a much larger Monte Carlo study, with multiple \(n\), \(T\), \(\lambda\), \(\rho_\delta\) and \(\ell\) combinations as described in the paper. The runnable code here deliberately uses a smaller, transparent dataset suitable for participants unfamiliar with R. Its purpose is to demonstrate the mechanics of simulation, fitting, forecasting and MSPE, not to regenerate the exact paper figures or claim that a short GLM calculation is the original Bayesian GMRF experiment.

References

1. Jaya, I. G. N. M. & Folmer, H. (2020). Bayesian spatiotemporal mapping of relative dengue disease risk in Bandung, Indonesia. Journal of Geographical Systems, 22, 105–142. DOI: 10.1007/s10109-019-00311-4.
2. Jaya, I. G. N. M. & Folmer, H. (2021). Identifying spatiotemporal clusters by means of agglomerative hierarchical clustering and Bayesian regression analysis with spatiotemporally varying coefficients. Geographical Analysis, 53, 767–817. DOI: 10.1111/gean.12264.
3. Jaya, I. G. N. M. & Folmer, H. (2025). Spatiotemporal forecasting models with and without a confounded covariate. Journal of Geographical Systems, 27, 113–146. DOI: 10.1007/s10109-024-00454-z.
4. Rue, H., Martino, S. & Chopin, N. (2009). Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. JRSS B, 71, 319–392.
5. Knorr-Held, L. (2000). Bayesian modelling of inseparable space-time variation in disease risk. Statistics in Medicine, 19, 2555–2567.
6. R-INLA Project. https://www.r-inla.org/download/.
7. Gómez-Rubio, V. (2020). Bayesian Inference with INLA. Chapter 1, “Introduction to Bayesian Inference.” https://becarioprecario.bitbucket.io/inla-gitbook/ch-intro.html. The workshop draws on its concise explanation of Bayes’ theorem, normal conjugacy, posterior marginals, and the introductory role of INLA; the Bandung calculations and explanations are adapted specifically for this course.