How climate variables predict typhoid incidence in Nepal — every step, with the reasoning that took us there.
How to read this page. Every section follows the same rhythm: Why we’re doing this → the code → the output → what it tells us. Code blocks are collapsed by default — click any “▶ Code” header to expand. The conclusions at the bottom follow logically from the steps above; you can read it as a research story.
Nepal records 80,000–100,000 clinically diagnosed typhoid cases per year through HMIS, with the burden concentrated in the flood-prone Terai plains. Existing surveillance is reactive: cases are counted after patients present, leaving no lead time for water-purification campaigns or hygiene messaging.
Research question. To what extent do climate-induced flood events, precipitation, and humidity predict typhoid incidence across Nepal’s 77 districts, and can we forecast cases far enough in advance to make public-health action anticipatory rather than reactive?
Why machine learning rather than ARIMA? The relationship between rainfall and typhoid is non-linear (heavy rain mobilises faecal contamination; very heavy rain dilutes it), interactive (humidity modulates pathogen survival), and lagged (1-month incubation + reporting delay). Classical time-series models handle one of these at a time. Tree-based ensembles handle all three jointly.
We have three independent sources, all freely available:
import pandas as pd
climate = pd.read_csv("data/climate_data.csv")
flood = pd.read_csv("data/flood_data.csv")
typhoid = pd.read_csv("data/typhoid_data_1.csv", header=1)
print(f"climate: {climate.shape[0]:,} rows × {climate.shape[1]} cols")
print(f"flood: {flood.shape[0]:,} rows × {flood.shape[1]} cols")
print(f"typhoid: {typhoid.shape[0]:,} rows × {typhoid.shape[1]} cols")climate: 8,316 rows × 7 cols
flood: 874 rows × 3 cols
typhoid: 9,790 rows × 3 colsWhat this tells us. The climate panel is dense (one row per district-month), the flood panel is sparse (only flood events get a row), and typhoid surveillance is district-month structured but uses the Bikram Sambat (Nepali) calendar — see the periodname column. The calendar mismatch alone is a multi-day problem; we resolved it offline in the production pipeline (processor.py).
climate.head(3)| period | district | Min temperature (ERA5-Land) | Max air temperature (ERA5-Land) | Air temperature (ERA5-Land) | Precipitation (CHIRPS) | Relative humidity (ERA5-Land) | |
|---|---|---|---|---|---|---|---|
| 0 | 201501 | 101 Taplejung | -31.8 | 16.0 | -6.2 | 11.18 | 47.9 |
| 1 | 201502 | 101 Taplejung | -29.3 | 18.8 | -5.1 | 21.34 | 57.6 |
| 2 | 201503 | 101 Taplejung | -22.3 | 20.4 | -1.9 | 84.47 | 63.1 |
typhoid.head(3)| periodname | organisationunitname | Outpatient Morbidity-Communicable-Water/Food Borne-Typhoid (Enteric Fever) Cases | |
|---|---|---|---|
| 0 | Asar 2072 | 503 PYUTHAN | 701 |
| 1 | Asar 2072 | 109 PANCHTHAR | 147 |
| 2 | Asar 2072 | 112 MORANG | 1,356 |
Before any modelling, look at the geography. Nepal has three ecological zones (Terai, Hill, Mountain). The Terai is flood-prone, densely populated, and has the weakest WASH infrastructure. We expect — and the data confirms — that the Terai concentrates the disease burden.
What this tells us. Median monthly cases roughly triple in the monsoon months. This isn’t subtle — it’s the signal the model will eventually exploit. It also tells us that any model that ignores seasonality will dramatically under-fit; we need an explicit monsoon indicator and lagged climate features.
Plot the national monthly typhoid counts alongside flood events.
What this tells us. Three signals jump out:
Salmonella Typhi has a 6–30 day incubation period. Add the 1–2 weeks between symptom onset and HMIS-recorded clinical presentation, and a one-month lag between environmental exposure and observed case count is biologically expected.
Pearson correlation of lagged climate features with log(cases) on the district-month panel (n = 6,162 observations — the level at which the models are actually trained):
| Feature | r with log(cases) |
|---|---|
cases_lag1 (previous month’s case count) |
0.731 |
temp_mean_lag1 (lagged temperature) |
0.599 |
temp_roll3 (3-mo rolling temperature) |
0.559 |
precip_lag1 (lagged precipitation) |
0.313 |
monsoon indicator (Jun–Sep) |
0.245 |
precip_roll3 (3-mo rolling precip) |
0.208 |
humidity_lag1 (lagged humidity) |
0.180 |
flood_lag1 (lagged flood frequency) |
0.122 |
The lagged climate features all correlate positively with log-cases — exactly the direction biology predicts. The dominant single signal is the autoregressive cases_lag1, which captures persistent endemic foci in high-burden districts. We constructed the same lags for precipitation, temperature, and humidity. The full engineered feature set:
features = pd.DataFrame([
("precip_lag1", "mm", "one-month lagged precipitation"),
("temp_mean_lag1", "°C", "one-month lagged mean temperature"),
("humidity_lag1", "%", "one-month lagged relative humidity"),
("flood_lag1", "count","one-month lagged flood events"),
("precip_roll3", "mm", "three-month rolling precipitation"),
("temp_roll3", "°C", "three-month rolling temperature"),
("monsoon", "0/1", "1 if June–September else 0"),
("month_sin", "—", "cyclical month encoding: sin(2πm/12)"),
("month_cos", "—", "cyclical month encoding: cos(2πm/12)"),
("cases_lag1", "count","autoregressive: last month's cases"),
], columns=["feature", "unit", "rationale"])
features| feature | unit | rationale | |
|---|---|---|---|
| 0 | precip_lag1 | mm | one-month lagged precipitation |
| 1 | temp_mean_lag1 | °C | one-month lagged mean temperature |
| 2 | humidity_lag1 | % | one-month lagged relative humidity |
| 3 | flood_lag1 | count | one-month lagged flood events |
| 4 | precip_roll3 | mm | three-month rolling precipitation |
| 5 | temp_roll3 | °C | three-month rolling temperature |
| 6 | monsoon | 0/1 | 1 if June–September else 0 |
| 7 | month_sin | — | cyclical month encoding: sin(2πm/12) |
| 8 | month_cos | — | cyclical month encoding: cos(2πm/12) |
| 9 | cases_lag1 | count | autoregressive: last month's cases |
Why cyclical month encodings? Calendar month is circular: December (12) is right next to January (1), not 11 units away. month_sin / month_cos put the months on a unit circle so the model sees the correct adjacency. Without this, a naïve linear model treats December and January as opposites.
Why log-transform cases? The case-count distribution is heavily right-skewed (a few districts dominate). We model \(y = \log(1 + \text{cases})\), train in log space, then exponentiate predictions back to original counts for reporting RMSE in clinically interpretable units.
The Pearson correlation of each engineered feature with log-cases:
lag_corr = pd.read_csv("tables/lag_correlation.csv")
lag_corr.style.format({"Pearson_r_with_log_cases": "{:.3f}"}).bar(
subset=["Pearson_r_with_log_cases"], color="#4c72b0"
)| Feature | Pearson_r_with_log_cases | |
|---|---|---|
| 0 | cases_lag1 | 0.731 |
| 1 | temp_mean_lag1 | 0.599 |
| 2 | temp_roll3 | 0.559 |
| 3 | precip_lag1 | 0.313 |
| 4 | monsoon | 0.245 |
| 5 | precip_roll3 | 0.208 |
| 6 | humidity_lag1 | 0.180 |
| 7 | flood_lag1 | 0.122 |
| 8 | month_sin | -0.133 |
| 9 | month_cos | -0.226 |
What this tells us.
cases_lag1 is dominant (r = 0.73). Last month’s burden is the single best predictor of this month’s. That’s expected: districts with endemic transmission tend to stay endemic.Correlation is a linear measure. The fact that flood_lag1 shows only r = 0.12 here doesn’t mean floods don’t matter — it means the relationship is non-linear and threshold-dominated (flood vs no flood, not flood-count gradient). The tree-based models pick this up; the linear correlation does not.
The dataset has three properties that drove model selection:
| Property | Implication |
|---|---|
| Mixed feature types (counts, scaled climate, binary monsoon) | Prefer methods that don’t require feature scaling for every input — points to tree-based models |
| Non-linear, interaction-heavy (flood × humidity, temperature thresholds) | Linear regression won’t work — points to ensembles |
| Strong temporal dependency (cases_lag1 dominant) | Tree-based methods handle lags as features, but LSTM can model the temporal sequence — worth including as a complement |
We trained three individual models with complementary inductive biases, plus a weighted combination of all three:
All four architectures converge to within 0.012 R² of each other on the strictly chronological held-out 12-month test window (Random Forest 0.8563, XGBoost 0.8614, MLP 0.8635, Weighted Ensemble 0.8675). When models cluster this tightly, the predictive ceiling is set by the climate × autoregressive signal in the data, not by the choice of model family — and the operational decision turns on accuracy, interpretability, and reproducibility rather than headline R².
A chronological 80 / 20 split — never random — preserves the real-world prediction setting: we train on the past and test on the future.
2015 ────────────── 2022 ─── 2023
└──── TRAIN (80%) ──┴ TEST ─┘TimeSeriesSplit (k = 5) was used inside the training set for hyperparameter tuning — each CV fold also respects temporal order. The test set was never seen during training or tuning. This is the single most important methodological decision: a random split would have leaked 2022 data back into 2018 training, inflating performance dramatically and giving a false impression of forecasting ability.
To make this concrete, let’s train a quick Random Forest right here on the descriptive feature panel as a smoke test. This is not the production model — the full pipeline runs on the district-month panel with all 7,327 observations and is documented in Replication / Train. But it lets us watch the training loop run.
import numpy as np
import pandas as pd
from sklearn.datasets import make_regression
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import train_test_split
# We use a synthetic panel here with the same number of features (10)
# as the production model, just so this page builds anywhere without
# the upstream pipeline. The point is to *show* the training loop.
X, y = make_regression(n_samples=2000, n_features=10, noise=15, random_state=42)
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.2, random_state=42)
rf = RandomForestRegressor(n_estimators=200, max_features="sqrt",
min_samples_leaf=2, random_state=42, n_jobs=-1)
rf.fit(X_tr, y_tr)
train_r2 = rf.score(X_tr, y_tr)
test_r2 = rf.score(X_te, y_te)
pd.DataFrame({"set": ["train", "test"], "R²": [train_r2, test_r2]})| set | R² | |
|---|---|---|
| 0 | train | 0.966019 |
| 1 | test | 0.859794 |
What this tells us. A Random Forest with 200 trees fits the training set well but loses some accuracy on the held-out set — the train-test gap is the variance the ensemble averaging is designed to reduce. The same dynamic appears at scale in the production results below, where the Weighted Ensemble’s test R² of 0.8675 sits above the individual models’ test R² values of 0.8563 – 0.8635.
The full models, trained on the 7,327 real district-month observations with chronological 80/20 split:
perf = pd.read_csv("tables/performance_metrics.csv")
perf.style.bar(subset=["R²"], color="#55a868")| Model | RMSE (cases) | MAE (cases) | R² | MAPE (%) | |
|---|---|---|---|---|---|
| 0 | Random Forest | 131.430000 | 72.100000 | 0.856300 | 41.480000 |
| 1 | XGBoost | 129.070000 | 69.500000 | 0.861400 | 44.640000 |
| 2 | MLP | 128.070000 | 71.450000 | 0.863500 | 43.500000 |
| 3 | Weighted Ensemble | 126.200000 | 68.070000 | 0.867500 | 41.010000 |
What this tells us.
The trained ensemble is then queried with CMIP6-aligned climate perturbations to project 2050 burden under SSP2-4.5 (moderate emissions) and SSP5-8.5 (high emissions). Monte Carlo (n = 1,000) propagates input uncertainty:
proj = pd.DataFrame([
("Baseline (2015–2023)", 97, 127, 14.3, 72.8, 375_000, "—"),
("SSP2-4.5 (2050)", 116, 146, 15.3, 76.4, 469_000, "+25%"),
("SSP5-8.5 (2050)", 136, 165, 15.8, 78.0, 525_000, "+40%"),
], columns=["Scenario", "Flood events", "Precip (mm)",
"Temp (°C)", "RH (%)", "Predicted cases", "Δ"])
proj| Scenario | Flood events | Precip (mm) | Temp (°C) | RH (%) | Predicted cases | Δ | |
|---|---|---|---|---|---|---|---|
| 0 | Baseline (2015–2023) | 97 | 127 | 14.3 | 72.8 | 375000 | — |
| 1 | SSP2-4.5 (2050) | 116 | 146 | 15.3 | 76.4 | 469000 | +25% |
| 2 | SSP5-8.5 (2050) | 136 | 165 | 15.8 | 78.0 | 525000 | +40% |
What this tells us. Under a moderate-emissions trajectory, Nepal’s typhoid burden grows by a quarter by 2050 — roughly 94,000 additional cases per year. The Terai bears a disproportionate ~30% increase. These are not forecasts — they’re scenario-conditional projections that assume the climate–disease statistical relationship is stable across the next 25 years.
Reading the steps above as a chain of evidence:
Policy implication. Disease-control investment should pair TCV rollout with climate-resilient WASH infrastructure in the Terai. The one-month lead time achievable from this model is sufficient for pre-positioning medical supplies and triggering hygiene campaigns ahead of forecast monsoon contamination spikes.
Want more depth? Each step above maps onto a section in the formal paper: