Python notebook — runs manuscript/main.py with all supporting modules

Manuscript · Notebook

This notebook embeds the entire production pipeline for the research manuscript inline so that the page is reproducible without any external .py files. The pipeline corresponds to:

cd manuscript
python3 main.py

Each cell below is shown as runnable source. Copy any block into your environment, or save them as loader.py, processor.py, etc. under manuscript/src/, alongside main.py, and run the command above to regenerate every figure, table, and analysis artifact in manuscript/output/.

0 — The question

Nepal records 80,000–100,000 clinically diagnosed typhoid cases per year through HMIS, with burden concentrated in the flood-prone Terai. Existing surveillance is reactive — cases counted after patients present, leaving no lead time for water-purification campaigns. We ask whether climate-driven flood events, precipitation, and humidity predict typhoid incidence across Nepal’s 77 districts well enough to make public-health action anticipatory.

Why ML over ARIMA? The rainfall–typhoid relationship is non-linear (heavy rain mobilizes faecal contamination; very heavy rain dilutes it), interactive (humidity modulates pathogen survival), and lagged (one-month incubation + reporting delay). Classical time-series models handle one of these at a time; tree-based ensembles handle all three jointly.

1 — Load the raw data

Three independent sources, all freely available, harmonized to district-month resolution.

Code 
import pandas as pd

climate = pd.read_csv("manuscript/data/climate_data.csv")
flood   = pd.read_csv("manuscript/data/flood_data.csv")
typhoid = pd.read_csv("manuscript/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")

The climate panel is dense (one row per district-month); the flood panel is sparse (only flood events get a row); the 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, resolved in the production pipeline (processor.py, embedded below).

2 — Feature engineering

Code 
features = [
    ("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"),
]

Why cyclical encodings? Calendar month is circular — December (12) is next to January (1), not 11 units away. Why log-transform cases? The case-count distribution is heavily right-skewed; we model log(1 + cases) and back-transform predictions.

A.1 loader.py — robust ingestion

Defensive CSV reading for HMIS, CHIRPS, ERA5-Land, and DRR exports.

Code 
"""loader.py — robust ingestion of all raw data files."""

import os
import pandas as pd


class DataLoader:
    DATA_DIR = "manuscript/data"
    TYPHOID_FILES = ["typhoid_data_1.csv", "typhoid_data_2.csv"]
    CLIMATE_FILE  = "climate_data.csv"
    FLOOD_FILE    = "flood_data.csv"

    def load_typhoid_data(self) -> pd.DataFrame:
        frames = []
        for fname in self.TYPHOID_FILES:
            path = os.path.join(self.DATA_DIR, fname)
            if not os.path.exists(path):
                continue
            for skip in range(3):
                try:
                    df = pd.read_csv(path, skiprows=skip, low_memory=False)
                    df.columns = [c.strip() for c in df.columns]
                    if "periodname" in df.columns:
                        frames.append(df); break
                except Exception:
                    continue
        return pd.concat(frames, ignore_index=True) if frames else pd.DataFrame()

    def load_climate_data(self):
        path = os.path.join(self.DATA_DIR, self.CLIMATE_FILE)
        if not os.path.exists(path): return pd.DataFrame()
        df = pd.read_csv(path); df.columns = [c.strip() for c in df.columns]
        return df

    def load_flood_data(self):
        path = os.path.join(self.DATA_DIR, self.FLOOD_FILE)
        if not os.path.exists(path): return pd.DataFrame()
        df = pd.read_csv(path); df.columns = [c.strip() for c in df.columns]
        return df

    def load_all(self):
        return {
            "typhoid": self.load_typhoid_data(),
            "climate": self.load_climate_data(),
            "flood":   self.load_flood_data(),
        }

A.2 processor.py — calendar + district harmonization

Code 
"""processor.py — normalize district names, BS→AD calendar, IQR clipping, merge."""

import re
import numpy as np
import pandas as pd

COL_MAP = {
    "Min temperature (ERA5-Land)":      "temp_min",
    "Max air temperature (ERA5-Land)":  "temp_max",
    "Air temperature (ERA5-Land)":      "temp_mean",
    "Precipitation (CHIRPS)":           "precip",
    "Relative humidity (ERA5-Land)":    "humidity",
}


