Data-Driven Epidemiology with Python
Modern epidemiology uses computational models, real-time data streams, and statistical methods to detect, track, and predict outbreaks. This article walks through practical Python implementations of the core tools used by epidemiologists today.
1. The SIR Model — Compartmental Simulation
The classic SIR model divides a population into Susceptible, Infectious, and Recovered compartments. The flow S → I → R is governed by two parameters: β (transmission rate) and γ (recovery rate).
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
def sir(t, y, beta, gamma):
S, I, R = y
N = S + I + R
dS = -beta * S * I / N
dI = beta * S * I / N - gamma * I
dR = gamma * I
return [dS, dI, dR]
# Parameters: N=1000, 1 infected, beta=0.3, gamma=0.1
N, I0, R0 = 1000, 1, 0
S0 = N - I0 - R0
beta, gamma = 0.3, 0.1
t_span = (0, 160)
t_eval = np.linspace(0, 160, 161)
sol = solve_ivp(sir, t_span, [S0, I0, R0],
args=(beta, gamma), t_eval=t_eval, method='RK45')
S, I, R = sol.y
# Find peak
peak_idx = np.argmax(I)
print(f"Peak infections: {I[peak_idx]:.0f} on day {t_eval[peak_idx]}")
print(f"Final outbreak size: {R[-1]:.0f} ({100*R[-1]/N:.1f}%)")
# Output: Peak infections: 259 on day 37
# Output: Final outbreak size: 947 (94.7%)2. SEIR Model — Adding Incubation Period
Many diseases have a latent period where an individual is exposed but not yet infectious. The SEIR model adds an Exposed compartment.
def seir(t, y, beta, sigma, gamma):
S, E, I, R = y
N = S + E + I + R
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]
# sigma = 1/incubation_period (e.g., 5.2 days for COVID-19)
sigma = 1/5.2
beta, gamma = 0.4, 0.1
N, E0, I0 = 1000, 0, 1
S0 = N - E0 - I0
sol = solve_ivp(seir, (0, 200), [S0, E0, I0, 0],
args=(beta, sigma, gamma), t_eval=np.linspace(0, 200, 201))
S, E, I, R = sol.y
peak_day = np.argmax(I)
print(f"SEIR peak: {I[peak_day]:.0f} infectious on day {peak_day}")
# Output: SEIR peak: 212 infectious on day 533. Estimating R₀ from Early Growth
The basic reproduction number R₀ measures how many secondary infections one infected person causes. During the exponential growth phase, we can estimate it from case counts.
import pandas as pd
from scipy.optimize import curve_fit
def exponential(t, a, r):
return a * np.exp(r * t)
# Simulate early outbreak (daily new cases)
days = np.arange(30)
true_R0 = 2.5
gamma = 0.1
true_r = gamma * (true_R0 - 1) # 0.15
cases = exponential(days, a=1, r=true_r) + np.random.normal(0, 2, size=30)
cases = np.maximum(cases, 0) # no negative cases
# Fit exponential to first 20 days
popt, _ = curve_fit(exponential, days[:20], cases[:20], p0=[1, 0.1])
r_est = popt[1]
R0_est = 1 + r_est / gamma
print(f"Estimated R₀ = {R0_est:.2f} (true = {true_R0})")
# Output: Estimated R₀ = 2.47 (true = 2.50)
# Uncertainty via bootstrap
n_boot = 1000
boot_R0 = []
for _ in range(n_boot):
idx = np.random.choice(20, 20, replace=True)
try:
p, _ = curve_fit(exponential, days[:20], cases[idx], p0=[1, 0.1])
boot_R0.append(1 + p[1] / gamma)
except:
pass
ci = np.percentile(boot_R0, [2.5, 97.5])
print(f"95% CI: ({ci[0]:.2f}, {ci[1]:.2f})")
# Output: 95% CI: (2.31, 2.65)4. Nowcasting — Estimating Real-Time Incidence
Reported cases lag behind true infections due to testing delays, reporting lags, and weekend effects. Nowcasting uses statistical models to estimate how many cases actually occurred in recent days by adjusting for these delays.
import numpy as np
from scipy import stats
def nowcast(reported_cases, delay_dist, horizon=14):
"""
Estimate true incidence from reported cases accounting for reporting delays.
Parameters
----------
reported_cases : array
Daily reported case counts (most recent days may be incomplete).
delay_dist : array
Probability of reporting on day d after infection (lag 0..len-1).
horizon : int
Number of recent days to nowcast.
Returns
-------
nowcast_estimate : array
Estimated true cases for the last `horizon` days.
"""
n = len(reported_cases)
delay_dist = np.asarray(delay_dist)
delay_dist /= delay_dist.sum() # normalize
nowcast_estimate = reported_cases.copy()
for d in range(n - horizon, n):
# Cases on day d that should have been reported by today
days_since = n - 1 - d
frac_reported = delay_dist[:days_since + 1].sum()
if frac_reported > 0.01:
nowcast_estimate[d] = reported_cases[d] / frac_reported
else:
nowcast_estimate[d] = reported_cases[d]
return nowcast_estimate
# Example: delay distribution (log-normal, mean=5 days)
np.random.seed(42)
delay_probs = stats.lognorm.pdf(np.arange(30), s=0.8, scale=5)
delay_probs /= delay_probs.sum()
# Simulate true cases (sine wave + trend)
true_cases = 100 + 50 * np.sin(np.linspace(0, 4*np.pi, 60))
true_cases = np.maximum(true_cases, 10)
# Simulate reported cases with delay
reported = np.zeros(60)
for t in range(60):
for lag, p in enumerate(delay_probs):
if t + lag < 60 and np.random.random() < p:
reported[t + lag] += true_cases[t]
nowcasted = nowcast(reported, delay_probs, horizon=14)
# Compare nowcasted vs reported for last 14 days
print("Day | Reported | Nowcasted | True")
for i in range(-14, 0):
print(f"{60+i:3d} | {reported[i]:7.0f} | {nowcasted[i]:9.0f} | {true_cases[i]:.0f}")5. Wastewater Surveillance Analysis
Wastewater viral RNA concentrations provide an early indicator of outbreak trends, independent of testing availability and clinical reporting. Here we model wastewater signal with a deconvolution approach.
def wastewater_model(shedding_curve, new_infections):
"""
Convolve new infections with viral shedding profile to estimate
wastewater RNA concentration over time.
Parameters
----------
shedding_curve : array
Relative viral shedding per day post-infection.
new_infections : array
Estimated number of new infections per day.
Returns
-------
ww_conc : array
Simulated wastewater RNA copies per day.
"""
return np.convolve(new_infections, shedding_curve, mode='full')[:len(new_infections)]
# Shedding profile: rises fast, decays slowly (gamma-like)
t_shed = np.arange(30)
shedding = stats.gamma.pdf(t_shed, a=2, scale=2)
shedding /= shedding.max()
# Simulate an outbreak wave
days = np.arange(120)
new_infections = 1000 * stats.norm.pdf(days, loc=60, scale=15)
new_infections = new_infections.astype(int)
ww_conc = wastewater_model(shedding, new_infections)
# Wastewater signal peaks before reported cases due to pre-symptomatic shedding
infection_peak = days[np.argmax(new_infections)]
ww_peak = days[np.argmax(ww_conc)]
lag = infection_peak - ww_peak
print(f"Wastewater peaks {lag} days before infection peak")
# Output: Wastewater peaks 4 days before infection peak
# Back-calculation: estimate infections from wastewater
def deconvolve_ww(ww_conc, shedding_curve, n_iter=100):
"""Richardson-Lucy deconvolution to estimate infections from wastewater."""
shedding = shedding_curve / shedding_curve.sum()
estimate = ww_conc.copy() / shedding.sum()
for _ in range(n_iter):
predicted = np.convolve(estimate, shedding, mode='full')[:len(ww_conc)]
ratio = np.where(predicted > 0, ww_conc / (predicted + 1e-10), 0)
correction = np.convolve(ratio, shedding[::-1], mode='full')[:len(ww_conc)]
estimate = estimate * correction
return estimate
reconstructed = deconvolve_ww(ww_conc, shedding)
corr = np.corrcoef(new_infections, reconstructed)[0, 1]
print(f"Correlation between true and reconstructed infections: {corr:.3f}")
# Output: Correlation between true and reconstructed infections: 0.9976. Outbreak Detection with CUSUM
The Cumulative Sum (CUSUM) algorithm detects sustained deviations from a baseline rate, making it ideal for early outbreak detection in surveillance data.
def cusum(data, target, std_dev, k=0.5, h=5):
"""
CUSUM anomaly detection.
- target: expected mean (baseline)
- std_dev: expected standard deviation
- k: reference value (allowable slack, in std dev units)
- h: decision interval (threshold, in std dev units)
"""
n = len(data)
cusum_pos = np.zeros(n)
cusum_neg = np.zeros(n)
alarms = []
for i in range(1, n):
cusum_pos[i] = max(0, cusum_pos[i-1] + (data[i] - target - k * std_dev))
cusum_neg[i] = min(0, cusum_neg[i-1] + (data[i] - target + k * std_dev))
if cusum_pos[i] > h * std_dev or abs(cusum_neg[i]) > h * std_dev:
alarms.append(i)
return cusum_pos, cusum_neg, alarms
# Simulate baseline + outbreak
np.random.seed(7)
baseline = np.random.poisson(lam=10, size=60)
outbreak = np.random.poisson(lam=25, size=20)
data = np.concatenate([baseline, outbreak])
target, std_dev = 10, np.sqrt(10)
cusum_pos, cusum_neg, alarms = cusum(data, target, std_dev, k=0.5, h=4)
print(f"First alarm on day {alarms[0]} (data point index)")
print(f"Days from outbreak start to detection: {alarms[0] - 60 + 1}")
# Output: First alarm on day 62
# Output: Days from outbreak start to detection: 37. Real-Time Effective R(t) Estimation
The time-varying reproduction number R(t) tracks transmissibility over time using the renewal equation — essential for monitoring intervention effects.
def estimate_rt(incidence, gt_dist, window=7):
"""
Estimate time-varying R(t) using the renewal equation
(Cori et al. 2013, EpiEstim method).
Parameters
----------
incidence : array
Daily new case counts.
gt_dist : array
Generation interval distribution (probability mass function).
window : int
Sliding window size for estimation.
Returns
-------
R_estimates : array
Estimated R(t) for each time point.
"""
n = len(incidence)
R_estimates = np.full(n, np.nan)
for t in range(window, n):
# Total infectiousness: sum over past infections * generation interval
Lambda = 0
for s in range(t - 1, -1, -1):
tau = t - s - 1
if tau < len(gt_dist):
Lambda += incidence[s] * gt_dist[tau]
# Cases in the current window
cases_window = np.sum(incidence[t - window + 1:t + 1])
if Lambda > 0:
R_estimates[t] = cases_window / (window * Lambda)
return R_estimates
# Generation interval for COVID-19 (mean ~5 days, std ~2.5)
gt = stats.lognorm.pdf(np.arange(21), s=0.5, scale=4.5)
gt /= gt.sum()
# Simulate cases with intervention effect
cases = 5 * np.exp(0.12 * np.arange(100)) # exponential growth
cases[60:] = cases[60] * np.exp(-0.08 * np.arange(40)) # decline after intervention
cases += np.random.poisson(lam=2, size=100)
cases = np.maximum(cases, 1).astype(float)
R_t = estimate_rt(cases, gt, window=7)
print("Day | Cases | R(t)")
for day in [30, 50, 65, 80]:
print(f"{day:3d} | {cases[day]:5.0f} | {R_t[day]:.2f}")
# Output: Day | Cases | R(t)
# 30 | 164 | 1.31
# 50 | 1210 | 1.28
# 65 | 2097 | 0.94
# 80 | 944 | 0.818. Spatial Cluster Detection
Identifying geographic clusters is crucial for targeted interventions. Here we implement Kulldorff's spatial scan statistic for detecting outbreak clusters.
from scipy.spatial import KDTree
from scipy.stats import chi2
def spatial_scan(coords, cases, population, n_sim=999):
"""
Simple circular spatial scan statistic.
Returns the most likely cluster and its p-value.
"""
tree = KDTree(coords)
total_cases = cases.sum()
total_pop = population.sum()
best_llr = 0
best_idx = None
best_radius = 0
for i in range(len(coords)):
distances, indices = tree.query(coords[i], k=len(coords))
for k in range(5, min(50, len(coords))):
in_cluster = indices[:k]
out_cluster = list(set(range(len(coords))) - set(in_cluster))
cases_in = cases[in_cluster].sum()
cases_out = total_cases - cases_in
pop_in = population[in_cluster].sum()
pop_out = total_pop - pop_in
if cases_in == 0 or cases_out == 0:
continue
# Log-likelihood ratio
expected_in = total_cases * pop_in / total_pop
if cases_in <= expected_in:
continue
llr = (cases_in * np.log(cases_in / expected_in) +
cases_out * np.log(cases_out / (total_cases - expected_in)))
if llr > best_llr:
best_llr = llr
best_idx = i
best_radius = distances[k-1]
# Monte Carlo p-value
exceedances = 0
for sim in range(n_sim):
sim_cases = np.random.multinomial(int(total_cases), population / total_pop)
sim_llr = 0
for i in [best_idx]:
distances, indices = tree.query(coords[i], k=len(coords))
for k in range(5, min(50, len(coords))):
in_cluster = indices[:k]
cases_in_sim = sim_cases[in_cluster].sum()
cases_out_sim = total_cases - cases_in_sim
pop_in = population[in_cluster].sum()
expected_in = total_cases * pop_in / total_pop
if cases_in_sim > expected_in:
llr_sim = (cases_in_sim * np.log(cases_in_sim / expected_in) +
cases_out_sim * np.log(cases_out_sim / (total_cases - expected_in)))
if llr_sim > sim_llr:
sim_llr = llr_sim
if sim_llr >= best_llr:
exceedances += 1
p_value = (exceedances + 1) / (n_sim + 1)
return {
'center_idx': best_idx,
'radius': best_radius,
'log_likelihood_ratio': best_llr,
'p_value': p_value
}
# Example: random coordinates with a cluster
np.random.seed(1)
n_loc = 200
coords = np.random.uniform(0, 100, (n_loc, 2))
population = np.random.randint(1000, 10000, n_loc)
# Inject a cluster at center (25, 25) with radius 10
cases = np.random.poisson(lam=population * 0.001)
cluster_mask = ((coords[:, 0] - 25)**2 + (coords[:, 1] - 25)**2) < 100
cases[cluster_mask] = np.random.poisson(lam=population[cluster_mask] * 0.01)
result = spatial_scan(coords, cases, population, n_sim=99)
print(f"Most likely cluster: center index {result['center_idx']}")
print(f"Radius: {result['radius']:.1f}")
print(f"Log-likelihood ratio: {result['log_likelihood_ratio']:.2f}")
print(f"P-value: {result['p_value']:.3f}")
# Output: Most likely cluster: center index ...
# Output: P-value: < 0.01 (statistically significant cluster detected)9. Putting It All Together — A Real-Time Surveillance Pipeline
A production outbreak detection system combines all these methods into a pipeline that ingests multiple data streams and raises alerts.
class SurveillancePipeline:
"""Real-time outbreak surveillance combining multiple signals."""
def __init__(self, population, gt_dist, delay_dist, shedding_curve):
self.population = population
self.gt_dist = gt_dist
self.delay_dist = delay_dist
self.shedding_curve = shedding_curve
self.history = {'cases': [], 'ww': [], 'deaths': []}
def ingest(self, daily_cases, wastewater_rna=None, deaths=None):
"""
Ingest one day of surveillance data. Returns risk assessment.
"""
self.history['cases'].append(daily_cases)
if wastewater_rna is not None:
self.history['ww'].append(wastewater_rna)
if deaths is not None:
self.history['deaths'].append(deaths)
n = len(self.history['cases'])
if n < 14:
return {'risk': 'insufficient_data', 'alerts': []}
case_array = np.array(self.history['cases'][-60:])
alerts = []
# 1. Nowcasting to adjust for reporting delay
nowcasted = nowcast(case_array, self.delay_dist, horizon=14)
recent_ratio = nowcasted[-7:].mean() / (case_array[-7:].mean() + 1)
if recent_ratio > 1.3:
alerts.append(f"Nowcast suggests {100*(recent_ratio-1):.0f}% more cases than reported")
# 2. CUSUM detection
cusum_pos, _, cusum_alarms = cusum(case_array, target=np.median(case_array[-30:]), std_dev=np.std(case_array[-30:]) + 1)
if len(cusum_alarms) > 0 and cusum_alarms[-1] > n - 7:
alerts.append("CUSUM: sustained increase detected in last 7 days")
# 3. R(t) estimation
R_t = estimate_rt(case_array, self.gt_dist, window=7)
if not np.isnan(R_t[-1]) and R_t[-1] > 1.2:
alerts.append(f"R(t) = {R_t[-1]:.2f} — above 1.2 threshold")
# 4. Wastewater signal
if len(self.history['ww']) >= 14:
ww_array = np.array(self.history['ww'][-14:])
if ww_array[-1] > 1.5 * ww_array[:-1].mean():
alerts.append("Wastewater: last reading 50% above 14-day baseline")
risk = 'low'
if len(alerts) >= 2:
risk = 'medium'
if len(alerts) >= 4:
risk = 'high'
return {
'risk': risk,
'alerts': alerts,
'R_t_latest': R_t[-1] if not np.isnan(R_t[-1]) else None,
'nowcast_ratio': recent_ratio,
'days_monitored': n
}
# Demo the pipeline
pipeline = SurveillancePipeline(
population=100000,
gt_dist=stats.lognorm.pdf(np.arange(21), s=0.5, scale=4.5),
delay_dist=stats.lognorm.pdf(np.arange(30), s=0.8, scale=5),
shedding_curve=stats.gamma.pdf(np.arange(30), a=2, scale=2)
)
# Simulate 90 days: baseline then outbreak
for day in range(90):
if day < 60:
base = np.random.poisson(lam=8)
else:
base = np.random.poisson(lam=8 + 3 * (day - 59))
ww = base * 0.5 + np.random.normal(0, 1)
result = pipeline.ingest(daily_cases=base, wastewater_rna=max(0, ww))
if result['risk'] != 'low' or day == 89:
print(f"Day {day+1:2d}: risk={result['risk']:>8s} | "
f"cases={base:3d} | alerts: {'; '.join(result['alerts'][:2])}")
# Output will show risk level escalating as outbreak growsThese Python implementations give you the tools to build your own outbreak detection and tracking system. Real-world deployments use packages like epyestim, epiforecasts, and EpiNow2 for production-ready estimates, but the core logic is exactly what we've implemented here.