Epidemic Forecasting: Methods, Evaluation, and AI

Epidemic forecasting estimates future disease incidence, hospital demand, or other outbreak outcomes over a defined horizon. Mechanistic, statistical, machine-learning, and ensemble models answer different questions, but every forecast should report its target, horizon, uncertainty, data cutoff, and assumptions.

Epidemic Forecasting Defined

Epidemic forecasting is the probabilistic estimation of future outbreak outcomes from surveillance data and explicit model assumptions. It differs from outbreak detection, which asks whether transmission is occurring now, and from scenario modeling, which asks what could happen under specified conditions. For operational use, forecasts should be judged by calibration and proper scoring rules, not by a single point estimate.

Learning Objectives

This chapter distinguishes outbreak detection from epidemic forecasting. You will learn to:

  • Understand why epidemic forecasting is harder than weather forecasting
  • Implement mechanistic compartmental models (SIR, SEIR) and their assumptions
  • Apply statistical time series methods (ARIMA) for short-term prediction
  • Build LSTM neural networks for complex epidemic patterns
  • Create ensemble approaches combining multiple forecasting methods
  • Quantify uncertainty with prediction intervals and scenario analyses
  • Evaluate probabilistic predictions using proper scoring rules (CRPS, Brier)
  • Address ethical obligations when predictions influence behavior

Prerequisites: Disease Surveillance and Outbreak Detection recommended, Machine Learning Fundamentals, The Data Problem.

The Big Picture: Epidemic forecasting is fundamentally harder than weather forecasting because humans change behavior in response to predictions (self-fulfilling/self-defeating prophecies), viruses evolve, and small errors cascade into large divergence. COVID-19 exposed these limits: even with unprecedented data, many models degraded during variant waves and behavioral shifts.

Surveillance vs. Forecasting (Critical Distinction):

  • Surveillance: “An outbreak is happening NOW” → based on current data, relatively objective, verifiable immediately
  • Forecasting: “This is what will happen in 2 weeks” → predicting future, depends on assumptions, cannot verify until future arrives
  • Stakes: Get surveillance wrong → miss outbreak. Get forecasting wrong → misallocate resources, lose public trust, fail to prevent catastrophe

Why Forecasting Is So Hard:

  1. Compounding uncertainty: Small parameter errors can produce large trajectory differences over time. A reduction that brings the effective reproduction number below 1 changes growth to decline; changing 2.5 to 2.3 is an 8% reduction and leaves transmission above that threshold.
  2. Self-Defeating Prophecy: Dire forecast → behavior change → forecast does not happen → public loses trust (“models were wrong!”)
  3. Human Adaptation: People change behavior when they see forecasts, breaking model assumptions
  4. Viral Evolution: Omicron produced a distribution shift that degraded models calibrated during Delta circulation
  5. Data Quality Issues: Reporting delays, testing bias, changing definitions (see The Data Problem)

Three Forecasting Approaches:

  1. Mechanistic Models (SIR, SEIR)
    • Based on epidemiological theory (compartments: Susceptible → Exposed → Infected → Recovered)
    • Strengths: Interpretable, captures disease biology, scenario analysis
    • Weaknesses: Strong assumptions (homogeneous mixing, constant parameters), parameter estimation challenges
    • When to use: Understanding transmission dynamics, evaluating intervention scenarios
  2. Statistical Time Series (ARIMA, Prophet)
    • Learn patterns from historical data without biological assumptions
    • Strengths: Good short-term predictions (1–2 weeks), minimal assumptions
    • Weaknesses: No mechanism, fails when patterns change, cannot handle interventions
    • When to use: Short-term operational forecasting with stable patterns
  3. Machine Learning (LSTMs, Gradient Boosting)
    • Neural networks capture complex temporal patterns
    • Strengths: Flexible, handles non-linear relationships, integrates multiple data streams
    • Weaknesses: “Black box,” requires large training data, overfits to specific outbreak contexts
    • When to use: Rich feature sets, complex patterns, adequate training data

Ensemble Methods (A Strong Default, Not a Rule):

Combine models only when their errors and assumptions are sufficiently diverse, then compare the ensemble against individual models and a simple baseline. CDC FluSight demonstrates strong season-wide performance, but its ensemble has also missed rapid rises and peak changes. Equal weighting is a defensible operational default when estimates for learned weights are unstable.

Evaluation (Beyond Accuracy):

  • Proper Scoring Rules: CRPS (Continuous Ranked Probability Score) and the Brier score reward probabilistic calibration, not just point estimates
  • Prediction Intervals: Quantify uncertainty. A 95% interval should contain the observed value 95% of the time
  • Scenario Analysis: “If R₀ drops by 20%, what happens?” A range of possibilities rather than a single prediction
  • Short-term vs. long-term: Forecast error generally increases with horizon; longer-horizon outputs require wider uncertainty and greater caution

COVID-19 Forecasting Lessons:

  • March 2020 public projections were conditional scenarios for a changing policy and behavioral environment, not fixed multi-year forecasts
  • As of September 28, 2021, the U.S. Forecast Hub ensemble reported case forecasts only one week ahead
  • Why models failed: Behavioral responses not captured, viral evolution, testing bias confusing disease growth with testing expansion, parameter assumptions broke
  • Measured horizon effect: In the U.S. Forecast Hub evaluation, the ensemble had relative weighted interval score of 0.61 versus the baseline across 1- to 4-week horizons. Four-week forecasts had about 76% more weighted interval score than 1-week forecasts, indicating worse probabilistic accuracy at longer horizons (Cramer et al., 2022).

Ethical Obligations:

Forecasts influence the outcomes they predict:

  • Transparency: Communicate uncertainty honestly, not false precision
  • Humility: Acknowledge limitations publicly
  • Avoid Crying Wolf: Repeated failed predictions destroy credibility for genuine threats
  • Equity: Ensure forecasts do not worsen disparities (resource allocation based on predictions)

When to Use Which Approach:

  • Mechanistic (SIR/SEIR): Policy scenario analysis, understanding transmission, resource planning
  • Statistical (ARIMA): Short-term operational forecasting (hospital bed planning 1–2 weeks)
  • ML (LSTMs): Complex multi-source data, non-linear patterns, adequate training data
  • Ensemble: When diverse components and prospective evaluation show a benefit over simple baselines

The Takeaway for Public Health Practitioners:

Epidemic forecasting will always be uncertain: quantify the uncertainty and communicate it. Do not chase false precision. Use forecasts for scenario planning (“If this happens, then that”), not as certain predictions. Mechanistic models help explain biology; statistical models work for short-term operations; ML adds flexibility but loses interpretability. Treat ensembles as strong candidates, then verify them against individual models and simple baselines. Most importantly: forecasts change behavior, which changes outcomes; that feedback is the purpose of public health forecasting. Failed forecasts that motivated effective action are successes, not failures.

Introduction: Moving from Description to Prediction

March 2020, Washington, D.C.: The White House Coronavirus Task Force presented conditional epidemic scenarios: 1.5–2.2 million U.S. deaths without mitigation and 100,000–240,000 with mitigation. Officials explicitly described the range as dependent on interventions and evolving evidence (White House briefing, March 31, 2020). Policies, behavior, treatment, vaccination, variants, and the forecast horizon subsequently changed.

January 2023: The United States passed 1.1 million COVID-19 deaths. That cumulative total is not a test of the March 2020 scenarios.

Models struggled with incomplete data, behavioral assumptions proved fragile, and epidemic forecasting faces fundamental limits that no algorithm can overcome.

January 2022: Models calibrated during Delta circulation degraded during Omicron. The virus evolved, immune landscapes changed, behavior shifted, and key assumptions no longer matched the data.

These are modeling failures and also lessons about the fundamental limits of prediction in complex adaptive systems.

The Critical Distinction

Surveillance (Disease Surveillance): “An outbreak is happening NOW”

  • Based on current/recent data
  • Relatively objective
  • Can be verified immediately

Forecasting (this chapter): “Here’s what will happen in 2 weeks”

  • Predicting the future
  • Depends on assumptions
  • Cannot be verified until the future arrives

The difference matters:

For resource allocation:

  • How many hospital beds will we need?
  • Should we order more ventilators?
  • Do we need to open field hospitals?

For policy decisions:

  • Should we implement lockdowns?
  • When can schools safely reopen?
  • Should we mandate masks?

For public communication:

  • Should people cancel travel plans?
  • Is it safe to visit elderly relatives?
  • Do we need to stock up on supplies?

Get surveillance wrong → Miss an outbreak

Get forecasting wrong → Misallocate resources, lose public trust, or fail to prevent catastrophe

Why Forecasting Is Harder Than Detection

The Butterfly Effect in Epidemiology

In chaos theory, the butterfly effect describes how small changes in initial conditions lead to vastly different outcomes. Epidemics exhibit this behavior in extreme form:

Small changes in behavior → Massive changes in outcomes

A reduction that brings the effective reproduction number below 1 changes expected growth to decline. Smaller reductions can still change peak timing and burden, but their magnitude depends on the model, immunity, contact structure, and starting conditions.

The Self-Fulfilling/Self-Defeating Prophecy Problem

Unlike weather forecasting, epidemic forecasts influence the thing being forecasted.

Scenario 1: Self-defeating prophecy

  • Model predicts catastrophic outbreak
  • Public panics, implements strict measures
  • Outbreak is mild
  • Conclusion: “Model was wrong! Overreaction!”

Scenario 2: Self-fulfilling prophecy

  • Model predicts mild outbreak
  • Public relaxes precautions
  • Outbreak explodes
  • Conclusion: “Model was wrong! Should have prepared!”

This creates an impossible situation: Success makes forecasters look wrong.

The Michael Fish Problem

In October 1987, BBC weatherman Michael Fish famously dismissed viewer concerns about an approaching hurricane: “Earlier on today, apparently, a woman rang the BBC and said she heard there was a hurricane on the way. Well, if you’re watching, don’t worry, there isn’t.”

That night, the Great Storm of 1987 hit the UK, the worst storm in 300 years. Eighteen people died. Millions of trees fell. The public never forgot.

Lesson for epidemic forecasting: Reassuring forecasts that turn out wrong destroy trust permanently. Conservative (alarming) forecasts that do not materialize also destroy trust, but more slowly.

The safest strategy: Communicate uncertainty honestly. Never publish a single-number forecast without ranges.

The Horizon Problem

Forecast accuracy degrades rapidly with time:

Forecast Horizon Typical Use
1 week ahead Operational planning
2 weeks ahead Resource allocation
4 weeks ahead Planning with wide uncertainty
3+ months ahead Scenario planning, not forecasting

Accuracy degrades for multiple reasons: unpredictable interventions (lockdowns, vaccine campaigns), behavior changes as people react to case counts, viral evolution introducing new variants, reporting changes from testing policy shifts, weather and seasonality effects, and random superspreading events.

The 2-week wall: Most successful forecasting systems focus on 1–2 week horizons. Anything beyond that is scenario planning, not prediction.

The COVID-19 Forecast Hub evaluation is in Case Study: COVID-19 Forecasting Hub.

The Mechanistic Approach: Compartmental Models

Before machine learning existed, epidemiologists developed mechanistic models, mathematical frameworks based on how diseases actually spread.

These models remain essential today because they are interpretable, biologically grounded, and enable scenario planning.

SIR Model: The Foundation

The SIR model, introduced by Kermack and McKendrick in 1927, divides a population into three compartments:

  • S: Susceptible (can get infected)
  • I: Infected (currently infectious)
  • R: Recovered/Removed (immune or dead)

The flow:

S → I → R
graph LR
 S[S: Susceptible] -->|β: Transmission rate| I[I: Infected]
 I -->|γ: Recovery rate| R[R: Recovered/Removed]

 style S fill:#e1f5dd
 style I fill:#ffe1e1
 style R fill:#e1e8ff
Figure 8.1: SIR compartmental model showing the flow from Susceptible to Infected to Recovered populations. The transmission rate β governs movement from S to I, while the recovery rate γ governs movement from I to R.

Governed by two parameters:

  • β (beta): Transmission rate (contacts per day × probability of transmission)
  • γ (gamma): Recovery rate (1/infectious period)

The basic reproduction number: \[R_0 = \frac{\beta}{\gamma}\]

\(R_0\) is the average number of people one infected person infects in a fully susceptible population.

  • \(R_0 < 1\): Epidemic dies out
  • \(R_0 = 1\): Threshold (each case replaces itself on average)
  • \(R_0 > 1\): Epidemic grows

The mathematics:

\[\frac{dS}{dt} = -\beta \frac{S \cdot I}{N}\]

\[\frac{dI}{dt} = \beta \frac{S \cdot I}{N} - \gamma I\]

\[\frac{dR}{dt} = \gamma I\]

Where \(N = S + I + R\) (total population).

Implementation:

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

def sir_model(y, t, beta, gamma, N):
 """
 Basic SIR model differential equations

 Parameters:
 - y: [S, I, R] current state
 - t: time
 - beta: transmission rate
 - gamma: recovery rate
 - N: total population

 Returns:
 - [dS/dt, dI/dt, dR/dt]
 """
 S, I, R = y

 dS_dt = -beta * S * I / N
 dI_dt = beta * S * I / N - gamma * I
 dR_dt = gamma * I

 return dS_dt, dI_dt, dR_dt

