Getting started with npfseir

The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.

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.

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.