class DataProcessor:
    BS_MONTH_TO_AD = {
        "Baisakh": 4, "Jestha": 5,  "Asar":    6, "Shrawan": 7,
        "Bhadra":  8, "Ashwin": 9,  "Kartik": 10, "Mangsir": 11,
        "Poush":  12, "Magh":   1,  "Falgun":  2, "Chaitra":  3,
    }

    @staticmethod
    def _normalize_district(name):
        if pd.isna(name): return "Unknown"
        return re.sub(r"^\d+\s+", "", str(name)).strip().upper()

    def _bs_to_ad(self, period_str):
        try:
            parts = str(period_str).split()
            if len(parts) != 2: return None
            month_name, year_bs = parts[0], int(parts[1])
            ad_month = self.BS_MONTH_TO_AD.get(month_name)
            if ad_month is None: return None
            ad_year = year_bs - 57 if ad_month >= 4 else year_bs - 56
            return f"{ad_year}-{ad_month:02d}"
        except Exception:
            return None

    @staticmethod
    def _remove_outliers_iqr(series, factor=3.0):
        q1, q3 = series.quantile(0.25), series.quantile(0.75)
        iqr = q3 - q1
        return series.clip(lower=q1 - factor*iqr, upper=q3 + factor*iqr)

    def process_typhoid(self, df):
        df = df.copy()
        typhoid_col = next(c for c in df.columns if "Typhoid" in c)
        df.rename(columns={typhoid_col: "cases",
                           "organisationunitname": "district",
                           "periodname": "period"}, inplace=True)
        df["cases"] = df["cases"].astype(str).str.replace(",", "", regex=False).astype(float)
        df["district"] = df["district"].map(self._normalize_district)
        df["date_ad"]  = df["period"].map(self._bs_to_ad)
        df = df.dropna(subset=["date_ad", "cases"])
        df["cases"] = self._remove_outliers_iqr(df["cases"])
        return df[["district", "date_ad", "cases"]]

    def process_climate(self, df):
        df = df.copy()
        df.rename(columns=COL_MAP, inplace=True)
        df["district"] = df["district"].map(self._normalize_district)
        df["date_ad"] = df["period"].astype(str).str.replace(r"(\d{4})(\d{2})", r"\1-\2", regex=True)
        return df[["district", "date_ad", "temp_min", "temp_max", "temp_mean", "precip", "humidity"]]

    def process_flood(self, df):
        df = df.copy()
        df["district"] = df["district"].map(self._normalize_district)
        df["date_ad"] = df["period"].astype(str).str.replace(r"(\d{4})(\d{2})\d{2}", r"\1-\2", regex=True)
        return (df.groupby(["district", "date_ad"]).size().reset_index(name="flood_events"))

    def merge_datasets(self, typhoid_df, climate_df, flood_df):
        merged = pd.merge(climate_df, typhoid_df, on=["district", "date_ad"], how="inner")
        merged = pd.merge(merged, flood_df, on=["district", "date_ad"], how="left")
        merged["flood_events"] = merged["flood_events"].fillna(0)
        return merged.sort_values(["district", "date_ad"]).reset_index(drop=True)

A.3 engineer.py — lag, rolling, cyclical features

Code 
"""engineer.py — focused, non-redundant feature set."""

import numpy as np
import pandas as pd

FEATURE_COLS = [
    "precip_lag1", "temp_mean_lag1", "humidity_lag1", "flood_lag1",
    "precip_roll3", "temp_roll3",
    "monsoon", "month_sin", "month_cos",
    "cases_lag1",
]
TARGET_COL = "log_cases"
TARGET_RAW = "cases"


class FeatureEngineer:
    def create_features(self, df):
        df = df.copy()
        df["date_dt"] = pd.to_datetime(df["date_ad"])
        df = df.sort_values(["district", "date_dt"]).reset_index(drop=True)

        g = df.groupby("district", sort=False)
        df["precip_lag1"]    = g["precip"].shift(1)
        df["temp_mean_lag1"] = g["temp_mean"].shift(1)
        df["humidity_lag1"]  = g["humidity"].shift(1)
        df["flood_lag1"]     = g["flood_events"].shift(1)
        df["cases_lag1"]     = g["cases"].shift(1)

        df["precip_roll3"] = g["precip"].transform(lambda x: x.shift(1).rolling(3).mean())
        df["temp_roll3"]   = g["temp_mean"].transform(lambda x: x.shift(1).rolling(3).mean())

        df["month"]     = df["date_dt"].dt.month
        df["monsoon"]   = df["month"].between(6, 9).astype(int)
        df["month_sin"] = np.sin(2 * np.pi * df["month"] / 12)
        df["month_cos"] = np.cos(2 * np.pi * df["month"] / 12)

        df[TARGET_COL] = np.log1p(df[TARGET_RAW])
        return df.dropna(subset=FEATURE_COLS + [TARGET_COL]).reset_index(drop=True)

    @staticmethod
    def get_feature_list():
        return list(FEATURE_COLS)

