Day 2: Demand Forecasting for Antibiotics

WISE Workshop | Addis Ababa, Feb 2026

In this notebook, you’ll build machine learning models to forecast demand for three classes of antibiotics across 100 community health centers. This scenario demonstrates real-world supply chain forecasting with:

• Multiple product classes with different demand patterns
• Facility-level variation based on population served
• Temporal forecasting with a 12-month horizon
• Seasonal effects from disease patterns

Setup

# Import packages
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

from sklearn.model_selection import cross_val_score
from sklearn.preprocessing import LabelEncoder
from sklearn.linear_model import LinearRegression, Ridge
from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score

import warnings
warnings.filterwarnings('ignore')

# Set random seed for reproducibility
np.random.seed(42)

print("Packages loaded!")

Part 1: Generate Synthetic Data

We’ll create realistic synthetic data representing: - 100 community health centers (CHC_001 to CHC_100) - 60 months of historical data (Jan 2019 - Dec 2023) - 3 antibiotic classes with distinct demand patterns: 1. Penicillins (e.g., Amoxicillin) - high volume, seasonal 2. Macrolides (e.g., Azithromycin) - medium volume, respiratory season peaks 3. Fluoroquinolones (e.g., Ciprofloxacin) - lower volume, more stable

# Define parameters
n_facilities = 100
n_months = 60  # Jan 2019 - Dec 2023
antibiotic_classes = ['Penicillins', 'Macrolides', 'Fluoroquinolones']

# Create facility IDs
facility_ids = [f'CHC_{i:03d}' for i in range(1, n_facilities + 1)]

# Assign population served to each facility (varies from 5,000 to 50,000)
facility_populations = {fid: np.random.randint(5000, 50001) for fid in facility_ids}

# Create date range
dates = pd.date_range('2019-01-01', periods=n_months, freq='MS')

print(f"Facilities: {n_facilities}")
print(f"Time period: {dates[0].strftime('%b %Y')} - {dates[-1].strftime('%b %Y')}")
print(f"Antibiotic classes: {antibiotic_classes}")

# Define antibiotic class characteristics
class_params = {
    'Penicillins': {
        'base_demand_per_1000': 15,  # High volume
        'seasonal_amplitude': 0.3,   # Strong seasonality
        'peak_month': 1,             # Peak in January (cold/flu season)
        'trend': 0.02                # 2% annual growth
    },
    'Macrolides': {
        'base_demand_per_1000': 8,   # Medium volume
        'seasonal_amplitude': 0.4,   # Very seasonal (respiratory)
        'peak_month': 12,            # Peak in December
        'trend': 0.03                # 3% annual growth
    },
    'Fluoroquinolones': {
        'base_demand_per_1000': 4,   # Lower volume
        'seasonal_amplitude': 0.1,   # More stable
        'peak_month': 7,             # Slight peak in rainy season
        'trend': 0.01                # 1% annual growth
    }
}

# Define rainy season months (June-September in Ethiopia)
rainy_months = [6, 7, 8, 9]

print("Antibiotic class parameters defined")

# Generate the dataset
data_rows = []

for facility_id in facility_ids:
    population = facility_populations[facility_id]
    
    for date in dates:
        month = date.month
        year = date.year
        
        # Determine season
        season = 'Rainy' if month in rainy_months else 'Dry'
        
        # Simulate occasional disease outbreaks (5% chance per facility-month)
        outbreak = np.random.random() < 0.05
        
        for abx_class in antibiotic_classes:
            params = class_params[abx_class]
            
            # Base demand scaled by population
            base = params['base_demand_per_1000'] * (population / 1000)
            
            # Seasonal effect (sinusoidal with class-specific peak)
            seasonal = params['seasonal_amplitude'] * np.sin(
                2 * np.pi * (month - params['peak_month']) / 12
            )
            
            # Year-over-year trend
            years_from_start = (year - 2019) + (month - 1) / 12
            trend = 1 + params['trend'] * years_from_start
            
            # Outbreak effect (doubles demand for Penicillins and Macrolides)
            outbreak_multiplier = 1.0
            if outbreak and abx_class in ['Penicillins', 'Macrolides']:
                outbreak_multiplier = 2.0
            
            # Calculate demand with noise
            demand = base * (1 + seasonal) * trend * outbreak_multiplier
            noise = np.random.normal(0, demand * 0.15)  # 15% noise
            demand = max(1, int(demand + noise))  # Ensure positive integer
            
            data_rows.append({
                'facility_id': facility_id,
                'date': date,
                'year': year,
                'month': month,
                'antibiotic_class': abx_class,
                'population_served': population,
                'season': season,
                'disease_outbreak': int(outbreak),
                'demand': demand
            })