# Illustrative parameters (R₀ = β/γ = 5)
N = 1_000_000 # Total population (e.g., a city)
I0 = 100  # Initial infected
R0 = 0   # Initial recovered (no prior immunity)
S0 = N - I0 - R0 # Initial susceptible

beta = 0.5  # Transmission rate (contacts per day × transmission probability)
gamma = 0.1 # Recovery rate (1/10 day infectious period)

# Time vector (days)
t = np.linspace(0, 200, 200)

# Initial conditions
y0 = S0, I0, R0

# Solve ODE system
solution = odeint(sir_model, y0, t, args=(beta, gamma, N))
S, I, R = solution.T

# Calculate R₀
R0_value = beta / gamma
print(f"Basic Reproduction Number (R₀): {R0_value:.2f}")

# Visualize
fig, axes = plt.subplots(1, 2, figsize=(16, 6))

# Left panel: Epidemic curves
axes[0].plot(t, S, label='Susceptible', color='blue', linewidth=2)
axes[0].plot(t, I, label='Infected', color='red', linewidth=3)
axes[0].plot(t, R, label='Recovered', color='green', linewidth=2)

# Mark peak
peak_idx = np.argmax(I)
axes[0].axvline(t[peak_idx], color='red', linestyle='--', alpha=0.5,
    label=f'Peak: Day {t[peak_idx]:.0f}')
axes[0].axhline(I[peak_idx], color='red', linestyle='--', alpha=0.5)

axes[0].set_xlabel('Days', fontsize=12)
axes[0].set_ylabel('Population', fontsize=12)
axes[0].set_title(f'SIR Model Epidemic Curve (R₀={R0_value:.1f})', fontsize=14)
axes[0].legend(fontsize=11)
axes[0].grid(True, alpha=0.3)

# Right panel: Phase portrait (S vs I)
axes[1].plot(S, I, 'k-', linewidth=2)
axes[1].scatter([S0], [I0], color='red', s=200, zorder=5,
    marker='*', label='Start', edgecolors='black', linewidth=2)
axes[1].scatter([S[-1]], [I[-1]], color='green', s=200, zorder=5,
    marker='*', label='End', edgecolors='black', linewidth=2)

# Threshold line: infections peak when S falls to N / R0
# (equivalently, when the immune fraction reaches 1 - 1/R0)
s_threshold = N / R0_value
axes[1].axvline(s_threshold, color='orange', linestyle='--',
    linewidth=2, label=f'S = N/R₀ (herd immunity {(1 - 1/R0_value)*100:.0f}%)')

axes[1].set_xlabel('Susceptible Population', fontsize=12)
axes[1].set_ylabel('Infected Population', fontsize=12)
axes[1].set_title('Phase Portrait: S vs I', fontsize=14)
axes[1].legend(fontsize=11)
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('sir_basic_model.png', dpi=300)
plt.show()

# Calculate key metrics
peak_infected_idx = np.argmax(I)
peak_day = t[peak_infected_idx]
peak_cases = I[peak_infected_idx]
final_attack_rate = (N - S[-1]) / N * 100
herd_immunity = (1 - 1/R0_value) * 100

print("\n" + "="*60)
print("SIR MODEL RESULTS")
print("="*60)
print(f"\nR₀ (Basic Reproduction Number): {R0_value:.2f}")
print(f"Herd Immunity Threshold: {herd_immunity:.1f}% of population")
print(f"\nPeak Infections: {peak_cases:,.0f} people")
print(f"Peak Day: {peak_day:.0f}")
print(f"\nFinal Attack Rate: {final_attack_rate:.1f}% of population infected")
print(f"Final Susceptible: {S[-1]:,.0f} people ({S[-1]/N*100:.1f}%)")
print(f"Total Infected: {R[-1]:,.0f} people")

# Demonstrate impact of interventions
print("\n" + "="*60)
print("SCENARIO ANALYSIS: Impact of Interventions")
print("="*60)

scenarios = {
 'No intervention (R₀=2.5)': 2.5,
 'Masks + distancing (R₀=2.0)': 2.0,
 'Partial lockdown (R₀=1.5)': 1.5,
 'Full lockdown (R₀=0.9)': 0.9
}

for scenario_name, R0_scenario in scenarios.items():
 beta_scenario = R0_scenario * gamma
 solution_scenario = odeint(sir_model, y0, t, args=(beta_scenario, gamma, N))
 S_s, I_s, R_s = solution_scenario.T

 peak_infections = I_s.max()
 total_infected = R_s[-1]
 attack_rate = total_infected / N * 100

 print(f"\n{scenario_name}")
 print(f" Peak infections: {peak_infections:,.0f}")
 print(f" Total infected: {total_infected:,.0f} ({attack_rate:.1f}%)")
Understanding R₀

R₀ is one of the most important numbers in epidemiology, and it is often misunderstood.

What it is:

  • Number of secondary infections caused by one infected person in a fully susceptible population
  • Determines whether an epidemic grows or dies out
  • Threshold for herd immunity: \(1 - 1/R_0\)

What it is NOT:

  • Not the number of people infected per day (that is incidence)
  • Not constant over time (behavior, interventions, immunity all change it)
  • Not inherent to a pathogen alone (depends on population behavior and contact patterns)

Examples:

  • Measles: \(R_0 \approx\) 12-18 (extremely contagious)
  • Influenza: \(R_0 \approx\) 1.3-1.8 (moderately contagious)
  • SARS-CoV-2 estimates vary by variant, population immunity, behavior, setting, time period, and estimation method. Variant-specific values require variant-specific evidence rather than transfer from a pre-pandemic review (Omicron review).

For a review of the reproduction number’s definition and interpretation, see Delamater et al., 2019, Emerging Infectious Diseases.

SEIR Model: Adding the Latent Period

The SIR model assumes people become infectious immediately after infection. For most diseases, this is wrong: there is an incubation period.

SEIR adds the Exposed compartment:

S → E → I → R
graph LR
 S[S: Susceptible] -->|β: Transmission rate| E[E: Exposed]
 E -->|σ: Incubation rate| I[I: Infected]
 I -->|γ: Recovery rate| R[R: Recovered/Removed]

 style S fill:#e1f5dd
 style E fill:#fff4e1
 style I fill:#ffe1e1
 style R fill:#e1e8ff
Figure 8.2: SEIR compartmental model extending SIR with an Exposed (E) compartment representing the incubation/latent period. Exposed individuals (E) are infected but not yet infectious, making this model more realistic for diseases like COVID-19, influenza, and Ebola that have significant incubation periods.
  • E: Exposed (infected but not yet infectious)

New parameter:

  • σ (sigma): Rate of progression from E to I (1/incubation period)

The mathematics:

\[\frac{dS}{dt} = -\beta \frac{S \cdot I}{N}\]

\[\frac{dE}{dt} = \beta \frac{S \cdot I}{N} - \sigma E\]

\[\frac{dI}{dt} = \sigma E - \gamma I\]

\[\frac{dR}{dt} = \gamma I\]

Why this matters:

For diseases like COVID-19 with long incubation periods (5–6 days), the SEIR model produces more realistic epidemic curves with delayed and broader peaks.

def seir_model(y, t, beta, sigma, gamma, N):
 """
 SEIR model with exposed/latent period

 Parameters:
 - y: [S, E, I, R] current state
 - sigma: rate of progression from E to I (1/incubation period)
 """
 S, E, I, R = y

 dS_dt = -beta * S * I / N
 dE_dt = beta * S * I / N - sigma * E
 dI_dt = sigma * E - gamma * I
 dR_dt = gamma * I

 return dS_dt, dE_dt, dI_dt, dR_dt

# Parameters for COVID-19-like disease
sigma = 1/5.5 # 5.5-day incubation period
gamma = 1/10 # 10-day infectious period
beta = 0.5

E0 = 50 # Initial exposed
I0 = 10 # Initial infected
R0 = 0
S0 = N - E0 - I0 - R0

# Solve SEIR
y0_seir = S0, E0, I0, R0
solution_seir = odeint(seir_model, y0_seir, t, args=(beta, sigma, gamma, N))
S_seir, E_seir, I_seir, R_seir = solution_seir.T

# Solve SIR for comparison (combine E and I into I)
y0_sir = S0, E0 + I0, R0
solution_sir = odeint(sir_model, y0_sir, t, args=(beta, gamma, N))
S_sir, I_sir, R_sir = solution_sir.T

# Compare
fig, axes = plt.subplots(2, 1, figsize=(14, 10))

# Top panel: SEIR compartments
axes[0].plot(t, S_seir, label='Susceptible', color='blue', linewidth=2)
axes[0].plot(t, E_seir, label='Exposed (not yet infectious)', color='orange', linewidth=2)
axes[0].plot(t, I_seir, label='Infected (infectious)', color='red', linewidth=3)
axes[0].plot(t, R_seir, label='Recovered', color='green', linewidth=2)

axes[0].set_ylabel('Population', fontsize=12)
axes[0].set_title('SEIR Model: Four Compartments', fontsize=14)
axes[0].legend(fontsize=11)
axes[0].grid(True, alpha=0.3)

# Bottom panel: Compare infectious populations (SIR vs SEIR)
axes[1].plot(t, I_sir, '--', label='SIR: Infected (no latency)',
   color='pink', linewidth=3, alpha=0.7)
axes[1].plot(t, I_seir, '-', label='SEIR: Infected (with latent period)',
   color='red', linewidth=3)
axes[1].plot(t, E_seir, '-', label='SEIR: Exposed (not infectious)',
   color='orange', linewidth=2, alpha=0.7)

# Mark peaks
sir_peak_day = t[np.argmax(I_sir)]
seir_peak_day = t[np.argmax(I_seir)]
delay = seir_peak_day - sir_peak_day

axes[1].axvline(sir_peak_day, color='pink', linestyle=':', alpha=0.5)
axes[1].axvline(seir_peak_day, color='red', linestyle=':', alpha=0.5)

axes[1].annotate(f'Peak delay: {delay:.0f} days',
    xy=(seir_peak_day, I_seir.max()),
    xytext=(seir_peak_day + 20, I_seir.max() * 0.8),
    arrowprops=dict(arrowstyle='->', color='black', lw=2),
    fontsize=12, bbox=dict(boxstyle='round', facecolor='wheat'))

axes[1].set_xlabel('Days', fontsize=12)
axes[1].set_ylabel('Population', fontsize=12)
axes[1].set_title('SIR vs SEIR: Impact of Latent Period', fontsize=14)
axes[1].legend(fontsize=11)
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('seir_vs_sir.png', dpi=300)
plt.show()

print(f"\nPeak timing difference: {delay:.1f} days")
print(f"SIR peak: {I_sir.max():,.0f} people on day {sir_peak_day:.0f}")
print(f"SEIR peak: {I_seir.max():,.0f} people on day {seir_peak_day:.0f}")
print(f"Peak reduction: {(1 - I_seir.max()/I_sir.max())*100:.1f}%")
Why the Latent Period Matters

The E compartment:

  1. Delays the epidemic peak (by ~1 incubation period)
  2. Reduces peak height slightly (infections spread out more)
  3. Makes interventions less effective if they only target symptomatic individuals (infected people can transmit before symptoms appear, which the basic SEIR model does not capture)
  4. Creates a “pre-symptomatic transmission” problem (critical for COVID-19)

For diseases with long incubation periods (COVID-19, Ebola), SEIR is essential. For diseases with short/negligible incubation (norovirus), SIR suffices.

Advanced Compartmental Models

Real epidemic modeling often requires more complexity:

SEIRD: Adding Death

\[S \to E \to I \to \begin{cases} R \text{ (recovered)} \\ D \text{ (dead)} \end{cases}\]

  • Split \(I \to R\) into \(I \to R\) (probability \(1-\text{CFR}\)) and \(I \to D\) (probability \(\text{CFR}\))
  • Essential for forecasting mortality burden

SEIRS: Waning Immunity

\[S \to E \to I \to R \to S\]

  • Immunity does not last forever (influenza, coronaviruses)
  • Creates endemic equilibrium (disease does not disappear)

Age-Structured Models

Different parameters for different age groups:

# Children (0-17)
beta_children = 0.3
gamma_children = 0.15

# Adults (18-64)
beta_adults = 0.5
gamma_adults = 0.1

# Elderly (65+)
beta_elderly = 0.4
gamma_elderly = 0.08
cfr_elderly = 0.15 # Higher case fatality rate

Critical for COVID-19 where age dramatically affects outcomes.

Spatial Models

Divide population into geographic regions with mobility between them:

# Movement between regions
mobility_matrix = [
 [0.95, 0.03, 0.02], # Region 1: 95% stay, 3% to R2, 2% to R3
 [0.02, 0.96, 0.02], # Region 2
 [0.01, 0.02, 0.97] # Region 3
]

For more on spatial models, see Keeling & Rohani, 2008, Modeling Infectious Diseases.

A 2026 Nature Medicine gravity-kernel analysis makes the same spatial-scale point with US school-entry MMR records rather than a three-region toy matrix. Chen and Bento assembled kindergarten coverage for more than 50,000 schools, 13,000 districts, and 3,000 counties (45 states plus DC, 2013–2025) and found that school-level measles Rᵥ crossed 1 in 2022–2023 while county-level averages stayed below threshold, as susceptibility roughly doubled from about 5% to 10% after the pandemic (Chen and Bento, 2026). That is modeled transmission potential under R₀ = 15, not a machine-learning paper, not a deployed school dashboard, and not evidence that finer surveillance reduced cases.