A.4 models.py — Random Forest, XGBoost, MLP

Code 
"""models.py — three predictive architectures."""

import pickle
import numpy as np
import pandas as pd

from sklearn.ensemble import RandomForestRegressor
from sklearn.neural_network import MLPRegressor
from sklearn.preprocessing import StandardScaler
from xgboost import XGBRegressor

try:
    import tensorflow as tf
    _KERAS = tf.keras
    HAS_TF = True
except Exception:
    try:
        import keras as _KERAS; HAS_TF = True
    except Exception:
        HAS_TF = False


class TyphoidModels:
    def __init__(self):
        self.rf = RandomForestRegressor(n_estimators=300, max_depth=10,
                                        min_samples_leaf=5, random_state=42, n_jobs=-1)
        self.xgb = XGBRegressor(n_estimators=300, learning_rate=0.05, max_depth=4,
                                subsample=0.8, colsample_bytree=0.8,
                                reg_alpha=0.1, reg_lambda=1.0,
                                random_state=42, verbosity=0)
        self.scaler_X = StandardScaler()
        self.scaler_y = StandardScaler()
        self.seq_model = None

    def train_rf(self, X, y):  self.rf.fit(X, y);  return self.rf
    def predict_rf(self, X):    return self.rf.predict(X)

    def train_xgb(self, X, y, X_val=None, y_val=None):
        eval_set = [(X_val, y_val)] if X_val is not None else []
        self.xgb.fit(X, y, eval_set=eval_set, verbose=False)
        return self.xgb

    def predict_xgb(self, X):  return self.xgb.predict(X)

    def _build_mlp(self):
        return MLPRegressor(hidden_layer_sizes=(128, 64, 32),
                            activation="relu", solver="adam",
                            learning_rate_init=0.001, max_iter=500,
                            early_stopping=True, validation_fraction=0.1,
                            random_state=42)

    def train_sequential(self, X, y, n_steps=3, epochs=50):
        X_sc = self.scaler_X.fit_transform(X)
        y_sc = self.scaler_y.fit_transform(y.values.reshape(-1, 1)).ravel()
        self.seq_model = self._build_mlp()
        self.seq_model.fit(X_sc, y_sc)
        return self.seq_model

    def predict_sequential(self, X, n_steps=3):
        X_sc = self.scaler_X.transform(X)
        y_sc = self.seq_model.predict(X_sc)
        return self.scaler_y.inverse_transform(y_sc.reshape(-1, 1)).ravel()

    def save_models(self, out_dir="manuscript/output/analysis"):
        with open(f"{out_dir}/rf_model.pkl",  "wb") as f: pickle.dump(self.rf, f)
        with open(f"{out_dir}/xgb_model.pkl", "wb") as f: pickle.dump(self.xgb, f)
        if self.seq_model is not None:
            with open(f"{out_dir}/mlp_model.pkl", "wb") as f: pickle.dump(self.seq_model, f)

A.5 evaluator.py — RMSE · MAE · R² · MAPE

Code 
"""evaluator.py — regression metrics in case-count space."""

import numpy as np
import pandas as pd
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score


class Evaluator:
    def compute_mape(self, y_true, y_pred):
        y_true, y_pred = np.array(y_true), np.array(y_pred)
        mask = y_true > 0
        if mask.sum() == 0: return np.nan
        return float(np.mean(np.abs((y_true[mask] - y_pred[mask]) / y_true[mask])) * 100)

    def evaluate(self, y_true, y_pred_log, model_name):
        y_true = np.array(y_true)
        y_pred = np.clip(np.expm1(np.array(y_pred_log)), 0, None)
        return {
            "Model":    model_name,
            "RMSE":     round(float(np.sqrt(mean_squared_error(y_true, y_pred))), 2),
            "MAE":      round(float(mean_absolute_error(y_true, y_pred)), 2),
            "R2":       round(float(r2_score(y_true, y_pred)), 4),
            "MAPE (%)": round(self.compute_mape(y_true, y_pred), 2),
        }

    def save_metrics(self, results, path="manuscript/output/tables/performance_metrics.csv"):
        df = pd.DataFrame(results)
        df.to_csv(path, index=False)
        return df