# Create DataFrame
df = pd.DataFrame(data_rows)

print(f"Dataset shape: {df.shape}")
print(f"Expected: {n_facilities} facilities × {n_months} months × {len(antibiotic_classes)} classes = {n_facilities * n_months * len(antibiotic_classes)} rows")
df.head(10)

# Summary statistics
print("Dataset Summary:")
print(f"\nDate range: {df['date'].min()} to {df['date'].max()}")
print(f"\nFacilities: {df['facility_id'].nunique()}")
print(f"\nAntibiotic classes: {df['antibiotic_class'].unique().tolist()}")
print(f"\nDemand statistics:")
df.groupby('antibiotic_class')['demand'].describe().round(1)

Part 2: Exploratory Data Analysis

# Demand distribution by antibiotic class
fig, axes = plt.subplots(1, 3, figsize=(14, 4))

for ax, abx_class in zip(axes, antibiotic_classes):
    class_data = df[df['antibiotic_class'] == abx_class]['demand']
    ax.hist(class_data, bins=50, edgecolor='black', alpha=0.7)
    ax.set_xlabel('Monthly Demand (units)')
    ax.set_ylabel('Frequency')
    ax.set_title(f'{abx_class}\nMean: {class_data.mean():.0f}, Median: {class_data.median():.0f}')
    ax.axvline(class_data.mean(), color='red', linestyle='--', label='Mean')

plt.tight_layout()
plt.show()

# Seasonal patterns by antibiotic class
monthly_demand = df.groupby(['month', 'antibiotic_class'])['demand'].mean().reset_index()

fig, ax = plt.subplots(figsize=(10, 5))

for abx_class in antibiotic_classes:
    class_data = monthly_demand[monthly_demand['antibiotic_class'] == abx_class]
    ax.plot(class_data['month'], class_data['demand'], marker='o', label=abx_class, linewidth=2)

ax.set_xlabel('Month')
ax.set_ylabel('Average Monthly Demand')
ax.set_title('Seasonal Patterns by Antibiotic Class')
ax.set_xticks(range(1, 13))
ax.set_xticklabels(['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 
                    'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'])
ax.legend()
ax.grid(True, alpha=0.3)

# Highlight rainy season
ax.axvspan(6, 9, alpha=0.2, color='blue', label='Rainy Season')

plt.tight_layout()
plt.show()

# Time series trends
monthly_total = df.groupby(['date', 'antibiotic_class'])['demand'].sum().reset_index()

fig, ax = plt.subplots(figsize=(12, 5))

for abx_class in antibiotic_classes:
    class_data = monthly_total[monthly_total['antibiotic_class'] == abx_class]
    ax.plot(class_data['date'], class_data['demand'], label=abx_class, linewidth=1.5, alpha=0.8)

ax.set_xlabel('Date')
ax.set_ylabel('Total Monthly Demand (all facilities)')
ax.set_title('Demand Trends Over Time')
ax.legend()
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Facility-level variation: demand vs population
facility_summary = df.groupby(['facility_id', 'antibiotic_class']).agg({
    'demand': 'mean',
    'population_served': 'first'
}).reset_index()

fig, ax = plt.subplots(figsize=(10, 6))

for abx_class in antibiotic_classes:
    class_data = facility_summary[facility_summary['antibiotic_class'] == abx_class]
    ax.scatter(class_data['population_served'], class_data['demand'], 
               alpha=0.6, label=abx_class, s=30)

ax.set_xlabel('Population Served')
ax.set_ylabel('Average Monthly Demand')
ax.set_title('Demand vs. Facility Population by Antibiotic Class')
ax.legend()
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

Part 3: Feature Engineering

# Create a copy for feature engineering
df_model = df.copy()

# Time features
df_model['quarter'] = df_model['date'].dt.quarter

# Cyclical encoding for month (captures circular nature of months)
df_model['month_sin'] = np.sin(2 * np.pi * df_model['month'] / 12)
df_model['month_cos'] = np.cos(2 * np.pi * df_model['month'] / 12)

# Encode categorical variables
le_facility = LabelEncoder()
le_class = LabelEncoder()
le_season = LabelEncoder()