Why Mechanistic Models Matter

Strengths:

  • Interpretable: Every parameter has biological meaning (transmission rate, recovery rate, etc.)
  • Scenario planning: Can simulate interventions (vaccination campaigns, social distancing, school closures)
  • Theory-grounded: Built on a century of epidemiological knowledge
  • Extrapolation: Can predict novel scenarios (new pathogen, untested interventions)
  • Policy-relevant: Policymakers understand “if R₀ drops to X, we get Y cases”

Limitations:

  • Simplistic assumptions: Homogeneous mixing (everyone contacts everyone equally, clearly false)
  • Constant parameters: \(\beta\) and \(\gamma\) change with behavior, seasons, immunity
  • No heterogeneity: Superspreaders, age structure, spatial clustering all ignored in basic models
  • Parameter uncertainty: Do not know true \(\beta\), \(\gamma\) until epidemic is over
  • Forecast degradation: Forecasts are most reliable at 1–2 weeks and become unreliable beyond that
The Fundamental Tension

Mechanistic models are:

  • Transparent and interpretable
  • Based on biological mechanisms
  • But require assumptions that are often violated in reality

Machine learning models (next section) are:

  • Flexible and data-driven
  • Do not require mechanistic assumptions
  • But opaque and hard to extrapolate beyond training data

Neither alone is sufficient. The best forecasting systems use both.

The Machine Learning Approach: Data-Driven Forecasting

Data-driven methods learn patterns directly from surveillance data instead of specifying transmission and recovery rates. They range from classical statistical time series models such as ARIMA, the usual baseline, to machine learning models such as LSTMs.

Time Series Forecasting with ARIMA

ARIMA (AutoRegressive Integrated Moving Average) is a classical statistical method for time series forecasting.

Why ARIMA for epidemics:

  • Captures temporal autocorrelation (yesterday’s cases predict today’s)
  • Handles seasonality (flu peaks in winter, dengue in rainy season)
  • No mechanistic assumptions required
  • Works with limited historical data

ARIMA components:

  1. AR (AutoRegressive): Use past values to predict future \[y_t = c + \phi_1 y_{t-1} + \phi_2 y_{t-2} + ... + \epsilon_t\]

  2. I (Integrated): Differencing to make series stationary \[\Delta y_t = y_t - y_{t-1}\]

  3. MA (Moving Average): Use past errors to improve prediction \[y_t = \mu + \epsilon_t + \theta_1 \epsilon_{t-1} + \theta_2 \epsilon_{t-2} + ...\]

ARIMA(p, d, q):

  • p: Number of autoregressive terms
  • d: Degree of differencing
  • q: Number of moving average terms

For seasonal data: SARIMA adds seasonal components.

Complete implementation for influenza forecasting:

import pandas as pd
import numpy as np
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.seasonal import seasonal_decompose
import matplotlib.pyplot as plt
from sklearn.metrics import mean_absolute_error, mean_squared_error

# Simulate synthetic influenza-like illness data (ILINet-style)
np.random.seed(42)
weeks = pd.date_range('2018-01-01', periods=260, freq='W-MON')
week_of_year = weeks.isocalendar().week

# Seasonal flu pattern (peaks in winter)
baseline = 2.0 + 2.5 * np.exp(-((week_of_year - 6) % 52)**2 / 50)

# Add year-to-year variation
year_effect = np.repeat([0, 0.3, -0.2, 0.5, 0.1], 52)[:len(weeks)]

# Random noise
noise = np.random.normal(0, 0.25, len(weeks))

# Combine
ili_pct = baseline + year_effect + noise
ili_pct = np.maximum(ili_pct, 0.5) # ILI never goes to zero

df = pd.DataFrame({'date': weeks, 'ili_pct': ili_pct})
df.set_index('date', inplace=True)

print(df.head(10))
print(f"\nTotal weeks: {len(df)}")
print(f"Date range: {df.index.min()} to {df.index.max()}")

# Train/test split (reserve last season for testing)
train_size = 208 # 4 years training
train = df.iloc[:train_size]
test = df.iloc[train_size:]

print(f"\nTraining weeks: {len(train)}")
print(f"Test weeks: {len(test)}")

# Exploratory analysis
fig, axes = plt.subplots(3, 1, figsize=(14, 12))

# Raw time series
axes[0].plot(train.index, train['ili_pct'], label='Training Data', color='blue')
axes[0].plot(test.index, test['ili_pct'], label='Test Data', color='orange')
axes[0].set_ylabel('ILI %', fontsize=12)
axes[0].set_title('Influenza-Like Illness Time Series', fontsize=14)
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# ACF (Autocorrelation Function)
plot_acf(train['ili_pct'], lags=52, ax=axes[1], alpha=0.05)
axes[1].set_title('Autocorrelation Function (ACF)', fontsize=14)
axes[1].set_xlabel('Lag (weeks)')

# PACF (Partial Autocorrelation Function)
plot_pacf(train['ili_pct'], lags=52, ax=axes[2], alpha=0.05)
axes[2].set_title('Partial Autocorrelation Function (PACF)', fontsize=14)
axes[2].set_xlabel('Lag (weeks)')

plt.tight_layout()
plt.savefig('arima_exploratory.png', dpi=300)
plt.show()

# Seasonal decomposition
decomposition = seasonal_decompose(train['ili_pct'], model='additive', period=52)

fig, axes = plt.subplots(4, 1, figsize=(14, 12))
decomposition.observed.plot(ax=axes[0], title='Observed')
decomposition.trend.plot(ax=axes[1], title='Trend')
decomposition.seasonal.plot(ax=axes[2], title='Seasonal')
decomposition.resid.plot(ax=axes[3], title='Residual')

for ax in axes:
 ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('seasonal_decomposition.png', dpi=300)
plt.show()

# Fit SARIMA model
# SARIMA(p, d, q)(P, D, Q, s)
# (2, 0, 1) for non-seasonal components
# (1, 0, 1, 52) for seasonal components with 52-week period

print("\nFitting SARIMA model...")
model = SARIMAX(
 train['ili_pct'],
 order=(2, 0, 1),    # (p, d, q)
 seasonal_order=(1, 0, 1, 52), # (P, D, Q, s)
 enforce_stationarity=False,
 enforce_invertibility=False
)

fit = model.fit(disp=False, maxiter=200)

print("\n" + "="*60)
print("MODEL SUMMARY")
print("="*60)
print(fit.summary())

# Forecast
forecast_steps = len(test)
forecast_result = fit.get_forecast(steps=forecast_steps)

# Get predictions and confidence intervals
forecast_mean = forecast_result.predicted_mean
forecast_ci = forecast_result.conf_int()

# Create forecast dataframe
forecast_df = pd.DataFrame({
 'date': test.index,
 'actual': test['ili_pct'].values,
 'forecast': forecast_mean.values,
 'lower_ci': forecast_ci.iloc[:, 0].values,
 'upper_ci': forecast_ci.iloc[:, 1].values
})

# Visualize forecast
fig, axes = plt.subplots(2, 1, figsize=(16, 10))

# Top: Full time series with forecast
axes[0].plot(train.index, train['ili_pct'], 'o-',
   label='Training Data', color='blue', alpha=0.6, markersize=4)
axes[0].plot(test.index, test['ili_pct'], 'o-',
   label='Actual (Test)', color='black', linewidth=2, markersize=6)
axes[0].plot(forecast_df['date'], forecast_df['forecast'], 's-',
   label='SARIMA Forecast', color='red', linewidth=2, markersize=6)

# Confidence interval
axes[0].fill_between(forecast_df['date'],
      forecast_df['lower_ci'],
      forecast_df['upper_ci'],
      alpha=0.3, color='red', label='95% Confidence Interval')

axes[0].axvline(train.index[-1], color='green', linestyle='--',
    linewidth=2, label='Train/Test Split')

axes[0].set_ylabel('ILI %', fontsize=12)
axes[0].set_title('SARIMA Forecast: Influenza-Like Illness', fontsize=14)
axes[0].legend(fontsize=10)
axes[0].grid(True, alpha=0.3)

# Bottom: Forecast errors
errors = forecast_df['actual'] - forecast_df['forecast']
colors = ['red' if e < 0 else 'green' for e in errors]

axes[1].bar(forecast_df['date'], errors, color=colors, alpha=0.6, edgecolor='black')
axes[1].axhline(0, color='black', linestyle='-', linewidth=1.5)

# Mark forecast horizons
for horizon, label in [(4, '1 month'), (13, '3 months')]:
 if horizon < len(errors):
  axes[1].axvline(forecast_df['date'].iloc[horizon],
      color='blue', linestyle=':', linewidth=2)
  axes[1].text(forecast_df['date'].iloc[horizon], errors.max() * 0.9,
     f'{label} ahead', rotation=90, fontsize=10)

axes[1].set_xlabel('Date', fontsize=12)
axes[1].set_ylabel('Forecast Error (Actual - Predicted)', fontsize=12)
axes[1].set_title('Forecast Errors Over Time', fontsize=14)
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('arima_flu_forecast.png', dpi=300)
plt.show()

# Evaluate performance at different horizons
horizons = [1, 4, 8, 13] # 1 week, 1 month, 2 months, 3 months

print("\n" + "="*60)
print("FORECAST PERFORMANCE BY HORIZON")
print("="*60)

for horizon in horizons:
 if horizon <= len(errors):
  subset_actual = forecast_df['actual'].iloc[:horizon]
  subset_pred = forecast_df['forecast'].iloc[:horizon]

  mae = mean_absolute_error(subset_actual, subset_pred)
  rmse = np.sqrt(mean_squared_error(subset_actual, subset_pred))
  mape = np.mean(np.abs((subset_actual - subset_pred) / subset_actual)) * 100

  # Coverage: % of actuals within CI
  subset_lower = forecast_df['lower_ci'].iloc[:horizon]
  subset_upper = forecast_df['upper_ci'].iloc[:horizon]
  coverage = ((subset_actual >= subset_lower) & (subset_actual <= subset_upper)).mean() * 100

  print(f"\n{horizon} week(s) ahead:")
  print(f" MAE: {mae:.3f}")
  print(f" RMSE: {rmse:.3f}")
  print(f" MAPE: {mape:.1f}%")
  print(f" 95% CI Coverage: {coverage:.0f}%")

# Compare to naive baseline (persistence forecast)
naive_forecast = train['ili_pct'].iloc[-1] # Last observed value
naive_errors = forecast_df['actual'] - naive_forecast
naive_mae = mean_absolute_error(forecast_df['actual'], [naive_forecast]*len(forecast_df))
arima_mae = mean_absolute_error(forecast_df['actual'], forecast_df['forecast'])

improvement = (1 - arima_mae / naive_mae) * 100

print("\n" + "="*60)
print("COMPARISON TO NAIVE BASELINE")
print("="*60)
print(f"Naive forecast (persistence): MAE = {naive_mae:.3f}")
print(f"SARIMA forecast: MAE = {arima_mae:.3f}")
print(f"Improvement: {improvement:.1f}%")
When ARIMA Works (and When It Does not)

ARIMA excels when:

  • Disease has regular seasonal patterns (influenza, RSV, dengue)
  • Reporting is stable and consistent
  • Short-term forecasts (1–4 weeks)
  • No major interventions or behavior changes expected
  • Historical data covers multiple seasons

ARIMA struggles when:

  • Novel outbreak (no historical pattern to learn)
  • Sudden interventions (lockdowns, vaccine rollouts)
  • Reporting systems change
  • Long-term forecasts (>4 weeks)
  • Nonlinear, chaotic dynamics

For comprehensive ARIMA tutorial, see Hyndman & Athanasopoulos, Forecasting: Principles and Practice.

A 2026 BMC Infectious Diseases retrospective analysis compared climate-informed Random Forest models with SHAP interpretation to a climate-naive SARIMA benchmark on a shared 52-week temporal holdout, using weekly Lassa fever surveillance from six endemic Nigerian states (2018–2025; n = 6,801 cases) and satellite-derived covariates (soil moisture, relative humidity, rainfall, temperature, and NDVI) (Tadema et al., 2026). Random Forest outperformed SARIMA in Edo, Ondo, and Plateau; SARIMA performed better in Bauchi, Ebonyi, and Taraba: a 3/3 geographic split, not a uniform “machine learning beats SARIMA” result, and the authors support locally tailored climate-informed surveillance rather than a single national forecasting approach. Soil moisture was the dominant SHAP climate driver in five of six states (rainfall in Bauchi), and in Ebonyi higher soil moisture was associated with lower weekly incidence (Negative Binomial IRR 0.836, 95% CI 0.797–0.876), an observational association rather than evidence of causal intervention. Retrospective model comparison on surveillance data; not a deployed early-warning system, not prospective alert evaluation, and the paper’s “early warning framework” title does not establish a live operational EWS.

Deep Learning: LSTMs for Epidemic Forecasting

LSTM (Long Short-Term Memory) networks are a type of recurrent neural network designed for sequential data.