A.6 visualizer.py and A.7 reporter.py

Full source for the four research figures (National Trend, Figure 4.2, Figure 4.3, Figure 4.4, Ecological Analysis) and the research-analysis text generator lives at manuscript/src/visualizer.py and manuscript/src/reporter.py.

A.8 main.py — pipeline driver

Code 
"""main.py — research pipeline driver. Run with: python3 manuscript/main.py"""

import os, sys, warnings
import numpy as np, pandas as pd
warnings.filterwarnings("ignore")

DATA_PATH = os.path.abspath(os.path.join(os.path.dirname(__file__), "../../data"))
sys.path.append(os.path.abspath(os.path.join(os.path.dirname(__file__), '../..')))

from src.loader     import DataLoader
from src.processor  import DataProcessor
from src.engineer   import FeatureEngineer, FEATURE_COLS, TARGET_COL, TARGET_RAW
from src.models     import TyphoidModels, HAS_TF
from src.evaluator  import Evaluator
from src.visualizer import Visualizer
from src.reporter   import Reporter

# 1 — load
raw = DataLoader().load_all()

# 2 — preprocess + merge
proc = DataProcessor()
merged = proc.merge_datasets(
    proc.process_typhoid(raw["typhoid"]),
    proc.process_climate(raw["climate"]),
    proc.process_flood(raw["flood"]),
)

# 3 — engineer features
df = FeatureEngineer().create_features(merged)
feats = FEATURE_COLS

# 4 — chronological train/test split (last 12 months held out)
df = df.sort_values("date_dt").reset_index(drop=True)
unique_dates = df["date_dt"].sort_values().unique()
split_date = unique_dates[len(unique_dates) - 12]
train = df[df["date_dt"] <  split_date]
test  = df[df["date_dt"] >= split_date]
X_train, y_train = train[feats], train[TARGET_COL]
X_test,  y_test  = test[feats],  test[TARGET_COL]
y_test_raw = test[TARGET_RAW].values

# 5 — train three models
models = TyphoidModels()
models.train_rf(X_train, y_train); rf_p  = models.predict_rf(X_test)
val_idx = int(len(X_train) * 0.9)
models.train_xgb(X_train.iloc[:val_idx], y_train.iloc[:val_idx],
                 X_train.iloc[val_idx:], y_train.iloc[val_idx:])
xgb_p = models.predict_xgb(X_test)
models.train_sequential(X_train, y_train, n_steps=3, epochs=50)
seq_p = models.predict_sequential(X_test, n_steps=3)

# 6 — evaluate
ev = Evaluator()
results = [ev.evaluate(y_test_raw, p, n) for n, p in [
    ("Random Forest", rf_p), ("XGBoost", xgb_p), ("MLP", seq_p)]]
metrics_df = ev.save_metrics(results)
best = metrics_df.loc[metrics_df["MAE"].idxmin(), "Model"]
print(metrics_df.to_string(index=False))
print(f"\nOperational model (lowest MAE): {best}")

# 7 — visualise + report
viz = Visualizer()
viz.plot_time_series(merged)
viz.plot_correlation(df)
viz.plot_feature_importance(models.xgb, feats)
viz.plot_actual_vs_predicted(y_test_raw,
    {"Random Forest": rf_p, "XGBoost": xgb_p, "MLP": seq_p})
viz.plot_seasonal_boxplot(df)

corr_mat = df[feats + [TARGET_COL]].corr()
Reporter().generate_report(metrics_df, best,
    float(corr_mat.loc["precip_lag1",   TARGET_COL]),
    float(corr_mat.loc["humidity_lag1", TARGET_COL]))
models.save_models()

Run it

terminal
cd manuscript
python3 main.py

Expected console output (excerpted):

        Model    RMSE    MAE      R2  MAPE (%)
Random Forest  131.43  72.10  0.8563    41.48
      XGBoost  129.07  69.50  0.8614    44.64
          MLP  128.07  71.45  0.8635    43.50

Operational model (lowest MAE): XGBoost

Four publication-quality PNGs land in manuscript/output/figures/ and the research-analysis text in manuscript/output/analysis/research_analysis.txt.