df_model['facility_encoded'] = le_facility.fit_transform(df_model['facility_id'])
df_model['class_encoded'] = le_class.fit_transform(df_model['antibiotic_class'])
df_model['season_encoded'] = le_season.fit_transform(df_model['season'])

print("Features created:")
print(df_model[['facility_id', 'month', 'month_sin', 'month_cos', 'quarter', 
                'antibiotic_class', 'class_encoded', 'season', 'season_encoded']].head())

# Create lag features (previous month demand)
# Sort by facility, class, and date first
df_model = df_model.sort_values(['facility_id', 'antibiotic_class', 'date'])

# Create lag features within each facility-class group
df_model['demand_lag1'] = df_model.groupby(['facility_id', 'antibiotic_class'])['demand'].shift(1)
df_model['demand_lag3'] = df_model.groupby(['facility_id', 'antibiotic_class'])['demand'].shift(3)

# Rolling average (3-month)
df_model['demand_rolling_3'] = df_model.groupby(['facility_id', 'antibiotic_class'])['demand'].transform(
    lambda x: x.shift(1).rolling(window=3, min_periods=1).mean()
)

# Drop rows with missing lag values (first few months per facility-class)
df_model = df_model.dropna(subset=['demand_lag1', 'demand_lag3'])

print(f"Shape after adding lag features: {df_model.shape}")
print(f"Rows removed due to lag: {len(df) - len(df_model)}")

# Define feature columns
feature_cols = [
    'month', 'month_sin', 'month_cos', 'quarter', 'year',
    'facility_encoded', 'class_encoded', 'season_encoded',
    'population_served', 'disease_outbreak',
    'demand_lag1', 'demand_lag3', 'demand_rolling_3'
]

print("Feature columns:")
for col in feature_cols:
    print(f"  - {col}")

Part 4: Train/Test Split (Temporal)

Important: For time series forecasting, we use a temporal split to prevent data leakage: - Training data: 2019-2022 (48 months) - Test data: 2023 (12 months) - simulates forecasting the next year

# Temporal train/test split
train_data = df_model[df_model['year'] <= 2022]
test_data = df_model[df_model['year'] == 2023]

X_train = train_data[feature_cols]
y_train = train_data['demand']

X_test = test_data[feature_cols]
y_test = test_data['demand']

print(f"Training period: {train_data['date'].min().strftime('%b %Y')} - {train_data['date'].max().strftime('%b %Y')}")
print(f"Training samples: {len(X_train):,}")
print(f"\nTest period: {test_data['date'].min().strftime('%b %Y')} - {test_data['date'].max().strftime('%b %Y')}")
print(f"Test samples: {len(X_test):,}")

Part 5: Train Models

# Initialize models
models = {
    'Linear Regression': LinearRegression(),
    'Ridge Regression': Ridge(alpha=1.0),
    'Random Forest': RandomForestRegressor(n_estimators=100, max_depth=15, random_state=42, n_jobs=-1),
    'Gradient Boosting': GradientBoostingRegressor(n_estimators=100, max_depth=5, random_state=42)
}

# Train and evaluate each model
results = []
predictions = {}

for name, model in models.items():
    print(f"Training {name}...")
    model.fit(X_train, y_train)
    
    # Predictions
    y_pred = model.predict(X_test)
    predictions[name] = y_pred
    
    # Metrics
    rmse = np.sqrt(mean_squared_error(y_test, y_pred))
    mae = mean_absolute_error(y_test, y_pred)
    r2 = r2_score(y_test, y_pred)
    
    results.append({
        'Model': name,
        'RMSE': rmse,
        'MAE': mae,
        'R²': r2
    })
    
    print(f"  RMSE: {rmse:.2f}, MAE: {mae:.2f}, R²: {r2:.3f}")

print("\nAll models trained!")

Part 6: Model Comparison

# Results summary
results_df = pd.DataFrame(results)
results_df = results_df.sort_values('R²', ascending=False)

print("Model Comparison (sorted by R²):")
display(results_df)

# Visualize model performance
fig, axes = plt.subplots(1, 3, figsize=(14, 4))

metrics = ['RMSE', 'MAE', 'R²']
colors = ['steelblue', 'darkorange', 'green']

for ax, metric, color in zip(axes, metrics, colors):
    bars = ax.barh(results_df['Model'], results_df[metric], color=color, alpha=0.7)
    ax.set_xlabel(metric)
    ax.set_title(f'{metric} by Model')
    
    # Add value labels
    for bar, val in zip(bars, results_df[metric]):
        ax.text(val, bar.get_y() + bar.get_height()/2, f'{val:.2f}', 
                va='center', ha='left', fontsize=9)