Why LSTMs for epidemics:

  • Can capture long-term dependencies (last month’s cases affect this month)
  • Handle multiple input features (cases + mobility + weather + interventions)
  • Do not require stationarity assumptions
  • Can learn nonlinear patterns

The basic architecture:

Input sequence → LSTM layers → Dense layers → Output (future cases)
 (past 14 days)        (next 7 days)

Complete implementation for COVID-19 forecasting:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.preprocessing import MinMaxScaler
from sklearn.metrics import mean_absolute_error, mean_squared_error
import tensorflow as tf
from tensorflow import keras
from tensorflow.keras import layers
import warnings
warnings.filterwarnings('ignore')

# Set random seeds for reproducibility
np.random.seed(42)
tf.random.set_seed(42)

# Simulate synthetic COVID-19 case data
days = np.arange(500)

# Multi-wave pattern (Original → Delta → Omicron-like)
wave1 = 500 * np.exp(-((days - 50)**2) / 300)
wave2 = 2000 * np.exp(-((days - 200)**2) / 400)
wave3 = 5000 * np.exp(-((days - 350)**2) / 500)

# Combine waves with baseline and noise
cases = 100 + wave1 + wave2 + wave3 + np.random.normal(0, 100, len(days))
cases = np.maximum(cases, 50) # Floor at 50 cases/day

# Add day-of-week effect (reporting artifacts)
day_of_week_effect = -200 * (days % 7 >= 5) # Lower on weekends
cases = cases + day_of_week_effect
cases = np.maximum(cases, 0)

# Create DataFrame
df = pd.DataFrame({
 'day': days,
 'cases': cases
})

print(f"Total days: {len(df)}")
print(f"Case range: {cases.min():.0f} to {cases.max():.0f}")

# Visualize raw data
plt.figure(figsize=(14, 6))
plt.plot(df['day'], df['cases'], 'o-', alpha=0.6, markersize=3)
plt.xlabel('Day', fontsize=12)
plt.ylabel('Daily Cases', fontsize=12)
plt.title('Simulated COVID-19 Daily Cases (Multi-Wave Pattern)', fontsize=14)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('lstm_data_raw.png', dpi=300)
plt.show()

# Prepare data for LSTM
def create_sequences(data, seq_length, forecast_horizon=1):
 """
 Create input sequences and targets for LSTM

 Parameters:
 - data: 1D array of values
 - seq_length: number of past time steps to use
 - forecast_horizon: number of future time steps to predict

 Returns:
 - X: [n_samples, seq_length] input sequences
 - y: [n_samples, forecast_horizon] target values
 """
 X, y = [], []

 for i in range(len(data) - seq_length - forecast_horizon + 1):
  # Input: seq_length days of history
  X.append(data[i:i+seq_length])

  # Target: next forecast_horizon days
  if forecast_horizon == 1:
   y.append(data[i+seq_length])
  else:
   y.append(data[i+seq_length:i+seq_length+forecast_horizon])

 return np.array(X), np.array(y)

# Normalize data (LSTM trains better with normalized inputs)
scaler = MinMaxScaler(feature_range=(0, 1))
cases_scaled = scaler.fit_transform(df['cases'].values.reshape(-1, 1)).flatten()

# Create sequences (use past 14 days to predict next day)
seq_length = 14
forecast_horizon = 1

X, y = create_sequences(cases_scaled, seq_length, forecast_horizon)

print(f"\nSequence creation:")
print(f"Input shape: {X.shape}") # (n_samples, seq_length)
print(f"Target shape: {y.shape}") # (n_samples,)

# Train/test split (80/20)
train_size = int(0.8 * len(X))
X_train, X_test = X[:train_size], X[train_size:]
y_train, y_test = y[:train_size], y[train_size:]

print(f"\nTrain samples: {len(X_train)}")
print(f"Test samples: {len(X_test)}")

# Reshape for LSTM [samples, time steps, features]
X_train = X_train.reshape((X_train.shape[0], X_train.shape[1], 1))
X_test = X_test.reshape((X_test.shape[0], X_test.shape[1], 1))

print(f"\nReshaped for LSTM:")
print(f"X_train shape: {X_train.shape}")
print(f"X_test shape: {X_test.shape}")

# Build LSTM model
model = keras.Sequential([
 # First LSTM layer (return sequences for stacking)
 layers.LSTM(64, activation='relu', return_sequences=True,
    input_shape=(seq_length, 1)),
 layers.Dropout(0.2),

 # Second LSTM layer
 layers.LSTM(32, activation='relu'),
 layers.Dropout(0.2),

 # Dense layers
 layers.Dense(16, activation='relu'),
 layers.Dense(forecast_horizon) # Output layer
])

model.compile(
 optimizer='adam',
 loss='mse',
 metrics=['mae']
)

print("\n" + "="*60)
print("LSTM MODEL ARCHITECTURE")
print("="*60)
model.summary()

# Train model
print("\n" + "="*60)
print("TRAINING MODEL")
print("="*60)

early_stopping = keras.callbacks.EarlyStopping(
 monitor='val_loss',
 patience=10,
 restore_best_weights=True
)

history = model.fit(
 X_train, y_train,
 epochs=100,
 batch_size=32,
 validation_split=0.2,
 callbacks=[early_stopping],
 verbose=1
)

# Make predictions
y_train_pred = model.predict(X_train, verbose=0)
y_test_pred = model.predict(X_test, verbose=0)

# Inverse transform to original scale
y_train_actual = scaler.inverse_transform(y_train.reshape(-1, 1)).flatten()
y_train_pred_original = scaler.inverse_transform(y_train_pred).flatten()

y_test_actual = scaler.inverse_transform(y_test.reshape(-1, 1)).flatten()
y_test_pred_original = scaler.inverse_transform(y_test_pred).flatten()

# Evaluate
train_mae = mean_absolute_error(y_train_actual, y_train_pred_original)
train_rmse = np.sqrt(mean_squared_error(y_train_actual, y_train_pred_original))

test_mae = mean_absolute_error(y_test_actual, y_test_pred_original)
test_rmse = np.sqrt(mean_squared_error(y_test_actual, y_test_pred_original))

print("\n" + "="*60)
print("MODEL PERFORMANCE")
print("="*60)
print(f"\nTraining Set:")
print(f" MAE: {train_mae:.1f} cases")
print(f" RMSE: {train_rmse:.1f} cases")

print(f"\nTest Set:")
print(f" MAE: {test_mae:.1f} cases")
print(f" RMSE: {test_rmse:.1f} cases")

# Visualize results
fig, axes = plt.subplots(3, 1, figsize=(16, 14))

# Top: Training history
axes[0].plot(history.history['loss'], label='Training Loss', linewidth=2)
axes[0].plot(history.history['val_loss'], label='Validation Loss', linewidth=2)
axes[0].set_xlabel('Epoch', fontsize=12)
axes[0].set_ylabel('Loss (MSE)', fontsize=12)
axes[0].set_title('LSTM Training History', fontsize=14)
axes[0].legend(fontsize=11)
axes[0].grid(True, alpha=0.3)

# Middle: Predictions vs actual (full dataset)
train_days = days[seq_length:seq_length+len(y_train_actual)]
test_days = days[seq_length+len(y_train_actual):seq_length+len(y_train_actual)+len(y_test_actual)]

axes[1].plot(train_days, y_train_actual, 'o-',
   label='Actual (Training)', color='blue', alpha=0.5, markersize=3)
axes[1].plot(test_days, y_test_actual, 'o-',
   label='Actual (Test)', color='black', linewidth=2, markersize=4)
axes[1].plot(test_days, y_test_pred_original, 's-',
   label='LSTM Forecast', color='red', linewidth=2, markersize=4)

axes[1].axvline(train_days[-1], color='green', linestyle='--',
    linewidth=2, label='Train/Test Split')

axes[1].set_xlabel('Day', fontsize=12)
axes[1].set_ylabel('Daily Cases', fontsize=12)
axes[1].set_title('LSTM Forecast: COVID-19 Daily Cases', fontsize=14)
axes[1].legend(fontsize=10)
axes[1].grid(True, alpha=0.3)

# Bottom: Forecast errors
errors = y_test_actual - y_test_pred_original
colors = ['red' if e < 0 else 'green' for e in errors]

axes[2].bar(test_days, errors, color=colors, alpha=0.6, edgecolor='black')
axes[2].axhline(0, color='black', linestyle='-', linewidth=1.5)

axes[2].set_xlabel('Day', fontsize=12)
axes[2].set_ylabel('Forecast Error (Actual - Predicted)', fontsize=12)
axes[2].set_title('LSTM Forecast Errors', fontsize=14)
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('lstm_covid_forecast.png', dpi=300)
plt.show()

# Multi-step ahead forecasting
print("\n" + "="*60)
print("MULTI-STEP AHEAD FORECASTING")
print("="*60)

# Use the last seq_length days as input
last_sequence = cases_scaled[-seq_length:]
current_sequence = last_sequence.reshape(1, seq_length, 1)

# Forecast next 14 days iteratively
forecast_days = 14
multi_step_forecast = []

for i in range(forecast_days):
 # Predict next day
 next_pred = model.predict(current_sequence, verbose=0)[0, 0]
 multi_step_forecast.append(next_pred)

 # Update sequence: drop oldest, add prediction
 current_sequence = np.append(current_sequence[0, 1:, 0], next_pred).reshape(1, seq_length, 1)

# Inverse transform
multi_step_forecast_original = scaler.inverse_transform(
 np.array(multi_step_forecast).reshape(-1, 1)
).flatten()

# Visualize multi-step forecast
fig, ax = plt.subplots(figsize=(14, 6))

# Historical data
ax.plot(days, cases, 'o-', label='Historical Cases',
  color='blue', alpha=0.6, markersize=3)

# Multi-step forecast
forecast_days_array = np.arange(len(days), len(days) + forecast_days)
ax.plot(forecast_days_array, multi_step_forecast_original, 's-',
  label=f'{forecast_days}-Day Forecast', color='red',
  linewidth=3, markersize=8)

ax.axvline(days[-1], color='green', linestyle='--', linewidth=2,
   label='Forecast Start')

ax.set_xlabel('Day', fontsize=12)
ax.set_ylabel('Daily Cases', fontsize=12)
ax.set_title(f'LSTM Multi-Step Ahead Forecast ({forecast_days} days)', fontsize=14)
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('lstm_multistep_forecast.png', dpi=300)
plt.show()

print(f"\n{forecast_days}-day forecast:")
for i, pred in enumerate(multi_step_forecast_original, 1):
 print(f" Day {i}: {pred:.0f} cases")
The LSTM Limitation: No Mechanistic Understanding

LSTMs learn patterns, not mechanisms.

Example scenario:

  • Training data includes natural epidemic curve + lockdown → cases drop
  • LSTM learns: “After cases reach X, they drop”
  • New epidemic: Cases reach X, but no lockdown this time
  • LSTM incorrectly predicts a drop

The problem: LSTM has no concept of “lockdown caused the drop.” It just sees a pattern.

When this matters:

  • Novel interventions (no historical data)
  • Changing policy responses
  • New variants with different characteristics
  • Long-term forecasts where context matters

Solution: Use LSTMs for short-term pattern recognition, mechanistic models for scenario planning.

A third approach constrains neural networks with epidemic structure rather than asking them to learn dynamics from scratch. Fujita and Akutsu combined a spatiotemporal neural network with an SIR module and reported more stable multi-region forecasts than baseline deep learning models while preserving interpretable epidemiologic parameters (Fujita & Akutsu, 2025, PLOS One).

A fourth approach feeds epidemiologic, policy, spatial, and genomic context into large language models rather than relying only on case-count histories. In retrospective U.S. COVID-19 experiments, PandemicLLM reported better 1- to 3-week forecasting performance than several CDC COVID-19 Forecast Hub baselines, especially during periods of rapid change, but evidence remains limited to retrospective evaluation in one disease setting (Du et al., 2025, Nature Computational Science).

For practical deep learning in epidemiology, see Chimmula & Zhang, 2020, Chaos, Solitons & Fractals.

Ensemble Methods: Combining Models

No single model is best for all situations. Ensemble methods combine multiple models to create more robust forecasts.

The CDC FluSight Challenge

The CDC’s FluSight challenge, running since 2013, invites teams to submit weekly influenza forecasts. Beginning with the 2021–2022 season, CDC shifted the primary forecasting target from outpatient influenza-like illness (ILI%) to weekly laboratory-confirmed influenza hospital admissions, reflecting improved data quality and clinical relevance (Mathis et al., 2024).

The setup:

  • Teams submit forecasts for the current week and 1–3 weeks ahead
  • Target: new laboratory-confirmed influenza hospital admissions at state and national levels
  • Evaluated on prediction interval coverage and point forecast accuracy

Key finding (Reich et al., 2019, PNAS):

Across the seasons and targets evaluated in this challenge, the ensemble was consistently competitive and often more accurate than individual submissions.

Recent operational evidence shows why that result should not become an unconditional rule. In CDC’s 2024–2025 FluSight evaluation, the FluSight ensemble had the best season-wide average relative weighted interval score among submitted models. It nevertheless failed to anticipate rapid hospitalization increases and changes around the peak. CDC first published the forecasts on March 5, 2025, after the ensemble’s prediction intervals began to capture future observations reliably.

