A synthetic WNV / Europe demonstration with an identifiability check
Purpose. A self-contained, reproducible R workflow that (i) simulates a realistic five-state West Nile virus emergence process across NUTS3-scale regions placed at real European coordinates, (ii) fits the multi-state survival model, and (iii) verifies that the estimated transition hazards recover the known simulated effects. This is the “identifiability pre-work” that answers the reviewer’s central technical question before any real-data analysis.
WNV biology. Humans are dead-end hosts, so human mobility is not an introduction driver for WNV. The introduction barrier (S2→S3) is driven by avian-flyway intensity, trade connectivity, and spatial diffusion (short-range bird/mosquito spread encoded by distance-to-focus and neighbourhood pressure). Human mobility fills that same structural slot only for Aedes-borne viruses.
pkgs <- c("ggplot2", "dplyr", "tidyr", "survival", "scales", "knitr")
need <- pkgs[!vapply(pkgs, requireNamespace, logical(1), quietly = TRUE)]
if (length(need)) install.packages(need, repos = "https://cloud.r-project.org")
invisible(lapply(pkgs, library, character.only = TRUE))
set.seed(20260721)~336 pseudo-NUTS3 regions are drawn as jittered clusters around real European country centroids (so the spread map reads as Europe). Coordinates are converted to kilometres by an equirectangular projection for distance calculations.
clusters <- data.frame(
lon = c(-8.2,-3.5, 2.3,-1.6,-8.0, 5.0,10.2, 9.6, 9.5,15.5,25.5,19.3,17.5,11.5,12.6,19.5,25.0,25.3,19.5,22.5,24.5),
lat = c(39.7,40.2,47.0,52.8,53.3,51.2,51.0,56.0,61.0,60.5,63.0,52.2,49.3,47.3,42.8,47.1,46.0,42.8,44.3,39.3,57.0),
n = c( 10, 24, 30, 22, 8, 12, 30, 8, 12, 16, 12, 20, 12, 12, 26, 10, 16, 10, 18, 12, 12))
lon <- unlist(Map(function(m,n) rnorm(n, m, 1.1), clusters$lon, clusters$n))
lat <- unlist(Map(function(m,n) rnorm(n, m, 0.8), clusters$lat, clusters$n))
N <- length(lon)
R_EARTH <- 6371; lat0 <- mean(lat) * pi/180
xk <- lon * pi/180 * R_EARTH * cos(lat0) # km east
yk <- lat * pi/180 * R_EARTH # km north
zc <- function(v) (v - mean(v)) / sd(v)
humidity_z <- rnorm(N); urban_z <- rnorm(N); density_z <- rnorm(N); socio_z <- rnorm(N)
avian_flyway <- zc(-0.05*(lat-46) + 0.03*(lon-12) + rnorm(N)) # higher in SE / Danube corridor
trade_z <- zc(-0.04*abs(lat-50) - 0.03*abs(lon-8) + rnorm(N)) # higher in central/western core
temp_base <- 27 - 0.48*(lat-36) + rnorm(N, 0, 0.8) # deg C: ~27 (south) .. 14 (north)
## TRUE coefficients (log-hazard) per transition -----------------------------
B1 <- c(intercept=-1.9, temp_c=0.22, humidity_z=0.15, urban_z=0.25) # ecological
B2 <- c(intercept=-2.4, avian_flyway=0.55, trade_z=0.30, nbr_pres=0.45, dist100=-0.58) # introduction
B3 <- c(intercept=-1.3, temp_c=0.16, nbr_pres=0.25, density_z=0.15) # transmission
B4 <- c(intercept=-1.6, temp_c=0.10, socio_z=0.30) # persistence
LAMBDA <- 250 # km distance-decay kernel
cloglog_p <- function(eta) 1 - exp(-exp(eta))
D <- as.matrix(dist(cbind(xk, yk))); diag(D) <- InfStates: 1 non-receptive → 2 ecological receptivity → 3 introduction (avian / trade / diffusion propagule pressure) → 4 local transmission → 5 persistence. Five seed regions start in state 4 in the historical WNV origin (Romania/Hungary).
state <- rep(1L, N)
seed <- order((lon-24)^2 + (lat-46)^2)[1:5]; state[seed] <- 4L # SE origin
panel <- list(); spread <- list()
for (yr in 2006:2024) {
temp_c <- (temp_base + 0.06*(yr-2006) + rnorm(N, 0, 0.4)) - 16
invaded <- state %in% c(4, 5)
if (any(invaded)) {
Dsub <- D[, invaded, drop = FALSE]
dist100 <- apply(Dsub, 1, min) / 100
nbr_pres <- rowSums(exp(-Dsub / LAMBDA))
} else { dist100 <- rep(3000/100, N); nbr_pres <- rep(0, N) }
ev <- integer(N); new_state <- state
for (i in seq_len(N)) {
s <- state[i]
eta <- switch(s,
B1["intercept"] + B1["temp_c"]*temp_c[i] + B1["humidity_z"]*humidity_z[i] + B1["urban_z"]*urban_z[i],
B2["intercept"] + B2["avian_flyway"]*avian_flyway[i] + B2["trade_z"]*trade_z[i] +
B2["nbr_pres"]*nbr_pres[i] + B2["dist100"]*dist100[i],
B3["intercept"] + B3["temp_c"]*temp_c[i] + B3["nbr_pres"]*nbr_pres[i] + B3["density_z"]*density_z[i],
B4["intercept"] + B4["temp_c"]*temp_c[i] + B4["socio_z"]*socio_z[i],
-Inf) # state 5 absorbing
if (is.finite(eta) && runif(1) < cloglog_p(eta)) { ev[i] <- 1L; new_state[i] <- s + 1L }
}
panel[[as.character(yr)]] <- data.frame(
region = seq_len(N), year = yr, state_start = state,
temp_c, humidity_z, urban_z, avian_flyway, trade_z, density_z, socio_z,
dist100, nbr_pres, event = ev)
state <- new_state
spread[[as.character(yr)]] <- state
}
panel <- do.call(rbind, panel)trans <- list(
"B1: S1 to S2 (Ecological)" = list(origin = 1, covs = c("temp_c","humidity_z","urban_z"), truth = B1),
"B2: S2 to S3 (Introduction)" = list(origin = 2, covs = c("avian_flyway","trade_z","nbr_pres","dist100"), truth = B2),
"B3: S3 to S4 (Transmission)" = list(origin = 3, covs = c("temp_c","nbr_pres","density_z"), truth = B3),
"B4: S4 to S5 (Persistence)" = list(origin = 4, covs = c("temp_c","socio_z"), truth = B4))
pretty <- c(temp_c="Temperature (C)", humidity_z="Humidity (SD)", urban_z="Urban land (SD)",
avian_flyway="Avian flyway (SD)", trade_z="Trade connectivity (SD)",
nbr_pres="Neighbourhood pressure", dist100="Distance to focus (/100 km)",
density_z="Population density (SD)", socio_z="Socioeconomic (SD)")The data-generating process is a discrete-time proportional-hazards model, so the exact estimator is a complementary-log-log GLM fitted to each transition’s at-risk rows. Together these transition-specific hazards constitute the Markov multi-state model.
est <- do.call(rbind, lapply(names(trans), function(tn) {
tt <- trans[[tn]]
sub <- panel[panel$state_start == tt$origin, ]
f <- as.formula(paste("event ~", paste(tt$covs, collapse = " + ")))
m <- glm(f, family = binomial("cloglog"), data = sub)
co <- summary(m)$coefficients
b <- co[tt$covs, "Estimate"]; se <- co[tt$covs, "Std. Error"]
data.frame(transition = tn, covariate = tt$covs, label = pretty[tt$covs],
true_beta = tt$truth[tt$covs], est_beta = b, se = se,
true_HR = exp(tt$truth[tt$covs]), est_HR = exp(b),
lo = exp(b - 1.96*se), hi = exp(b + 1.96*se), row.names = NULL)
}))
est$covered <- est$true_HR >= est$lo & est$true_HR <= est$hi
cat(sprintf("Mean |beta error| = %.3f | True HR within 95%% CI: %d/%d\n",
mean(abs(est$est_beta - est$true_beta)), sum(est$covered), nrow(est)))Mean |beta error| = 0.053 | True HR within 95% CI: 12/12knitr::kable(est[, c("transition","label","true_HR","est_HR","lo","hi","covered")],
digits = 2, caption = "Estimated vs. true hazard ratios by transition.")| transition | label | true_HR | est_HR | lo | hi | covered |
|---|---|---|---|---|---|---|
| B1: S1 to S2 (Ecological) | Temperature (C) | 1.25 | 1.29 | 1.24 | 1.34 | TRUE |
| B1: S1 to S2 (Ecological) | Humidity (SD) | 1.16 | 1.18 | 1.06 | 1.32 | TRUE |
| B1: S1 to S2 (Ecological) | Urban land (SD) | 1.28 | 1.28 | 1.14 | 1.43 | TRUE |
| B2: S2 to S3 (Introduction) | Avian flyway (SD) | 1.73 | 1.69 | 1.42 | 2.00 | TRUE |
| B2: S2 to S3 (Introduction) | Trade connectivity (SD) | 1.35 | 1.26 | 1.07 | 1.49 | TRUE |
| B2: S2 to S3 (Introduction) | Neighbourhood pressure | 1.57 | 1.56 | 1.42 | 1.71 | TRUE |
| B2: S2 to S3 (Introduction) | Distance to focus (/100 km) | 0.56 | 0.53 | 0.45 | 0.62 | TRUE |
| B3: S3 to S4 (Transmission) | Temperature (C) | 1.17 | 1.12 | 0.97 | 1.31 | TRUE |
| B3: S3 to S4 (Transmission) | Neighbourhood pressure | 1.28 | 1.41 | 1.17 | 1.71 | TRUE |
| B3: S3 to S4 (Transmission) | Population density (SD) | 1.16 | 1.44 | 1.01 | 2.05 | TRUE |
| B4: S4 to S5 (Persistence) | Temperature (C) | 1.11 | 1.12 | 1.03 | 1.21 | TRUE |
| B4: S4 to S5 (Persistence) | Socioeconomic (SD) | 1.35 | 1.43 | 1.22 | 1.67 | TRUE |
The same transitions as a counting-process Cox model (survival::coxph), the canonical multi-state-survival formulation, for cross-checking:
cp <- panel; cp$tstart <- cp$year - 2006; cp$tstop <- cp$tstart + 1
cox_tab <- do.call(rbind, lapply(names(trans), function(tn) {
tt <- trans[[tn]]; sub <- cp[cp$state_start == tt$origin, ]
f <- as.formula(paste("Surv(tstart, tstop, event) ~", paste(tt$covs, collapse = " + ")))
m <- coxph(f, data = sub, ties = "efron")
data.frame(transition = tn, covariate = tt$covs,
coxph_HR = exp(coef(m)[tt$covs]), true_HR = exp(tt$truth[tt$covs]), row.names = NULL)
}))
knitr::kable(cox_tab, digits = 2, caption = "Cox (counting-process) HRs agree with the cloglog GLM and the truth.")| transition | covariate | coxph_HR | true_HR |
|---|---|---|---|
| B1: S1 to S2 (Ecological) | temp_c | 1.27 | 1.25 |
| B1: S1 to S2 (Ecological) | humidity_z | 1.18 | 1.16 |
| B1: S1 to S2 (Ecological) | urban_z | 1.25 | 1.28 |
| B2: S2 to S3 (Introduction) | avian_flyway | 1.40 | 1.73 |
| B2: S2 to S3 (Introduction) | trade_z | 1.21 | 1.35 |
| B2: S2 to S3 (Introduction) | nbr_pres | 1.20 | 1.57 |
| B2: S2 to S3 (Introduction) | dist100 | 0.44 | 0.56 |
| B3: S3 to S4 (Transmission) | temp_c | 1.02 | 1.17 |
| B3: S3 to S4 (Transmission) | nbr_pres | 1.06 | 1.28 |
| B3: S3 to S4 (Transmission) | density_z | 1.04 | 1.16 |
| B4: S4 to S5 (Persistence) | temp_c | 1.11 | 1.11 |
| B4: S4 to S5 (Persistence) | socio_z | 1.41 | 1.35 |
statecols <- c("1"="#D9D9D9","2"="#F2B48A","3"="#E8833A","4"="#C0392B","5"="#7A1E15")
statelabs <- c("1"="S1 Non-receptive","2"="S2 Receptivity","3"="S3 Introduction",
"4"="S4 Local transmission","5"="S5 Persistence")h <- chull(lon, lat) # base-R convex hull => "coastline"
hp <- data.frame(lon = lon[h], lat = lat[h])
snap <- do.call(rbind, lapply(c(2008,2014,2020,2024), function(yr)
data.frame(lon = lon, lat = lat, state = factor(spread[[as.character(yr)]]), year = yr)))
ggplot(snap, aes(lon, lat)) +
geom_polygon(data = hp, aes(lon, lat), fill = "#EEF1F5", colour = "#C4CCD8", inherit.aes = FALSE) +
geom_point(aes(colour = state), size = 1.4) +
facet_wrap(~year, nrow = 1) +
scale_colour_manual(values = statecols, labels = statelabs, name = NULL) +
coord_fixed(ratio = 1/cos(mean(lat)*pi/180)) +
theme_void() +
theme(legend.position = "bottom", plot.title = element_text(face="bold", colour="#1F3864")) +
ggtitle("Simulated WNV emergence spreading from the SE origin across Europe")
occ <- do.call(rbind, lapply(2006:2024, function(yr) {
st <- spread[[as.character(yr)]]
data.frame(year = yr, state = factor(1:5), frac = vapply(1:5, function(s) mean(st==s), numeric(1)))
}))
ggplot(occ, aes(year, frac, fill = state)) + geom_area() +
scale_fill_manual(values = statecols, labels = statelabs, name = NULL) +
labs(x = "Year", y = "Fraction of regions",
title = "State-occupation probabilities over calendar time") +
theme_minimal() + theme(plot.title = element_text(face="bold", colour="#1F3864"))
est$label <- factor(est$label, levels = rev(unique(est$label)))
ggplot(est, aes(est_HR, label)) +
geom_vline(xintercept = 1, linetype = "dotted", colour = "grey60") +
geom_errorbarh(aes(xmin = lo, xmax = hi), height = .2, colour = "#7A8CA8") +
geom_point(aes(colour = "Estimated (95% CI)"), size = 2.4) +
geom_point(aes(x = true_HR, colour = "True (simulated)"), shape = 18, size = 3.6) +
facet_wrap(~transition, scales = "free") +
scale_colour_manual(values = c("Estimated (95% CI)"="#1F3864","True (simulated)"="#C0392B"), name = NULL) +
labs(x = "Hazard ratio", y = NULL,
title = "Multi-state estimates recover the known transition effects") +
theme_bw() + theme(plot.title = element_text(face="bold", colour="#1F3864"), legend.position = "bottom")
groups <- list(Climate = c("temp_c","humidity_z"),
`Avian / trade / diffusion` = c("avian_flyway","trade_z","nbr_pres","dist100"),
`Socioeconomic / urban` = c("urban_z","density_z","socio_z"))
barrier <- do.call(rbind, lapply(names(trans), function(tn) {
tt <- trans[[tn]]; sub <- panel[panel$state_start == tt$origin, ]; e <- est[est$transition == tn, ]
do.call(rbind, lapply(names(groups), function(g) {
cs <- intersect(groups[[g]], tt$covs); tot <- 0
for (c in cs) tot <- tot + abs(e$est_beta[e$covariate == c]) * sd(sub[[c]])
data.frame(transition = tn, driver = g, effect = tot)
}))
}))
ggplot(barrier, aes(transition, effect, fill = driver)) +
geom_col(position = "dodge") +
scale_fill_manual(values = c("Climate"="#E0791F","Avian / trade / diffusion"="#4E9B5B",
"Socioeconomic / urban"="#8A6FB0")) +
labs(x = NULL, y = expression("Standardised effect " * "|" * beta * "|" %.% "SD"),
title = "Which barrier is climate-limited vs diffusion-limited?") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 12, hjust = 1), plot.title = element_text(face="bold", colour="#1F3864"))
mstate transition probabilities (optional)library(mstate)
tmat <- transMat(x = list(c(2), c(3), c(4), c(5), c()),
names = c("S1","S2","S3","S4","S5"))
# msprep(...) region-level entry/exit times per state, then:
# ms <- msfit(cox_full, newdata, trans = tmat); pt <- probtrans(ms, predt = 0)
# plot(pt, from = 1)