EpiNova

# Install from GitHub
devtools::install_github("causalfragility-lab/EpiNova")

# Core dependencies (installed automatically)
# deSolve, ggplot2, dplyr, tidyr, splines, scales

# Optional — unlock extra capabilities
install.packages("DEoptim")   # faster global optimisation in fit_mle()
install.packages("EpiEstim")  # full Bayesian Rt estimation
install.packages("MASS")      # parametric bootstrap uncertainty
install.packages("numDeriv")  # Hessian-based confidence intervals
install.packages("cmdstanr")  # HMC backend (Stan)

library(EpiNova)

# --- Built-in Hubei COVID-19 data ---
props <- prep_proportions(hubei_covid)   # Y = active infected, R = removed

# --- Smooth spline intervention ---
pi_fn <- build_pi_spline(
  knot_times  = c(0, 10, 22, 60),
  knot_values = c(1, 0.9, 0.4, 0.25)
)

# --- Solve SEIRD ODE ---
traj <- solve_model(
  params = list(beta = 0.35, gamma = 0.07, sigma = 0.2, delta = 0.005),
  init   = c(S = 0.9999, E = 0.00005, I = 0.00005, R = 0, D = 0),
  times  = 0:200,
  model  = "SEIRD",
  pi_fn  = pi_fn
)
plot_trajectory(traj, obs_Y = props$Y, obs_R = props$R, T_obs_end = 29)

# --- Scenario comparison ---
sc <- do.call(rbind, lapply(
  list("No NPI"   = function(t) 1,
       "Lockdown" = build_pi_step(c(10), c(1, 0.4)),
       "Spline"   = pi_fn),
  function(pi) {
    tr <- solve_model(
      list(beta=0.35, gamma=0.07, sigma=0.2, delta=0.005),
      c(S=0.9999, E=5e-5, I=5e-5, R=0, D=0), 0:200, "SEIRD", pi_fn=pi)
    data.frame(time=tr$time, I_median=tr$I,
               I_lower=tr$I*0.8, I_upper=tr$I*1.2, scenario=deparse(substitute(pi)))
  }
))
plot_scenarios(sc)

# --- Real-time Rt estimation (no extra packages) ---
new_cases <- pmax(0L, diff(hubei_covid$NI))
Rt_df     <- estimate_Rt_simple(new_cases, mean_si = 5.2, window = 7L)
plot_Rt(Rt_df, change_times = c(10, 22))

# Fast MLE — no external binary needed
fit <- fit_mle(Y, R, N, model = "SEIRD",
               par_init = list(beta=0.3, gamma=0.1, sigma=0.2,
                               delta=0.003, I0=1e-4, E0=2e-4))

# Real-time Sequential Monte Carlo
fit <- fit_smc(Y, R, N, model = "SEIR",
               prior_fn = my_prior, n_particles = 2000)

# Dependency-free Rt estimation
Rt <- estimate_Rt_simple(new_cases, mean_si = 5.2, window = 7L)

# Full Bayesian Rt via EpiEstim (when installed)
Rt <- estimate_Rt(new_cases, mean_si = 5.2, sd_si = 2.8)

# Build mobility matrix from geography
M <- gravity_mobility(N_vec, dist_mat, kappa = 1e-7)

# Couple n patches with per-patch parameters and interventions
ode <- build_multipatch_SEIR(
  n_patches  = 3,
  M          = M,
  beta_vec   = c(0.35, 0.28, 0.22),
  pi_fn_list = list(strict_lockdown, mild_lockdown, function(t) 1)
)
mp_df <- solve_multipatch(ode, init_mat, times = 0:150, n_patches = 3)
plot_multipatch_snapshot(mp_df, t_snapshot = 30,
                          patch_names = c("Hubei", "Guangdong", "Beijing"))

# BMA across model types
ens <- ensemble_forecast(Y, R, N,
         models    = c("SEIR", "SEIRD", "SVEIRD"),
         par_inits = list(...), par_bounds = list(...))

# "What-if" scenario projection
sc <- project_scenarios(fit, scenarios = list(
  "No intervention" = function(t) 1,
  "Light lockdown"  = build_pi_step(c(10), c(1, 0.6)),
  "Strict lockdown" = build_pi_spline(c(0,10,20,60), c(1,0.8,0.2,0.2))
))

# Score your forecasts with proper scoring rules
score_forecast(ens$ensemble_forecast, actual_Y = Y_holdout)
# Returns: CRPS, 95% coverage, MAE

plot_trajectory(traj, obs_Y = Y, obs_R = R, T_obs_end = 29)
plot_forecast(forecast_df, obs_Y = Y)
plot_scenarios(sc_df)
plot_Rt(Rt_df, change_times = c(10, 22))
plot_multipatch_snapshot(mp_df, t_snapshot = 30)

@Manual{EpiNova2025,
  title  = {EpiNova: Flexible Extended State-Space Epidemiological Models
             with Modern Inference},
  author = {Hait, Subir},
  year   = {2026},
  note   = {R package version 0.1.0},
  url    = {https://github.com/causalfragility-lab/EpiNova}
}