Weighting also requires an evidence gate. A disease-agnostic learned-weight ensemble outperformed equal weighting on average across the outbreaks evaluated by Murph et al., but the authors found that equal weighting can be the safer real-time choice because learned weights have a higher ceiling and a lower floor (Murph et al., 2026, Nature Communications). Start with equal weighting when historical performance is sparse or unstable. Adopt learned weights only after time-split or prospective evaluation shows durable gains without loss of calibration.

AI-generated forecasting models are now entering the same benchmark arena. In a 2026 Nature study, Aygun et al. introduced Empirical Research Assistance, an LLM and tree-search system that writes scientific software to maximize a prespecified quality metric. In one epidemiology task, it generated 14 COVID-19 hospitalization forecasting models that outperformed the CDC ensemble and all individual models in the benchmark setting (Aygun et al., 2026). For public health agencies, the lesson is not to replace forecast hubs with autonomous code generation. It is that evaluation pipelines must be able to audit algorithmically generated models with the same calibration, uncertainty, and external-validity standards used for human-built submissions.

Why ensembles work:

  1. Different models capture different patterns
    • Mechanistic models: Understand transmission dynamics
    • ARIMA: Captures seasonal patterns
    • ML models: Learn complex nonlinear relationships
  2. Averaging reduces overfitting
    • Individual models may be overconfident
    • Ensemble spreads risk across models
  3. Robustness to model-specific failures
    • If one model fails (e.g., unexpected intervention), others compensate

Building an Ensemble Forecaster

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from statsmodels.tsa.statespace.sarimax import SARIMAX
from sklearn.metrics import mean_absolute_error
import warnings
warnings.filterwarnings('ignore')

# Generate synthetic epidemic data
np.random.seed(42)
days = np.arange(200)

# Simulated cases (mix of exponential growth + intervention)
phase1 = 100 * np.exp(0.03 * days[:80]) # Growth phase
phase2 = phase1[-1] * np.exp(-0.02 * (days[80:] - 80)) # Decline (intervention)

cases = np.concatenate([phase1, phase2])
cases = cases + np.random.normal(0, cases * 0.1) # Add noise
cases = np.maximum(cases, 10)

df = pd.DataFrame({'day': days, 'cases': cases})

# Split: train on first 150 days, forecast next 50
train_size = 150
train = df.iloc[:train_size]
test = df.iloc[train_size:]

print(f"Training days: {len(train)}")
print(f"Test days: {len(test)}")

# === MODEL 1: SEIR Mechanistic Model ===

def seir_forecast(train_data, forecast_steps, beta=0.03, sigma=0.2, gamma=0.1):
 """
 Simulate an SEIR model with fixed (unfitted) parameters and
 return the forecast window
 """
 # Population
 N = 1000000

 # Initial conditions (estimate from data)
 I0 = train_data['cases'].iloc[0]
 E0 = I0
 R0_init = 0
 S0 = N - E0 - I0 - R0_init

 def seir_model(y, t, beta, sigma, gamma, N):
  S, E, I, R = y
  dS = -beta * S * I / N
  dE = beta * S * I / N - sigma * E
  dI = sigma * E - gamma * I
  dR = gamma * I
  return dS, dE, dI, dR

 # Solve SEIR for training + forecast period
 t = np.arange(len(train_data) + forecast_steps)
 y0 = S0, E0, I0, R0_init

 solution = odeint(seir_model, y0, t, args=(beta, sigma, gamma, N))
 _, _, I, _ = solution.T

 # Return only forecast period
 forecast = I[len(train_data):]

 return forecast

seir_forecast_values = seir_forecast(train, len(test))

# === MODEL 2: ARIMA ===

def arima_forecast(train_data, forecast_steps):
 """
 Fit ARIMA model and forecast
 """
 try:
  model = SARIMAX(train_data['cases'], order=(2, 1, 1))
  fit = model.fit(disp=False)
  forecast = fit.forecast(steps=forecast_steps)
  return forecast.values
 except:
  # Fallback: return mean if ARIMA fails
  return np.full(forecast_steps, train_data['cases'].mean())

arima_forecast_values = arima_forecast(train, len(test))

# === MODEL 3: Polynomial Trend Extrapolation (simple baseline) ===

def ml_forecast(train_data, forecast_steps):
 """
 Simple ML: Fit polynomial to recent trend and extrapolate
 """
 from numpy.polynomial import Polynomial

 # Use last 30 days for trend
 recent = train_data.iloc[-30:]

 # Fit polynomial
 p = Polynomial.fit(recent['day'], recent['cases'], deg=2)

 # Forecast
 forecast_days = np.arange(train_data['day'].iloc[-1] + 1,
        train_data['day'].iloc[-1] + 1 + forecast_steps)
 forecast = p(forecast_days)

 return np.maximum(forecast, 0) # No negative cases

ml_forecast_values = ml_forecast(train, len(test))

# === ENSEMBLE: Combine Models ===

# Method 1: Simple average
ensemble_simple = (seir_forecast_values + arima_forecast_values + ml_forecast_values) / 3

# Method 2: Weighted average (based on domain knowledge or historical performance)
weights = {
 'seir': 0.4, # Higher weight for mechanistic model (interpretable)
 'arima': 0.3, # Medium weight (good for short-term)
 'ml': 0.3  # Medium weight (captures recent trends)
}

ensemble_weighted = (
 weights['seir'] * seir_forecast_values +
 weights['arima'] * arima_forecast_values +
 weights['ml'] * ml_forecast_values
)

# Method 3: Median ensemble (more robust to outliers)
all_forecasts = np.array([seir_forecast_values, arima_forecast_values, ml_forecast_values])
ensemble_median = np.median(all_forecasts, axis=0)

# Visualize
fig, axes = plt.subplots(2, 1, figsize=(16, 12))

# Top: All forecasts
axes[0].plot(train['day'], train['cases'], 'o-',
   label='Training Data', color='blue', alpha=0.5, markersize=3)
axes[0].plot(test['day'], test['cases'], 'o-',
   label='Actual (Test)', color='black', linewidth=3, markersize=6)

axes[0].plot(test['day'], seir_forecast_values, '^--',
   label='SEIR Model', alpha=0.7, markersize=6)
axes[0].plot(test['day'], arima_forecast_values, 's--',
   label='ARIMA Model', alpha=0.7, markersize=6)
axes[0].plot(test['day'], ml_forecast_values, 'd--',
   label='ML Model', alpha=0.7, markersize=6)

axes[0].plot(test['day'], ensemble_weighted, '*-',
   label='Weighted Ensemble', color='red', linewidth=3, markersize=10)

axes[0].axvline(train['day'].iloc[-1], color='green', linestyle='--',
    linewidth=2, label='Train/Test Split')

axes[0].set_ylabel('Daily Cases', fontsize=12)
axes[0].set_title('Ensemble Forecasting: Combining Multiple Models', fontsize=14)
axes[0].legend(fontsize=10, loc='upper right')
axes[0].grid(True, alpha=0.3)

# Bottom: Performance comparison
models = {
 'SEIR': seir_forecast_values,
 'ARIMA': arima_forecast_values,
 'ML': ml_forecast_values,
 'Simple Ensemble': ensemble_simple,
 'Weighted Ensemble': ensemble_weighted,
 'Median Ensemble': ensemble_median
}

performance = []
for name, forecast in models.items():
 mae = mean_absolute_error(test['cases'], forecast)
 rmse = np.sqrt(((test['cases'] - forecast) ** 2).mean())
 performance.append({'Model': name, 'MAE': mae, 'RMSE': rmse})

perf_df = pd.DataFrame(performance).sort_values('MAE')

x_pos = np.arange(len(perf_df))
axes[1].barh(x_pos, perf_df['MAE'], alpha=0.8, edgecolor='black')
axes[1].set_yticks(x_pos)
axes[1].set_yticklabels(perf_df['Model'])
axes[1].set_xlabel('Mean Absolute Error (MAE)', fontsize=12)
axes[1].set_title('Model Performance Comparison', fontsize=14)
axes[1].grid(True, alpha=0.3, axis='x')

# Highlight ensemble
for i, model in enumerate(perf_df['Model']):
 if 'Ensemble' in model:
  axes[1].get_children()[i].set_color('red')
  axes[1].get_children()[i].set_alpha(1.0)

plt.tight_layout()
plt.savefig('ensemble_forecast.png', dpi=300)
plt.show()

# Print performance
print("\n" + "="*60)
print("MODEL PERFORMANCE COMPARISON")
print("="*60)
print(perf_df.to_string(index=False))

# Best model
best_model = perf_df.iloc[0]
print(f"\nBest Model: {best_model['Model']}")
print(f" MAE: {best_model['MAE']:.1f}")
print(f" RMSE: {best_model['RMSE']:.1f}")

# Improvement of best ensemble over best individual model
best_individual = perf_df[~perf_df['Model'].str.contains('Ensemble')].iloc[0]
improvement = (1 - best_model['MAE'] / best_individual['MAE']) * 100

print(f"\nEnsemble improvement over best individual model: {improvement:.1f}%")
Ensemble Best Practices

Weight Selection:

  1. Equal weights (simple average)
    • Easiest to implement
    • Works surprisingly well
    • Default starting point
  2. Performance-based weights
    • Weight by historical accuracy
    • Requires validation set
    • Can overfit to past performance
  3. Domain-based weights
    • Mechanistic models: Higher weight for scenario planning
    • ML models: Higher weight for short-term forecasts
    • Time series models: Higher weight for seasonal diseases
  4. Adaptive weights
    • Change weights based on recent performance
    • More complex; learned weights can be less stable in real time
    • Requires careful validation

When ensemble does not help:

  • All models use the same data and make the same mistakes
  • One model is vastly superior to others (ensemble dilutes its performance)
  • Forecast horizon is so long that all models are guessing

For theoretical foundations, see Ray & Reich, 2018, PLOS Computational Biology.

Uncertainty Quantification: Beyond Point Forecasts

The cardinal sin of forecasting: Publishing a single number without uncertainty.

Example of bad forecasting: > “Model predicts 50,000 cases next week.”

If the real range is 30,000–80,000, decision-makers need to know. A point estimate without uncertainty bounds obscures exactly the information needed for planning.

Types of Uncertainty

1. Parameter Uncertainty

  • Do not know exact values of \(\beta\), \(\gamma\), etc.
  • Solution: Bayesian inference, confidence intervals

2. Structural Uncertainty

  • Model is wrong (all models are wrong)
  • SIR assumes homogeneous mixing (false)
  • Solution: Multiple models, ensemble

3. Stochastic Uncertainty

  • Random chance (superspreading events, who gets infected)
  • Solution: Monte Carlo simulation

4. Scenario Uncertainty

  • Future interventions unknown (for example, whether a lockdown will occur)
  • Solution: Scenario planning (best/worst/expected case)

Quantifying Forecast Uncertainty

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy.stats import norm
import seaborn as sns

# === METHOD 1: Bootstrap Confidence Intervals ===

def bootstrap_forecast(train_data, forecast_steps, n_bootstrap=1000):
 """
 Use bootstrap resampling to quantify uncertainty
 """
 forecasts = []

 for i in range(n_bootstrap):
  # Resample training data with replacement
  boot_sample = train_data.sample(n=len(train_data), replace=True)
  boot_sample = boot_sample.sort_values('day').reset_index(drop=True)

  # Fit model (simple exponential for speed)
  from scipy.optimize import curve_fit

  def exp_growth(x, a, b):
   return a * np.exp(b * x)

  try:
   params, _ = curve_fit(exp_growth, boot_sample['day'], boot_sample['cases'],
         p0=[100, 0.01], maxfev=5000)

   # Forecast
   forecast_days = np.arange(train_data['day'].iloc[-1] + 1,
          train_data['day'].iloc[-1] + 1 + forecast_steps)
   forecast = exp_growth(forecast_days, *params)
   forecasts.append(forecast)
  except:
   # If fit fails, use mean
   forecasts.append(np.full(forecast_steps, train_data['cases'].mean()))

 forecasts = np.array(forecasts)

 # Calculate percentiles
 forecast_mean = np.mean(forecasts, axis=0)
 forecast_lower = np.percentile(forecasts, 2.5, axis=0)
 forecast_upper = np.percentile(forecasts, 97.5, axis=0)

 return forecast_mean, forecast_lower, forecast_upper

# === METHOD 2: Bayesian Parameter Uncertainty ===

def bayesian_seir_forecast(train_data, forecast_steps, n_samples=500):
 """
 SEIR with Bayesian parameter uncertainty
 """
 N = 1000000

 # Illustrative prior distributions, chosen so R₀ = β/γ spans about 1.0–2.3
 beta_samples = np.random.uniform(0.12, 0.18, n_samples) # Transmission rate
 gamma_samples = np.random.uniform(0.08, 0.12, n_samples) # Recovery rate
 sigma_samples = np.random.uniform(0.15, 0.25, n_samples) # E→I rate

 forecasts = []

 for beta, gamma, sigma in zip(beta_samples, gamma_samples, sigma_samples):
  # Initial conditions
  I0 = train_data['cases'].iloc[0]
  E0 = I0
  R0_init = 0
  S0 = N - E0 - I0 - R0_init

  def seir(y, t, beta, sigma, gamma, N):
   S, E, I, R = y
   dS = -beta * S * I / N
   dE = beta * S * I / N - sigma * E
   dI = sigma * E - gamma * I
   dR = gamma * I
   return dS, dE, dI, dR

  # Solve
  t = np.arange(len(train_data) + forecast_steps)
  y0 = S0, E0, I0, R0_init

  solution = odeint(seir, y0, t, args=(beta, sigma, gamma, N))
  _, _, I, _ = solution.T

  forecasts.append(I[len(train_data):])

 forecasts = np.array(forecasts)

 forecast_mean = np.mean(forecasts, axis=0)
 forecast_lower = np.percentile(forecasts, 2.5, axis=0)
 forecast_upper = np.percentile(forecasts, 97.5, axis=0)

 return forecast_mean, forecast_lower, forecast_upper, forecasts