plt.tight_layout()
plt.show()

# Predictions vs Actual for each model
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
axes = axes.flatten()

for ax, (name, y_pred) in zip(axes, predictions.items()):
    ax.scatter(y_test, y_pred, alpha=0.3, s=10)
    
    # Perfect prediction line
    max_val = max(y_test.max(), max(y_pred))
    ax.plot([0, max_val], [0, max_val], 'r--', lw=2, label='Perfect Prediction')
    
    ax.set_xlabel('Actual Demand')
    ax.set_ylabel('Predicted Demand')
    ax.set_title(f'{name}')
    ax.legend(loc='upper left')

plt.tight_layout()
plt.show()

Part 7: Error Analysis by Antibiotic Class

# Use best model (Gradient Boosting typically)
best_model_name = results_df.iloc[0]['Model']
best_pred = predictions[best_model_name]

# Add predictions to test data
test_data_eval = test_data.copy()
test_data_eval['predicted'] = best_pred
test_data_eval['error'] = test_data_eval['predicted'] - test_data_eval['demand']
test_data_eval['abs_error'] = np.abs(test_data_eval['error'])
test_data_eval['pct_error'] = 100 * test_data_eval['abs_error'] / test_data_eval['demand']

# Metrics by antibiotic class
class_metrics = test_data_eval.groupby('antibiotic_class').agg({
    'demand': 'mean',
    'abs_error': 'mean',
    'pct_error': 'mean'
}).round(2)
class_metrics.columns = ['Avg Demand', 'MAE', 'MAPE (%)']

print(f"Error Analysis by Antibiotic Class ({best_model_name}):")
display(class_metrics)

# Error distribution by class
fig, axes = plt.subplots(1, 3, figsize=(14, 4))

for ax, abx_class in zip(axes, antibiotic_classes):
    class_errors = test_data_eval[test_data_eval['antibiotic_class'] == abx_class]['error']
    ax.hist(class_errors, bins=40, edgecolor='black', alpha=0.7)
    ax.axvline(0, color='red', linestyle='--', linewidth=2)
    ax.set_xlabel('Prediction Error')
    ax.set_ylabel('Frequency')
    ax.set_title(f'{abx_class}\nMean Error: {class_errors.mean():.1f}')

plt.tight_layout()
plt.show()

Part 8: Feature Importance

# Feature importance from Random Forest and Gradient Boosting
rf_importance = pd.DataFrame({
    'Feature': feature_cols,
    'Random Forest': models['Random Forest'].feature_importances_
})

gb_importance = pd.DataFrame({
    'Feature': feature_cols,
    'Gradient Boosting': models['Gradient Boosting'].feature_importances_
})

importance_df = rf_importance.merge(gb_importance, on='Feature')
importance_df['Average'] = (importance_df['Random Forest'] + importance_df['Gradient Boosting']) / 2
importance_df = importance_df.sort_values('Average', ascending=False)

print("Feature Importance:")
display(importance_df.round(4))

# Visualize feature importance
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# Random Forest
rf_sorted = importance_df.sort_values('Random Forest', ascending=True)
axes[0].barh(rf_sorted['Feature'], rf_sorted['Random Forest'], color='steelblue', alpha=0.7)
axes[0].set_xlabel('Importance')
axes[0].set_title('Random Forest Feature Importance')

# Gradient Boosting
gb_sorted = importance_df.sort_values('Gradient Boosting', ascending=True)
axes[1].barh(gb_sorted['Feature'], gb_sorted['Gradient Boosting'], color='darkorange', alpha=0.7)
axes[1].set_xlabel('Importance')
axes[1].set_title('Gradient Boosting Feature Importance')

plt.tight_layout()
plt.show()

Part 9: Forecast Visualization

# Aggregate predictions by month and antibiotic class
forecast_summary = test_data_eval.groupby(['date', 'antibiotic_class']).agg({
    'demand': 'sum',
    'predicted': 'sum'
}).reset_index()

# Also get historical data for context
historical_summary = train_data.groupby(['date', 'antibiotic_class'])['demand'].sum().reset_index()

print("Forecast summary prepared")

# Plot 12-month forecast vs actual for each antibiotic class
fig, axes = plt.subplots(3, 1, figsize=(12, 10))

