Getting started with npfseir

Overview

npfseir implements the online Bayesian inference framework of Feng and Wang (2025) for stochastic SEIR epidemic models with a time-varying transmission rate. The transmission rate \(\beta_t\) is modelled as a latent Ornstein–Uhlenbeck (OU) process in log-space, and inference is performed via the nested particle filter (NPF) of Crisan and Miguez (2018).

1. Simulate a synthetic epidemic

fixed <- list(
  eps        = 1/5,        # incubation rate (5-day mean)
  gamma      = 1/10,       # recovery rate (10-day mean)
  delta      = 1/(70*365), # background mortality
  b          = 1/(70*365), # birth rate
  q          = 0.30,       # case detection probability
  N          = 1e6,        # reference population size
  sigma_comp = 0.01        # compartment diffusion coefficient
)

x0   <- c(S = 1e6 - 100, E = 50, I = 50, R = 0, Psi = log(0.25))
traj <- seir_simulate(n_steps = 300, kappa = 0.05, sigma_psi = 0.10,
                      mu_psi = log(0.25), x0 = x0, fixed = fixed)

plot(traj$day, traj$obs, type = "l", col = "gray40",
     xlab = "Day", ylab = "Daily confirmed cases",
     main = "Simulated epidemic (300 days)")
lines(traj$day, exp(traj$Psi) * 3000, col = "steelblue", lty = 2)
legend("topright", c("Observed counts", "True beta_t (scaled)"),
       col = c("gray40","steelblue"), lty = c(1,2), bty = "n")

2. Run the nested particle filter

The npf_seir() function takes only the observation vector and the fixed parameter list. By default it uses \(K = 100\) outer particles and \(M = 200\) inner particles; we use smaller values here for speed.

obs <- traj$obs[-1]   # observations P_1, ..., P_T (exclude day-0 seed)

fit <- npf_seir(
  observations = obs,
  fixed        = fixed,
  K            = 50,    # outer (parameter) particles
  M            = 100,   # inner (state) particles
  seed         = 1
)

3. Inspect results

print(fit)
#> Nested Particle Filter - stochastic SEIR with latent OU transmission
#>   T = 300 observations   K = 50 outer   M = 100 inner particles
#>   Posterior means of OU hyperparameters (at final step):
#>     kappa     = 0.1722  (timescale ~5.8 days)
#>     sigma_Psi = 0.1811
#>     mu_Psi    = -1.3324   exp(mu_Psi) = 0.2638
summary(fit, burn = 30)
#> === npf_seir summary ===
#> 
#> Call:
#> npf_seir(observations = obs, fixed = fixed, K = 50, M = 100, 
#>     seed = 1)
#> 
#> Posterior parameter estimates (weighted mean at final step):
#>   kappa     = 0.1722  [alpha = 0.8418]
#>   sigma_Psi = 0.1811  [tau   = 0.1666]
#>   mu_Psi    = -1.3324  [exp(mu_Psi) = 0.2638]
#>   Implied reproduction scale: exp(mu_Psi)/gamma = 2.638
#> 
#> R_t range (steps 31-300): min=0.326  median=0.612  max=3.684
#> Days with R_t > 1: 83 / 270

4. Plot the R_t trajectory

plot(fit, burn = 30)
abline(h = 1, lty = 2, col = "red")

5. Forecast

fc <- predict(fit, horizon = 14, n_draws = 500)
cat("14-day ahead forecast (mean ± 95% PI):\n")
#> 14-day ahead forecast (mean <U+00B1> 95% PI):
print(round(data.frame(day = 1:14, mean = fc$mean,
                       lo = fc$lo, hi = fc$hi), 1))
#>    day mean  lo   hi
#> 1    1 14.2 7.0 23.0
#> 2    2 13.3 6.0 22.0
#> 3    3 13.2 7.0 22.5
#> 4    4 12.6 5.0 21.0
#> 5    5 12.0 5.0 21.0
#> 6    6 11.5 4.0 20.0
#> 7    7 11.2 4.5 20.5
#> 8    8 10.6 4.0 19.0
#> 9    9 10.5 4.0 19.5
#> 10  10 10.0 4.0 18.0
#> 11  11  9.7 3.0 19.0
#> 12  12  9.5 3.0 18.0
#> 13  13  9.2 3.0 18.0
#> 14  14  8.5 2.0 17.0

6. Cori-style R_t comparison (illustrative only)

cori_rt() implements an in-house Cori-style renewal estimator for illustrative comparison of estimands only. It is not the R EpiEstim package.

ct <- cori_rt(obs, tau = 7, si_mean = 5.5, si_sd = 2.5)
plot(ct$t, ct$Rt_mean, type = "l", col = "darkred",
     ylim = c(0, 4), xlab = "Day", ylab = expression(R[t]),
     main = "Cori-style Rt (illustrative only)")
abline(h = 1, lty = 2)

References

Crisan, D. and Miguez, J. (2018). Nested particle filters for online parameter estimation in discrete-time state-space Markov models. Bernoulli, 24(4A):3039–3086.

Feng, W. and Wang, W. (2025). Bayesian sequential inference for a stochastic SEIR model with latent transmission dynamics. Preprint.