# === METHOD 3: Ensemble Spread as Uncertainty ===

def ensemble_uncertainty(individual_forecasts):
 """
 Use spread across ensemble members as uncertainty measure
 """
 forecasts_array = np.array(individual_forecasts)

 forecast_mean = np.mean(forecasts_array, axis=0)
 forecast_std = np.std(forecasts_array, axis=0)

 # 95% CI assuming normal distribution
 forecast_lower = forecast_mean - 1.96 * forecast_std
 forecast_upper = forecast_mean + 1.96 * forecast_std

 return forecast_mean, forecast_lower, forecast_upper

# Generate example data
np.random.seed(42)
days = np.arange(150)
cases = 100 * np.exp(0.025 * days) + np.random.normal(0, 50, len(days))
cases = np.maximum(cases, 10)

df = pd.DataFrame({'day': days, 'cases': cases})

train_size = 120
train = df.iloc[:train_size]
test = df.iloc[train_size:]

# Generate forecasts with uncertainty
print("Generating forecasts with uncertainty quantification...")

# Bootstrap
print(" Running bootstrap...")
boot_mean, boot_lower, boot_upper = bootstrap_forecast(train, len(test), n_bootstrap=500)

# Bayesian
print(" Running Bayesian SEIR...")
bayes_mean, bayes_lower, bayes_upper, bayes_samples = bayesian_seir_forecast(
 train, len(test), n_samples=200
)

# Visualize
fig, axes = plt.subplots(3, 1, figsize=(16, 16))

# Top: Bootstrap uncertainty
axes[0].plot(train['day'], train['cases'], 'o-',
   label='Training Data', color='blue', alpha=0.5, markersize=3)
axes[0].plot(test['day'], test['cases'], 'o-',
   label='Actual (Test)', color='black', linewidth=2, markersize=6)

axes[0].plot(test['day'], boot_mean, 's-',
   label='Bootstrap Forecast', color='red', linewidth=2, markersize=6)
axes[0].fill_between(test['day'], boot_lower, boot_upper,
      alpha=0.3, color='red', label='95% Bootstrap CI')

axes[0].axvline(train['day'].iloc[-1], color='green', linestyle='--', linewidth=2)

axes[0].set_ylabel('Daily Cases', fontsize=12)
axes[0].set_title('Bootstrap Uncertainty Quantification', fontsize=14)
axes[0].legend(fontsize=10)
axes[0].grid(True, alpha=0.3)

# Middle: Bayesian uncertainty
axes[1].plot(train['day'], train['cases'], 'o-',
   label='Training Data', color='blue', alpha=0.5, markersize=3)
axes[1].plot(test['day'], test['cases'], 'o-',
   label='Actual (Test)', color='black', linewidth=2, markersize=6)

axes[1].plot(test['day'], bayes_mean, '^-',
   label='Bayesian SEIR Forecast', color='purple', linewidth=2, markersize=6)
axes[1].fill_between(test['day'], bayes_lower, bayes_upper,
      alpha=0.3, color='purple', label='95% Credible Interval')

axes[1].axvline(train['day'].iloc[-1], color='green', linestyle='--', linewidth=2)

axes[1].set_ylabel('Daily Cases', fontsize=12)
axes[1].set_title('Bayesian Parameter Uncertainty', fontsize=14)
axes[1].legend(fontsize=10)
axes[1].grid(True, alpha=0.3)

# Bottom: Distribution of forecasts (spaghetti plot)
for i in range(min(50, len(bayes_samples))): # Plot subset for clarity
 axes[2].plot(test['day'], bayes_samples[i], '-',
    color='purple', alpha=0.1, linewidth=0.5)

axes[2].plot(train['day'], train['cases'], 'o-',
   label='Training Data', color='blue', alpha=0.7, markersize=3)
axes[2].plot(test['day'], test['cases'], 'o-',
   label='Actual (Test)', color='black', linewidth=3, markersize=6, zorder=5)
axes[2].plot(test['day'], bayes_mean, '-',
   label='Mean Forecast', color='red', linewidth=3, zorder=4)

axes[2].axvline(train['day'].iloc[-1], color='green', linestyle='--', linewidth=2)

axes[2].set_xlabel('Day', fontsize=12)
axes[2].set_ylabel('Daily Cases', fontsize=12)
axes[2].set_title('Forecast Uncertainty: Distribution of Possible Futures (Spaghetti Plot)', fontsize=14)
axes[2].legend(fontsize=10)
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('uncertainty_quantification.png', dpi=300)
plt.show()

# Evaluate coverage
def calculate_coverage(actual, lower, upper):
 """Calculate what % of actual values fall within prediction interval"""
 within_interval = (actual >= lower) & (actual <= upper)
 coverage = within_interval.mean() * 100
 return coverage

boot_coverage = calculate_coverage(test['cases'].values, boot_lower, boot_upper)
bayes_coverage = calculate_coverage(test['cases'].values, bayes_lower, bayes_upper)

print("\n" + "="*60)
print("UNCERTAINTY QUANTIFICATION EVALUATION")
print("="*60)
print(f"\nBootstrap 95% CI:")
print(f" Coverage: {boot_coverage:.0f}% (should be ~95%)")

print(f"\nBayesian 95% Credible Interval:")
print(f" Coverage: {bayes_coverage:.0f}% (should be ~95%)")

# Interval width (narrower is better, if coverage is maintained)
boot_width = (boot_upper - boot_lower).mean()
bayes_width = (bayes_upper - bayes_lower).mean()

print(f"\nMean Interval Width:")
print(f" Bootstrap: {boot_width:.1f} cases")
print(f" Bayesian: {bayes_width:.1f} cases")
The Coverage vs. Sharpness Trade-Off

Coverage: % of actual values that fall within the prediction interval

  • Good: 95% coverage for a 95% interval
  • Bad: Only 70% coverage (interval too narrow)

Sharpness: How narrow the prediction interval is

  • Good: Narrow intervals (more informative)
  • Bad: Wide intervals (“cases will be between 10 and 10,000”)

The trade-off:

  • Wide intervals → High coverage, low utility
  • Narrow intervals → Low coverage, misleading

Goal: Maximize sharpness while maintaining nominal coverage.

For rigorous evaluation, see Gneiting & Raftery, 2007, Journal of the American Statistical Association.

Case Study: COVID-19 Forecasting Hub

The COVID-19 Forecast Hub collected real-time forecasts of COVID-19 hospitalizations, cases, and deaths for the United States from 2020 to 2024.

It launched in April 2020 (Cramer et al., 2022).

Setup

Coordinated by: Reich Lab at UMass Amherst

Participants: More than 90 academic, industry, and independent research groups

Frequency: Weekly forecasts

Targets:

  • Incident deaths (1–4 weeks ahead)
  • Incident hospitalizations
  • Incident cases (when data reliable)

Ensemble method: Equal-weight mean of submitted forecasts at each quantile through the week ending July 18, 2020, then the median from the week ending July 25, 2020 (Cramer et al., 2022)

Use: CDC and state health departments for decision-making

Results

Cramer et al., 2022, PNAS evaluated death forecasts from April 2020 through October 2021:

Key findings:

Across 1- to 4-week horizons:

  • The ensemble’s relative weighted interval score was 0.61 versus the baseline
  • Four-week forecasts had about 76% more weighted interval score than 1-week forecasts
  • Uncertainty widened and probabilistic accuracy degraded with horizon

Challenges encountered:

  • Novel variants (Alpha and Delta) produced distribution shifts that degraded many models simultaneously during the April 2020–October 2021 evaluation, which ended before Omicron
  • Changing testing made case forecasts unreliable
  • Policy changes (lockdowns, mask mandates) were unpredictable
  • Behavior changes as people reacted to case counts
  • Data quality issues including reporting delays and definition changes

Lessons Learned

What Worked

Short-term forecasts (1–2 weeks ahead):

  • Useful for operational planning (hospital beds, staffing)
  • Ensemble approach was robust
  • Simple models often competitive with complex models

Transparency and collaboration:

  • Open data, open code, open evaluation
  • Rapid iteration and model improvement
  • Community learned from collective failures

Ensemble method:

  • Strong average performance across evaluated forecasts
  • Greater robustness when component errors differ
  • No guarantee during shared failure modes or rapid change

For detailed analysis, see COVID-19 Forecast Hub research publications.

What Did Not Work

Forecasts at 5–20 weeks:

For the two teams that submitted 20-week forecasts for all 50 states, average weighted interval score at 20 weeks was 2.9–4 times the 1-week score, with the largest increase from underprediction penalties (Cramer et al., 2022).

Forecasting through major disruptions:

  • New variants degraded models calibrated to earlier transmission and immunity conditions
  • Policy changes (lockdowns) invalidated forecasts immediately
  • As of September 28, 2021, during Delta and before Omicron, persistent large inaccuracies led the Hub ensemble to report only one-week-ahead case forecasts and 14-day-ahead hospitalization forecasts (Forecast Hub weekly report, December 28, 2021)

Point forecasts without uncertainty:

  • Some early models provided single numbers
  • Dangerously misleading
  • Led to bad decisions

The fundamental problem: COVID-19 exhibited extreme non-stationarity: the data-generating process changed every few months (new variants, vaccines, behavior).

When the rules keep changing, forecasting beyond short horizons becomes unreliable.

Beyond the US Hub: International and LMIC Forecasting

The COVID-19 Forecast Hub was the largest but not the only multi-model forecasting effort. The European COVID-19 Forecast Hub (now part of RespiCast following its October 2024 transition to ECDC) collected forecasts from 48 models across 32 countries. Its ensemble was more accurate than the baseline and all individual models for one- and two-week-ahead death forecasts, while case forecasts and longer horizons were less reliable (Sherratt et al., 2023).

These efforts share a limitation: they depend on countries with robust surveillance infrastructure and consistent case reporting. Modeling across the 47 WHO African region countries estimated that only 1.4% of SARS-CoV-2 infections (1 in 71) were captured by reported data, with only one in three COVID-19 deaths reported (Cabore et al., 2022). When over 98% of outcomes are missing from training data, standard forecasting approaches fail. Adapted methods (partially observed Markov models, hierarchical reconciliation across geographic levels, integration of seroprevalence surveys) become necessary.

The pattern extends beyond COVID-19. A formal multi-team dengue forecasting challenge across Iquitos, Peru and San Juan, Puerto Rico found that multimodel ensembles outperformed individual models, but forecast skill was lowest during high-incidence seasons, precisely when forecasts matter most (Johansson et al., 2019). Counterintuitively, models incorporating biological mechanisms had lower average skill than simpler statistical approaches, reinforcing that mechanistic complexity does not guarantee predictive accuracy.

Two 2026 retrospective dengue comparisons make the same point: across Brazil’s 27 capital cities at 1- to 4-week horizons, CatBoost outperformed a more complex recurrent network (GRU) in most capitals, reaching R² 0.92 in Goiânia, yet still had negative R² in Florianópolis and Porto Alegre (Silva et al., 2026), and across 11 Bangladeshi districts no single model performed best at every 1- to 6-month horizon, with SARIMAX best for one-step outbreak detection (F1 0.854) and Prophet best for pooled magnitude at horizons 2 to 6 (Fuad, Milki & Aziz, 2026).

A systematic review of 67 AI-based early warning studies confirmed this gap: no pooled sensitivity benchmark exists across systems, and concerns about data quality, bias, and deployment readiness remain unresolved (Villanueva-Miranda et al., 2025).

A different kind of input is a fixed numerical summary of each place, learned once from search, mobility, built-environment, and weather data and added to the model a team already runs. A 5 October 2026 preprint by Google researchers with partners at WHO AFRO, Oxford, NYU, Mount Sinai, and the University of Washington tested Google’s Population Dynamics Foundation Model (PDFM) this way, without retraining it for each task (Muslim et al., 2026; Google Research summary). In the Democratic Republic of the Congo, models were trained on surveillance data before November 2023 and tested over the next 89 weeks across 403 health zones. Adding the embeddings to case history did not help at 1 to 2 weeks, when recent cases carry most of the signal, but at 8 weeks it raised the share of each week’s five highest-ranked zones that went on to have a cholera emergence event from 0.36 to 0.42 (about 2.1 rather than 1.8 correct picks per week). Low internet use in the DRC meant search data had to be pooled over six months, which limits how quickly the embeddings can reflect an acute shock. In Mexico, adding the embeddings to the TimesFM 2.0 forecasting model improved one-month dengue forecasts overall, yet only 47.6% of 2,457 municipalities improved and the median municipality did slightly worse; gains concentrated where dengue was already circulating (64.6% to 72.1% of municipalities with active transmission, against 38.9% of those with no cases). These are unreviewed, developer-led evaluations on held-out surveillance data, not evidence that the shortlists changed vaccine, supply, or vector-control decisions or reduced cases.