for ax, abx_class in zip(axes, antibiotic_classes):
    # Historical data (last 2 years for context)
    hist_class = historical_summary[
        (historical_summary['antibiotic_class'] == abx_class) & 
        (historical_summary['date'] >= '2021-01-01')
    ]
    
    # Forecast period
    forecast_class = forecast_summary[forecast_summary['antibiotic_class'] == abx_class]
    
    # Plot historical
    ax.plot(hist_class['date'], hist_class['demand'], 
            color='steelblue', linewidth=2, label='Historical')
    
    # Plot actual 2023
    ax.plot(forecast_class['date'], forecast_class['demand'], 
            color='green', linewidth=2, label='Actual (2023)', marker='o')
    
    # Plot predicted 2023
    ax.plot(forecast_class['date'], forecast_class['predicted'], 
            color='red', linewidth=2, linestyle='--', label='Forecast', marker='s')
    
    ax.set_ylabel('Total Demand')
    ax.set_title(f'{abx_class}')
    ax.legend(loc='upper left')
    ax.grid(True, alpha=0.3)
    
    # Add vertical line at forecast start
    ax.axvline(pd.Timestamp('2023-01-01'), color='gray', linestyle=':', alpha=0.7)

axes[-1].set_xlabel('Date')
plt.suptitle(f'12-Month Demand Forecast vs Actual ({best_model_name})', fontsize=14, y=1.02)
plt.tight_layout()
plt.show()

# Monthly forecast accuracy
monthly_accuracy = forecast_summary.copy()
monthly_accuracy['error'] = monthly_accuracy['predicted'] - monthly_accuracy['demand']
monthly_accuracy['pct_error'] = 100 * monthly_accuracy['error'] / monthly_accuracy['demand']
monthly_accuracy['month'] = monthly_accuracy['date'].dt.strftime('%b')

# Pivot for display
accuracy_pivot = monthly_accuracy.pivot(index='month', columns='antibiotic_class', values='pct_error')

# Reorder months
month_order = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']
accuracy_pivot = accuracy_pivot.reindex(month_order)

print("Monthly Forecast Error (% over/under actual):")
display(accuracy_pivot.round(1))

Part 10: Summary & Key Takeaways

# Final summary
print("="*60)
print("DEMAND FORECASTING SUMMARY")
print("="*60)
print(f"\nDataset:")
print(f"  - Facilities: 100 Community Health Centers")
print(f"  - Time period: Jan 2019 - Dec 2023 (60 months)")
print(f"  - Antibiotic classes: Penicillins, Macrolides, Fluoroquinolones")
print(f"  - Total records: {len(df):,}")

print(f"\nBest Performing Model: {best_model_name}")
best_results = results_df[results_df['Model'] == best_model_name].iloc[0]
print(f"  - RMSE: {best_results['RMSE']:.2f}")
print(f"  - MAE: {best_results['MAE']:.2f}")
print(f"  - R²: {best_results['R²']:.3f}")

print(f"\nTop 5 Most Important Features:")
for i, row in importance_df.head(5).iterrows():
    print(f"  {row['Feature']}: {row['Average']:.3f}")

print(f"\nKey Insights:")
print("  - Lag features (previous demand) are strong predictors")
print("  - Population served drives facility-level variation")
print("  - Seasonal patterns differ by antibiotic class")
print("  - Ensemble methods (RF, GB) outperform linear models")

# Practical implications
print("\n" + "="*60)
print("PRACTICAL IMPLICATIONS FOR SUPPLY CHAIN PLANNING")
print("="*60)

print("""
1. INVENTORY PLANNING
   - Use 12-month forecasts to plan annual procurement
   - Account for seasonal peaks (especially for Penicillins and Macrolides)
   - Buffer stock during outbreak-prone periods

2. FACILITY-LEVEL ALLOCATION
   - Scale distribution by population served
   - Use rolling 3-month averages to smooth allocation
   - Monitor facilities with high forecast errors

3. CLASS-SPECIFIC STRATEGIES
   - Penicillins: High volume, plan for January peak
   - Macrolides: Strong seasonality, respiratory season focus
   - Fluoroquinolones: More stable, can use longer procurement cycles

4. MODEL UPDATES
   - Retrain models quarterly with new data
   - Monitor for emerging resistance patterns
   - Incorporate disease surveillance data when available
""")

---

Next Steps: - Experiment with additional features (weather, disease surveillance) - Try more advanced models (XGBoost, LSTM for time series) - Build confidence intervals for forecasts - Deploy model for real-time demand monitoring