See also Global Health and Equity.

When NOT to Forecast: The Ethics of Publishing Predictions

Forecasts have consequences. Sometimes, not forecasting is the ethical choice.

The Self-Fulfilling/Self-Defeating Problem Revisited

Scenario planning disguised as prediction:

Many “forecasts” during COVID-19 were actually scenarios:

  • “If we do nothing, 2 million will die”
  • “If we lock down, 100,000 will die”

These are conditional statements about policy choices, not predictions.

The problem: Media reports often presented them as predictions, which confused the public.

When to Say “We Don’t Know”

Do not publish forecasts when:

Novel pathogen with no historical data

  • Early COVID-19 (January 2020): Every model was pure speculation
  • Better: “Here’s a range of scenarios from mild to catastrophic”

Major intervention pending

  • Government about to announce lockdown
  • Forecasting what happens “if current trends continue” is misleading

Data quality too poor

  • Testing availability changing daily
  • Reporting systems overwhelmed
  • Better: Focus on qualitative assessment

Model disagreement too high

  • If ensemble predictions range from 10,000 to 500,000 cases, the forecast offers little decision value
  • Better: “Current uncertainty is too high for actionable forecasts”

Horizon too long

  • Forecasting 3+ months ahead for rapidly changing epidemic
  • Better: Provide scenarios, not forecasts
The Responsibility of Forecasters

Ioannidis, Cripps, and Tanner argued that many early COVID-19 forecasts failed and recommended “careful modeling of predictive distributions rather than focusing on point estimates” and “continuously reappraising models based on their validated performance” (Ioannidis et al., 2020).

Best practices:

  • Always provide uncertainty quantification
  • Clearly distinguish forecasts from scenarios
  • Update forecasts as data changes
  • Acknowledge model limitations publicly
  • Archive forecasts for retrospective evaluation
  • Resist pressure to forecast when uncertainty is unquantifiable

Practical Guidance: Building Your First Forecasting System

Pre-Flight Checklist

Before you start forecasting:

☐ Do you have sufficient historical data?

  • Time series: ≥2 years for seasonal diseases
  • SEIR models: At least early epidemic data to estimate parameters

☐ Is your data quality acceptable?

☐ What decisions will this forecast inform?

  • Hospital resource allocation → 1–2 week forecasts sufficient
  • Long-term policy planning → Need scenarios, not forecasts
  • Public communication → Uncertainty is essential

☐ What is the forecast horizon?

  • 1–2 weeks: Achievable with most methods
  • 3–4 weeks: Difficult, high uncertainty

  • 4 weeks: Scenario planning, not forecasting

Model Selection Strategy

Start simple, add complexity only when justified:

1. Naive baseline (persistence: tomorrow = today)
 ↓
2. Time series model (ARIMA, Prophet)
 ↓ (if you need scenario planning)
3. Mechanistic model (SEIR)
 ↓ (if you have large training data)
4. Machine learning (LSTM, ensemble)

Never skip the naive baseline. A complex model that does not beat “tomorrow will be like today” adds no value.

Evaluation Framework

Proper time series cross-validation:

# WRONG: Random train/test split
X_train, X_test = train_test_split(X, test_size=0.2) # Don't do this

# CORRECT: Respect temporal order
train_size = int(0.8 * len(X))
X_train = X[:train_size]
X_test = X[train_size:] # Correct approach

Report multiple metrics:

  • MAE (Mean Absolute Error)
  • RMSE (Root Mean Squared Error)
  • Coverage (95% CI should contain 95% of actuals)
  • Sharpness (average interval width)

Track performance over time:

  • Forecast accuracy degrades → Retrain models
  • Sudden drop in accuracy → Something changed (investigate)

Implementation Template

class EpidemicForecaster:
 """
 Template for an epidemic forecasting system
 """
 def __init__(self):
  self.models = {}
  self.performance_history = []

 def add_model(self, name, model):
  """Add a model to the ensemble"""
  self.models[name] = model

 def fit(self, train_data):
  """Fit all models"""
  for name, model in self.models.items():
   model.fit(train_data)

 def forecast(self, horizon, method='weighted_ensemble'):
  """Generate ensemble forecast"""
  # Get forecasts from all models
  forecasts = {}
  for name, model in self.models.items():
   forecasts[name] = model.predict(horizon)

  # Combine
  if method == 'simple_average':
   ensemble = np.mean(list(forecasts.values()), axis=0)
  elif method == 'weighted_ensemble':
   # Weight by recent performance
   weights = self._calculate_weights()
   ensemble = sum(w * f for w, f in zip(weights, forecasts.values()))

  return ensemble, forecasts

 def evaluate(self, actual, forecast):
  """Evaluate forecast and log performance"""
  mae = mean_absolute_error(actual, forecast)

  self.performance_history.append({
   'date': pd.Timestamp.now(),
   'mae': mae
  })

  return mae

 def _calculate_weights(self):
  """Calculate model weights based on recent performance"""
  # Implementation depends on your needs
  pass

# Usage
forecaster = EpidemicForecaster()
forecaster.add_model('seir', SEIRModel())
forecaster.add_model('arima', ARIMAModel())
forecaster.add_model('lstm', LSTMModel())

forecaster.fit(training_data)
forecast, individual_forecasts = forecaster.forecast(horizon=14)

Key Takeaways

  1. Forecasting ≠ Surveillance. Detection shows what is happening now; forecasting estimates what will happen. The latter is substantially harder.

  2. Error increases with horizon. Short-horizon forecasts can support operations, while longer horizons require wider uncertainty and comparison with appropriate baselines.

  3. Mechanistic models provide interpretable scenario planning. SEIR models can encode assumptions such as a vaccination scenario. Purely autoregressive ML models do not identify intervention effects without causal or mechanistic assumptions.

  4. Machine learning excels at pattern recognition in stable conditions. It breaks down when conditions change (new variants, interventions).

  5. Ensembles are a strong default, not a guarantee. Use diverse components, compare against simple baselines, and require prospective evidence at each operational horizon.

  6. Uncertainty quantification is mandatory. Never publish a single-number forecast. Always provide confidence/prediction intervals.

  7. Evaluation is essential. Track forecast performance over time. If accuracy degrades, investigate why and retrain.

  8. Sometimes, do not forecast. When uncertainty is unquantifiable or forecasts would be misleading, say “we don’t know.”

  9. Forecasts influence outcomes. The act of forecasting changes behavior, creating self-fulfilling or self-defeating prophecies.

  10. COVID-19 Forecast Hub demonstrated both ensemble value and shared failure modes. Novel variants can degrade many components at once.

  11. Most validated forecasting infrastructure is concentrated in high-income countries. FluSight, COVID-19 Forecast Hub, and the methods validated in this chapter rely on robust reporting systems that many countries lack. Forecasting in data-scarce settings (incomplete case reporting, limited sequencing capacity, comorbidity burdens that alter disease dynamics) requires adapted approaches, including hierarchical models that reconcile sparse subnational data with regional trends. See also Global Health and Equity.

Practice Exercises

Exercise 1: Implement and Compare Models

Build SIR, SEIR, ARIMA, and LSTM models for the same epidemic dataset. Which performs best at 1-week, 2-week, and 4-week horizons? Why?

Exercise 2: Sensitivity Analysis

For an SEIR model, vary \(\beta\) by ±20%. How much do peak timing and magnitude change? What does this tell you about parameter uncertainty?

Exercise 3: Ensemble Optimization

Given forecasts from 5 models, find the optimal weights that minimize MAE on a validation set. Do these weights generalize to the test set?

Exercise 4: Uncertainty Calibration

Generate forecasts with 95% prediction intervals. What percentage of actual values fall within your intervals? If not 95%, how can you calibrate?

Exercise 5: Scenario Planning

Using an SEIR model, create three scenarios: (1) No intervention, (2) 50% transmission reduction, (3) 80% transmission reduction. How different are the outcomes?

Check Your Understanding

Test your knowledge of epidemic forecasting methods and their appropriate applications.

Question 1

In a basic SIR compartmental model, you calculate R₀ = 2.5. What does this mean for epidemic dynamics, and what happens when R₀ drops to 0.9 due to interventions?

  1. R₀ = 2.5 means 2.5% of the population gets infected; R₀ = 0.9 means 0.9% gets infected
  2. R₀ = 2.5 means the epidemic grows exponentially; R₀ = 0.9 means the epidemic dies out
  3. R₀ values only matter for novel pathogens and do not change with interventions
  4. R₀ = 2.5 means each infected person has 2.5 contacts per day

Correct Answer: b) R₀ = 2.5 means the epidemic grows exponentially; R₀ = 0.9 means the epidemic dies out

The basic reproduction number (R₀) is the most important parameter in epidemic modeling:

What R₀ means:

  • R₀ = number of secondary infections caused by one infected person in a fully susceptible population
  • R₀ = 2.5: Each infected person infects 2.5 others → epidemic grows exponentially
  • R₀ = 0.9: Each infected person infects <1 person → epidemic dies out
  • R₀ = 1: Threshold (each infected person replaces themselves on average)

Critical threshold: R₀ > 1 → epidemic grows; R₀ < 1 → epidemic dies out

How interventions affect R₀:

  • R₀ = β/γ where β = transmission rate, γ = recovery rate
  • Social distancing reduces β → lowers R₀
  • Faster isolation reduces infectious period → increases γ → lowers R₀
  • A change from 2.5 to 2.3 is an 8% reduction and remains above 1; crossing below 1 changes expected growth to decline

Examples from the chapter’s SIR scenario code (N = 1,000,000, γ = 0.1):

  • R₀ = 2.5 (no intervention) → Peak: about 234,000 people infected at once
  • R₀ = 2.0 (distancing) → Peak: about 154,000 people infected at once
  • R₀ = 0.9 (lockdown) → No peak, epidemic dies out

Why other answers are wrong:

  • R₀ ≠ percentage infected (that is the attack rate)
  • R₀ ≠ contacts per day (it is contact rate × transmission probability × infectious period)
  • R₀ in this model changes with interventions, which is the purpose of intervening
Question 2

The COVID-19 Forecast Hub found that probabilistic error increased with horizon, with 4-week forecasts having about 76% more weighted interval score than 1-week forecasts. What fundamental challenge explains this horizon effect?

  1. Longer forecasts require more computational power that was not available
  2. Epidemics exhibit chaotic dynamics where small changes in initial conditions, behaviors, and interventions compound over time
  3. Four weeks is beyond the maximum incubation period so the data becomes unreliable
  4. Forecasters deliberately make long-term forecasts less accurate to avoid liability

Correct Answer: b) Epidemics exhibit chaotic dynamics where small changes in initial conditions, behaviors, and interventions compound over time

This is the butterfly effect in epidemiology, a fundamental characteristic of complex adaptive systems:

Why accuracy degrades with forecast horizon:

  1. Compounding uncertainty: Small errors in parameter estimation (β, γ) grow exponentially
  2. Unpredictable interventions: Lockdowns, mask mandates, vaccine campaigns cannot be predicted
  3. Behavior changes: People react to case counts, changing transmission dynamics
  4. Viral evolution: New variants with different R₀ values break model assumptions
  5. Reporting changes: Testing policy shifts invalidate surveillance data patterns
  6. Random events: Superspreading events can dramatically shift trajectories

From the chapter:

Forecast Horizon Typical Accuracy Usefulness
1 week ahead ±10-20% High - operational planning
2 weeks ahead ±20-40% Moderate - resource allocation
4 weeks ahead ±40-100% Low - scenario planning only
3+ months ahead Worse than baseline Essentially guesswork

The self-fulfilling/self-defeating prophecy:

  • Forecasts influence behavior, which changes outcomes
  • Unlike weather forecasting (rain does not care about predictions), epidemic forecasts affect the thing being forecasted
  • Long horizons amplify this feedback loop

Real-world example: As of September 28, 2021, during Delta and before Omicron, persistent inaccuracies led the U.S. Forecast Hub to stop publishing case forecasts beyond one week (Forecast Hub weekly report, December 28, 2021).

The lesson: Focus on short-term forecasts (1–2 weeks) for operational decisions. Treat longer horizons as scenario planning, not prediction.

Question 3

You are comparing mechanistic SEIR models versus LSTM deep learning models for epidemic forecasting. Which statement best describes their complementary strengths and appropriate use cases?

  1. LSTM models are always superior because they use more advanced AI technology
  2. SEIR models enable interpretable scenario planning (what if we vaccinate 50%?), while LSTMs excel at short-term pattern recognition in stable conditions
  3. SEIR models are outdated; modern forecasting should only use deep learning
  4. LSTM models understand causal mechanisms better than SEIR models

Correct Answer: b) SEIR models enable interpretable scenario planning (what if we vaccinate 50%?), while LSTMs excel at short-term pattern recognition in stable conditions

This captures the fundamental tension between mechanistic and machine learning approaches:

SEIR/Mechanistic Models:

Strengths:

  • Interpretable: Every parameter (β, γ, σ) has biological meaning
  • Scenario planning: Can simulate “what if we reduce transmission by 30%?”
  • Theory-grounded: Built on a century of epidemiological knowledge
  • Extrapolation: Can predict novel scenarios without historical data
  • Policy-relevant: Stakeholders understand “if R₀ drops to X, we get Y cases”

Limitations:

  • Simplistic assumptions (homogeneous mixing, constant parameters)
  • Do not capture heterogeneity (superspreaders, age structure)
  • Require knowing parameters (β, γ) that are hard to estimate

LSTM/Machine Learning Models:

Strengths:

  • Data-driven: Learn patterns without mechanistic assumptions
  • Flexible: Capture complex nonlinear relationships
  • Short-term accuracy: Excel at 1–2 week forecasts when conditions are stable
  • Multiple inputs: Can integrate cases + mobility + weather + interventions

Limitations:

  • No mechanistic understanding: Learn patterns, not causes
  • Break on distribution shift: New variants, interventions invalidate training data
  • Opaque: Cannot explain why predictions changed
  • Cannot answer counterfactuals: “What if we vaccinate?” → no historical data for that

The LSTM limitation example from the chapter:

  • Training data: epidemic curve + lockdown → cases drop
  • LSTM learns: “After cases reach X, they drop”
  • New epidemic: cases reach X, but no lockdown
  • LSTM incorrectly predicts a drop
  • Problem: No concept that lockdown caused the drop

Best practice (from COVID-19 Forecast Hub):

  • Use both in an ensemble
  • SEIR for scenario planning and understanding mechanisms
  • LSTM for short-term pattern recognition
  • Ensemble combines their strengths and averages out individual failures
Question 4

Why can an ensemble of forecasting models be more robust than an individual model?

  1. Ensemble methods use more computational resources so they are inherently more accurate
  2. Different models capture different patterns; averaging reduces overfitting and provides robustness when individual models fail
  3. Ensemble methods average out all the errors so they always perform better mathematically
  4. The ensemble was weighted toward the best model, which is why it performed well

Correct Answer: b) Different models capture different patterns; averaging reduces overfitting and provides robustness when individual models fail

Ensemble forecasting is one of the most important lessons from both the CDC FluSight Challenge and COVID-19 Forecast Hub, but the benefit must be measured rather than assumed:

Why ensembles work:

  1. Different models capture different patterns:
    • Mechanistic models (SEIR): Understand transmission dynamics
    • Time series models (ARIMA): Capture seasonal patterns
    • ML models (LSTM): Learn complex nonlinear relationships
    • Each model has blind spots; ensembles fill in gaps
  2. Averaging reduces overfitting:
    • Individual models may be overconfident in their predictions
    • Ensemble spreads risk across multiple approaches
    • Extreme predictions get averaged out
  3. Robustness to model-specific failures:
    • If one model fails (e.g., LSTM breaks on new variant), others compensate
    • No need to identify “best” model in advance
    • Diversity helps only when component errors are not identical

Evidence from forecast hubs:

  • The COVID-19 Forecast Hub ensemble had strong average performance across locations, targets, and horizons (Cramer et al., 2022, PNAS)
  • CDC’s 2024–2025 FluSight ensemble led its season-wide primary metric but missed rapid rises and peak changes
  • Equal weighting remains a defensible real-time default when evidence for learned weights is sparse or unstable

Ensemble methods:

  1. Simple average: Equal weight to all models, useful as a stable starting point
  2. Weighted average: Weight by historical performance
  3. Median ensemble: More robust to outliers
  4. Adaptive weights: Change based on recent performance

When ensemble does not help:

  • All models use same data and make same mistakes
  • One model is durably superior in prospective evaluation
  • Forecast horizon so long all models are guessing

From the chapter example: Models compared:

  • SEIR, ARIMA, ML → best individual MAE
  • Simple ensemble → compare with the best individual and a naive baseline
  • Weighted ensemble → retain only if validation shows a durable improvement

The lesson: An ensemble is an evaluated forecast, not a default claim of superiority. Retain it only when its calibration and proper scores justify operational use.

Question 5

You publish a forecast predicting “severe outbreak in 2 weeks” with wide confidence intervals (10,000–100,000 cases). The public panics and implements strict measures. The outbreak is mild (15,000 cases). What forecasting challenge does this illustrate?

  1. Your model had high bias and should have been regularized better
  2. The self-defeating prophecy problem: forecasts influence behavior, which changes outcomes, making the forecast appear wrong even though it drove the right action
  3. Confidence intervals were too wide and should have been narrower for credibility
  4. The model should have predicted the intervention and adjusted for it

Correct Answer: b) The self-defeating prophecy problem: forecasts influence behavior, which changes outcomes, making the forecast appear wrong even though it drove the right action

This is one of the most profound challenges in epidemic forecasting, unique to predictions that influence human behavior:

The paradox: Success makes forecasters look wrong

Scenario 1: Self-defeating prophecy (this question)

  • Model predicts catastrophic outbreak
  • Public panics, implements strict measures
  • Outbreak is mild
  • Conclusion: “Model was wrong! Overreaction!”
  • Reality: Model was right about unmitigated scenario; intervention changed trajectory

Scenario 2: Self-fulfilling prophecy

  • Model predicts mild outbreak
  • Public relaxes precautions
  • Outbreak explodes
  • Conclusion: “Model was wrong! Should have prepared!”
  • Reality: Model was right for baseline scenario; behavior change amplified outbreak

Unlike weather forecasting:

  • Rain forecast does not affect whether it rains
  • Epidemic forecasts DO affect transmission (people change behavior)
  • Creates an impossible situation for forecasters

The Michael Fish Problem (from the chapter):

  • 1987: BBC weatherman dismissed hurricane concerns
  • Hurricane hit anyway → 18 died, permanent trust damage
  • Lesson: Reassuring forecasts that turn out wrong destroy trust permanently
  • But alarming forecasts that do not materialize also erode trust (more slowly)

Ethical best practices:

Clearly distinguish forecasts from scenarios:

  • “If current trends continue…” (forecast)
  • “If we do nothing…” (scenario)
  • “If we implement lockdown…” (scenario)

Communicate uncertainty honestly:

  • Wide confidence intervals (10K–100K) show genuine uncertainty
  • Better than false precision (53,472 cases)

Frame conditionally:

  • “Without intervention: 50,000–100,000 cases”
  • “With strong measures: 10,000–30,000 cases”

Update as conditions change:

  • “Our forecast assumed no intervention; behavior changed, so we update”

From Ioannidis et al., 2020: Model predictive distributions rather than point estimates, and continuously reappraise models against validated performance.

When NOT to forecast:

  • Novel pathogen (no data)
  • Major intervention pending
  • Model disagreement too high
  • Better to say: “Current uncertainty is too high for actionable forecasts”
Question 6

You are building an epidemic forecasting system for influenza. Your ARIMA model achieves 12% MAPE, while a “naive” persistence model (tomorrow = today) achieves 15% MAPE. Your LSTM achieves 11% MAPE but requires 100x more computational resources. What should you deploy?

  1. LSTM model because it has the lowest error
  2. ARIMA model because it balances accuracy and simplicity, and beats the naive baseline
  3. Naive model because it is simplest and the differences are too small to matter
  4. An ensemble of all three models weighted by their inverse MAPE

Correct Answer: b) ARIMA model because it balances accuracy and simplicity, and beats the naive baseline

This question tests practical judgment about model selection and deployment tradeoffs:

Why ARIMA is the best choice here:

  1. Beats naive baseline (12% vs 15% MAPE)
    • Rule from the chapter: “A complex model that does not beat ‘tomorrow will be like today’ adds no value”
    • 20% improvement over baseline is meaningful
  2. Reasonable complexity
    • ARIMA is interpretable and well-understood
    • Can explain seasonal patterns, trends to stakeholders
    • Easier to diagnose when forecasts degrade
  3. Computational efficiency
    • Fast training and inference
    • Can update frequently as new data arrives
    • No specialized hardware needed
  4. Appropriate for influenza
    • Flu has strong seasonal patterns (ARIMA excels here)
    • Reporting is stable and consistent
    • Short-term forecasts (1–4 weeks) are the use case

Why NOT choose LSTM (11% MAPE):

  • Only 1 percentage point improvement over ARIMA (8% relative improvement)
  • 100x computational cost for marginal gain
  • Harder to maintain and debug
  • Opaque predictions harder to explain to health departments
  • Principle: Simplest model that beats baseline and meets requirements

Why NOT choose naive baseline:

  • 15% MAPE is worse than ARIMA
  • Does not capture seasonality
  • However, always compute it as a benchmark

Why ensemble might not be best:

  • More complex to maintain
  • Marginal improvement unlikely to justify complexity for weekly flu forecasts
  • Better for high-stakes scenarios (pandemic)

From the chapter’s model selection strategy:

1. Naive baseline (persistence: tomorrow = today)
 ↓ (beat this? Yes → continue)
2. Time series model (ARIMA, Prophet) ← this case
 ↓ (need scenario planning? No)
3. Mechanistic model (SEIR)
 ↓ (have large training data + need more accuracy? No)
4. Machine learning (LSTM, ensemble)

Best practice for production systems:

  1. Deploy ARIMA as primary model
  2. Monitor performance weekly
  3. Keep naive baseline for comparison
  4. If accuracy degrades, investigate and retrain
  5. Consider LSTM only if:
    • Accuracy requirements become critical
    • Computational resources available
    • Interpretability less important

Context matters:

  • Routine hospital bed planning (lower stakes): ARIMA is fine
  • Pandemic ventilator procurement (high stakes): Might justify LSTM or ensemble
  • Public communication: Need interpretable model (ARIMA)

The chapter emphasizes: “Never skip the naive baseline. A complex model that does not beat ‘tomorrow will be like today’ adds no value.”

Why is epidemic forecasting difficult?

Transmission depends on pathogen biology, immunity, behavior, policy, testing, reporting, and interventions that can change in response to the forecast itself. Small assumption errors compound with horizon, and a new variant or reporting regime can invalidate learned relationships. Forecasts should therefore be probabilistic, compared with a baseline, evaluated with proper scoring rules, and paired with explicit assumptions and update schedules.

What is the difference between SIR and SEIR models?

Both are compartmental models. SIR represents susceptible, infectious, and recovered states. SEIR adds an exposed state for people who are infected but not yet infectious. The additional compartment can better represent diseases with a meaningful latent period, but it also adds parameters and assumptions. Model choice should reflect the decision, data, disease process, and identifiability rather than a belief that more compartments automatically improve accuracy.

How should epidemic forecast accuracy be reported?

Report performance by target, horizon, location, and period, using scores that evaluate the full predictive distribution when uncertainty is provided. Point error alone can reward overconfident forecasts or hide poor interval coverage. The U.S. COVID-19 Forecast Hub evaluated forecasts with weighted interval score and relative measures rather than the invented universal percentage-error bands that once appeared in this chapter (Cramer et al., 2022).

When should ensembles be used?

Ensembles are useful when constituent models capture different plausible structures and are evaluated consistently. They can reduce dependence on one model’s failure mode and often perform well on average, but they cannot protect against a shared blind spot such as a new variant or broken input feed. The ensemble method, inclusion rules, weighting, missing forecasts, and version changes should be documented.

Can forecasts support scenario planning?

Yes, but a forecast and a scenario answer different questions. A forecast estimates what is likely under stated assumptions. A scenario examines what could happen under specified interventions or conditions. Mechanistic and causal assumptions are needed to interpret intervention scenarios. A purely autoregressive machine-learning model does not identify the effect of a policy merely because it predicts the next observation.

Discussion Questions

  1. SEIR models assume homogeneous mixing (everyone contacts everyone equally). How does this assumption break down in reality? How might you modify SEIR to be more realistic?

  2. The COVID-19 Forecast Hub ensemble outperformed individual models on average, and as of September 28, 2021, during Delta and before Omicron, persistent inaccuracies led the Hub to limit case forecasts to one week ahead. What does this tell us about the limits of forecasting? When should forecasts give way to scenario planning?

  3. Mechanistic models (SEIR) are interpretable but make strong assumptions. ML models (LSTMs) are flexible but opaque. For public health decision-making, which is more important: interpretability or accuracy? Does it depend on the decision?

  4. Imagine you are a state health department epidemiologist in March 2020. Multiple models predict vastly different COVID-19 death tolls (100K to 2M). How do you communicate this uncertainty to policymakers? To the public?

  5. Google Flu Trends failed because search behavior changed when people became aware of the model. Is this problem inherent to all epidemic forecasting, or specific to web-based surveillance?

  6. The COVID-19 Forecast Hub published forecasts weekly, even when model uncertainty was extremely high. Was this ethical? Should forecasters refuse to forecast when uncertainty is unquantifiable?

Further Resources

Books and Guides

Academic Papers

Forecasting Evaluation:

Machine Learning Methods:

Ethics and Limitations:

Tools and Platforms

  • COVID-19 Forecast Hub - Archive of all COVID-19 forecasts with open data and code
  • CDC FluSight - Annual influenza forecasting challenge (now targeting hospital admissions)
  • EpiModel R package - Network-based epidemic modeling framework
  • Delphi Epidata - CMU Delphi group’s Python API for epidemiological data and forecasting

Online Courses

Next: AI in Genomic Surveillance and Pathogen Analysis →