XGBoost Biologics Usage Forecasting Notebook

This notebook prepares monthly biologics usage data, engineers time-series features, trains XGBoost models for 1-, 3-, and 6-month forecast horizons, evaluates model performance against a naïve baseline, analyses residual bias, applies calibration, and demonstrates a February cut-off forecast.

Notebook flow: data cleaning → aggregation → feature engineering → model training → evaluation → diagnostics → calibration → practical cut-off demonstration.

# Import libraries

import shap
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.model_selection import TimeSeriesSplit, GridSearchCV
from sklearn.metrics import (
    mean_absolute_error,
    mean_squared_error,
    mean_absolute_percentage_error,
    r2_score,
    make_scorer
)

from xgboost import XGBRegressor
from dieboldmariano import dm_test

1. Data loading

# Read in the data

df_xdi = pd.read_csv("xdi_data.csv")

df_xdi.head(25)

	HB	issue_date	location	costcentre	quantity	issue_type	drugfull	drug_packsize	quantity_container	quantity_doseunits	issue_time	homecare_issue
0	CTM	2021-10-18	loc_fb19f988d57a	cc_67da5a161ab4	5 x 1 Vial Pack	ONESTOP	drug_7eae036aa176	1 Vial Pack	5	0	16:30:57	0
1	CTM	2021-10-06	loc_fb19f988d57a	cc_67da5a161ab4	4 x 1 Pre-Filled Syringe Pack	ONESTOP	drug_a3ffd5cb15c9	1 Pre-Filled Syringe Pack	4	0	14:41:07	0
2	CTM	2021-11-17	loc_fb19f988d57a	cc_67da5a161ab4	4 x 1 Pre-Filled Syringe Pack	ONESTOP	drug_a3ffd5cb15c9	1 Pre-Filled Syringe Pack	4	0	09:32:44	0
3	CTM	2021-11-24	loc_fb19f988d57a	cc_67da5a161ab4	4 x 1 Pre-Filled Syringe Pack	ONESTOP	drug_a3ffd5cb15c9	1 Pre-Filled Syringe Pack	4	0	08:56:09	0
4	CTM	2021-12-14	loc_fb19f988d57a	cc_67da5a161ab4	5 x 1 Vial Pack	ONESTOP	drug_7eae036aa176	1 Vial Pack	5	0	09:37:09	0
5	CTM	2021-12-20	loc_fb19f988d57a	cc_67da5a161ab4	4 x 1 Pre-Filled Syringe Pack	OUT	drug_a3ffd5cb15c9	1 Pre-Filled Syringe Pack	4	0	13:45:17	0
6	CTM	2021-12-06	loc_fb19f988d57a	cc_67da5a161ab4	5 x 1 Vial Pack	ONESTOP	drug_7eae036aa176	1 Vial Pack	5	0	15:25:53	0
7	CTM	2022-01-14	loc_fb19f988d57a	cc_67da5a161ab4	2 x 1 Pre-Filled Syringe Pack	ONESTOP	drug_a3ffd5cb15c9	1 Pre-Filled Syringe Pack	2	0	14:50:57	0
8	CTM	2022-02-11	loc_fb19f988d57a	cc_67da5a161ab4	6 x 1 Vial Pack	ONESTOP	drug_7eae036aa176	1 Vial Pack	6	0	10:05:10	0
9	CTM	2022-02-11	loc_fb19f988d57a	cc_67da5a161ab4	2 x 1 Pre-Filled Syringe Pack	ONESTOP	drug_a3ffd5cb15c9	1 Pre-Filled Syringe Pack	2	0	11:56:00	0
10	CTM	2022-03-14	loc_fb19f988d57a	cc_67da5a161ab4	6 x 1 Vial Pack	ONESTOP	drug_7eae036aa176	1 Vial Pack	6	0	10:11:00	0
11	CTM	2022-03-23	loc_fb19f988d57a	cc_67da5a161ab4	2 x 1 Pre-Filled Syringe Pack	ONESTOP	drug_a3ffd5cb15c9	1 Pre-Filled Syringe Pack	2	0	16:50:31	0
12	CTM	2021-09-28	loc_fb19f988d57a	cc_67da5a161ab4	4 x 1 Pre-Filled Syringe Pack	ONESTOP	drug_a3ffd5cb15c9	1 Pre-Filled Syringe Pack	4	0	10:30:49	0
13	CTM	2022-10-26	loc_fb19f988d57a	cc_67da5a161ab4	2 x 1 Pre-Filled Syringe Pack	ONESTOP	drug_a3ffd5cb15c9	1 Pre-Filled Syringe Pack	2	0	14:55:37	0
14	CTM	2022-10-28	loc_fb19f988d57a	cc_67da5a161ab4	12 x 1 Pre-Filled Syringe Pack	OUT	drug_a3ffd5cb15c9	1 Pre-Filled Syringe Pack	12	0	15:16:16	0
15	CTM	2022-11-18	loc_fb19f988d57a	cc_67da5a161ab4	2 x 1 Pre-Filled Syringe Pack	ONESTOP	drug_a3ffd5cb15c9	1 Pre-Filled Syringe Pack	2	0	09:23:06	0
16	CTM	2022-11-18	loc_fb19f988d57a	cc_67da5a161ab4	5 x 1 Vial Pack	ONESTOP	drug_7eae036aa176	1 Vial Pack	5	0	09:18:56	0
17	CTM	2022-12-28	loc_fb19f988d57a	cc_67da5a161ab4	2 x 1 Pre-Filled Syringe Pack	ONESTOP	drug_a3ffd5cb15c9	1 Pre-Filled Syringe Pack	2	0	10:32:14	0
18	CTM	2023-02-10	loc_fb19f988d57a	cc_67da5a161ab4	5 x 1 Vial Pack	ONESTOP	drug_7eae036aa176	1 Vial Pack	5	0	13:25:09	0
19	CTM	2023-02-28	loc_fb19f988d57a	cc_67da5a161ab4	2 x 1 Pre-Filled Syringe Pack	ONESTOP	drug_a3ffd5cb15c9	1 Pre-Filled Syringe Pack	2	0	15:54:10	0
20	CTM	2023-02-28	loc_fb19f988d57a	cc_67da5a161ab4	5 x 1 Vial Pack	ONESTOP	drug_7eae036aa176	1 Vial Pack	5	0	16:05:00	0
21	CTM	2023-03-13	loc_fb19f988d57a	cc_67da5a161ab4	14 x 1 Pre-Filled Syringe Pack	OUT	drug_a3ffd5cb15c9	1 Pre-Filled Syringe Pack	14	0	11:23:41	0
22	CTM	2023-03-23	loc_fb19f988d57a	cc_67da5a161ab4	4 x 1 Vial Pack	INP	drug_7eae036aa176	1 Vial Pack	4	0	10:14:37	0
23	CTM	2022-04-21	loc_fb19f988d57a	cc_67da5a161ab4	2 x 1 Pre-Filled Syringe Pack	ONESTOP	drug_a3ffd5cb15c9	1 Pre-Filled Syringe Pack	2	0	08:42:43	0
24	CTM	2022-04-05	loc_fb19f988d57a	cc_67da5a161ab4	6 x 1 Vial Pack	ONESTOP	drug_7eae036aa176	1 Vial Pack	6	0	08:53:04	0

2. Data cleaning and unit derivation

# Extract units per pack from drug_packsize
df_xdi["units_per_pack"] = (
    df_xdi["drug_packsize"]
    .astype(str)
    .str.extract(r"(\d+)")
    .astype(float)
)

# Clean quantity fields here
def clean_quantity(series):
    cleaned = (
        series
        .astype(str)
        .str.strip()
        .str.replace(",", "", regex=False)
        .str.extract(r"(-?\d+\.?\d*)")[0]
    )

    return pd.to_numeric(cleaned, errors="coerce")

df_xdi["quantity_container_raw"] = df_xdi["quantity_container"]
df_xdi["quantity_doseunits_raw"] = df_xdi["quantity_doseunits"]

df_xdi["quantity_container"] = clean_quantity(df_xdi["quantity_container"])
df_xdi["quantity_doseunits"] = clean_quantity(df_xdi["quantity_doseunits"])

df_xdi["quantity_container"] = df_xdi["quantity_container"].fillna(0)
df_xdi["quantity_doseunits"] = df_xdi["quantity_doseunits"].fillna(0)

df_xdi["total_units_used"] = (
    df_xdi["quantity_container"] * df_xdi["units_per_pack"]
    + df_xdi["quantity_doseunits"]
)

print("Quantity source counts:")
df_xdi["quantity_source"] = np.select(
    [
        (df_xdi["quantity_container"] > 0) & (df_xdi["quantity_doseunits"] == 0),
        (df_xdi["quantity_container"] == 0) & (df_xdi["quantity_doseunits"] > 0),
        (df_xdi["quantity_container"] > 0) & (df_xdi["quantity_doseunits"] > 0),
        (df_xdi["quantity_container"] == 0) & (df_xdi["quantity_doseunits"] == 0),
    ],
    [
        "container_only",
        "doseunit_only",
        "both",
        "neither"
    ],
    default="unknown"
)

print(df_xdi["quantity_source"].value_counts())
print("\nZero total_units_used percentage:")
print((df_xdi["total_units_used"] == 0).mean() * 100)

display(
    df_xdi[
        [
            "quantity_container_raw",
            "quantity_container",
            "quantity_doseunits_raw",
            "quantity_doseunits",
            "drug_packsize",
            "units_per_pack",
            "total_units_used",
            "quantity_source"
        ]
    ].head(100)
)

Quantity source counts:
quantity_source
both              8159
container_only    1290
doseunit_only      977
neither            145
Name: count, dtype: int64

Zero total_units_used percentage:
1.3716772301579794

	quantity_container_raw	quantity_container	quantity_doseunits_raw	quantity_doseunits	drug_packsize	units_per_pack	total_units_used	quantity_source
0	5	5	0	0	1 Vial Pack	1.0	5.0	container_only
1	4	4	0	0	1 Pre-Filled Syringe Pack	1.0	4.0	container_only
2	4	4	0	0	1 Pre-Filled Syringe Pack	1.0	4.0	container_only
3	4	4	0	0	1 Pre-Filled Syringe Pack	1.0	4.0	container_only
4	5	5	0	0	1 Vial Pack	1.0	5.0	container_only
...	...	...	...	...	...	...	...	...
95	2	2	0	0	1 Pre-Filled Syringe Pack	1.0	2.0	container_only
96	2	2	0	0	1 Pre-Filled Syringe Pack	1.0	2.0	container_only
97	2	2	0	0	1 Pre-Filled Syringe Pack	1.0	2.0	container_only
98	2	2	0	0	1 Pre-Filled Syringe Pack	1.0	2.0	container_only
99	2	2	0	0	1 Pre-Filled Syringe Pack	1.0	2.0	container_only

100 rows × 8 columns

3. Date handling and monthly aggregation

# Ensure issue_date is datetime

df_xdi["issue_date"] = pd.to_datetime(
    df_xdi["issue_date"],
    errors="coerce"
)

# Remove rows with missing dates
df_xdi = df_xdi.dropna(subset=["issue_date"])

# Create weekly period

df_xdi["month_period"] = (
    df_xdi["issue_date"]
    .dt.to_period("M")
    .dt.start_time
)

# Ensure homecare_issue is numeric

df_xdi["homecare_issue"] = pd.to_numeric(
    df_xdi["homecare_issue"],
    errors="coerce"
).fillna(0)

# Aggregate the data by drug and month

monthly_units = (
    df_xdi
    .groupby(["drugfull" , "month_period"], as_index=False)
    .agg(
        homecare_units=(
            "total_units_used",
            lambda x: x[df_xdi.loc[x.index, "homecare_issue"] == 1].sum()
        ),
        hospital_units=(
            "total_units_used",
            lambda x: x[df_xdi.loc[x.index, "homecare_issue"] == 0].sum()
        )
    )
)

monthly_units["total_units"] = (
    monthly_units["homecare_units"]
    + monthly_units["hospital_units"]
)

monthly_units.head(25)

	drugfull	month_period	homecare_units	hospital_units	total_units
0	drug_0050f02cadea	2021-09-01	33.0	0.0	33.0
1	drug_0050f02cadea	2021-10-01	52.0	0.0	52.0
2	drug_0050f02cadea	2021-11-01	44.0	0.0	44.0
3	drug_0050f02cadea	2021-12-01	64.0	0.0	64.0
4	drug_0050f02cadea	2022-01-01	40.0	0.0	40.0
5	drug_0050f02cadea	2022-02-01	32.0	0.0	32.0
6	drug_0050f02cadea	2022-03-01	48.0	0.0	48.0
7	drug_0050f02cadea	2022-04-01	48.0	0.0	48.0
8	drug_0050f02cadea	2022-05-01	48.0	0.0	48.0
9	drug_0050f02cadea	2022-06-01	30.0	0.0	30.0
10	drug_0050f02cadea	2022-07-01	40.0	0.0	40.0
11	drug_0050f02cadea	2022-08-01	56.0	0.0	56.0
12	drug_0050f02cadea	2022-09-01	84.0	0.0	84.0
13	drug_0050f02cadea	2022-10-01	4.0	0.0	4.0
14	drug_0050f02cadea	2022-11-01	66.0	0.0	66.0
15	drug_0050f02cadea	2022-12-01	86.0	0.0	86.0
16	drug_0050f02cadea	2023-01-01	40.0	0.0	40.0
17	drug_0050f02cadea	2023-02-01	44.0	0.0	44.0
18	drug_0050f02cadea	2023-03-01	72.0	0.0	72.0
19	drug_0050f02cadea	2023-04-01	62.0	0.0	62.0
20	drug_0050f02cadea	2023-05-01	83.0	0.0	83.0
21	drug_0050f02cadea	2023-06-01	92.0	0.0	92.0
22	drug_0050f02cadea	2023-07-01	57.0	0.0	57.0
23	drug_0050f02cadea	2023-08-01	97.0	0.0	97.0
24	drug_0050f02cadea	2023-09-01	104.0	0.0	104.0

4. Monthly panel creation

# Create complete monthly panel so missing drug/month combinations are treated as zero usage

all_months = pd.date_range(
    monthly_units["month_period"].min(),
    monthly_units["month_period"].max(),
    freq="MS"
)

drug_list = (
    monthly_units[["drugfull"]]
    .drop_duplicates()
)

complete_index = (
    drug_list
    .merge(
        pd.DataFrame({"month_period": all_months}),
        how="cross"
    )
)

monthly_units = (
    complete_index
    .merge(
        monthly_units,
        on=["drugfull", "month_period"],
        how="left"
    )
)

monthly_units["homecare_units"] = monthly_units["homecare_units"].fillna(0)
monthly_units["hospital_units"] = monthly_units["hospital_units"].fillna(0)
monthly_units["total_units"] = monthly_units["total_units"].fillna(0)

monthly_units = monthly_units.sort_values(
    ["drugfull", "month_period"]
).reset_index(drop=True)

# Sort data before creating time-series features

monthly_units = monthly_units.sort_values(
    ["drugfull", "month_period"]
).reset_index(drop=True)

5. Feature engineering

# Create lag features for previous demand

for lag in [1, 2, 3, 4, 5, 6, 12]:
    monthly_units[f"lag_{lag}"] = (
        monthly_units
        .groupby(["drugfull"])["total_units"]
        .shift(lag)
    )

# Add roling mean, std, min, max features

group_cols = ["drugfull"]

monthly_units["used_last_month"] = (
    monthly_units["lag_1"] > 0
).astype(int)

monthly_units["used_last_3m_count"] = (
    monthly_units
    .groupby(group_cols)["total_units"]
    .transform(lambda s: s.shift(1).gt(0).rolling(3).sum())
)

monthly_units["used_last_6m_count"] = (
    monthly_units
    .groupby(group_cols)["total_units"]
    .transform(lambda s: s.shift(1).gt(0).rolling(6).sum())
)

monthly_units["nonzero_rate_6m"] = (
    monthly_units["used_last_6m_count"] / 6
)


for window in [3, 6]:
    monthly_units[f"roll_mean_{window}"] = (
        monthly_units
        .groupby(group_cols)["total_units"]
        .transform(lambda s: s.shift(1).rolling(window).mean())
    )

    monthly_units[f"roll_std_{window}"] = (
        monthly_units
        .groupby(group_cols)["total_units"]
        .transform(lambda s: s.shift(1).rolling(window).std())
    )

    monthly_units[f"roll_min_{window}"] = (
        monthly_units
        .groupby(group_cols)["total_units"]
        .transform(lambda s: s.shift(1).rolling(window).min())
    )

    monthly_units[f"roll_max_{window}"] = (
        monthly_units
        .groupby(group_cols)["total_units"]
        .transform(lambda s: s.shift(1).rolling(window).max())
    )

group_cols = ["drugfull"]

monthly_units["recent_growth_3m"] = (
    monthly_units["roll_mean_3"] /
    monthly_units["roll_mean_6"].replace(0, np.nan)
)

monthly_units["recent_growth_3m"] = (
    monthly_units["recent_growth_3m"]
    .replace([np.inf, -np.inf], np.nan)
    .fillna(0)
)

monthly_units["lag_1_to_mean_6_ratio"] = (
    monthly_units["lag_1"] /
    monthly_units["roll_mean_6"].replace(0, np.nan)
)

monthly_units["lag_1_to_mean_6_ratio"] = (
    monthly_units["lag_1_to_mean_6_ratio"]
    .replace([np.inf, -np.inf], np.nan)
    .fillna(0)
)

display(monthly_units)

	drugfull	month_period	homecare_units	hospital_units	total_units	lag_1	lag_2	lag_3	lag_4	lag_5	...	roll_mean_3	roll_std_3	roll_min_3	roll_max_3	roll_mean_6	roll_std_6	roll_min_6	roll_max_6	recent_growth_3m	lag_1_to_mean_6_ratio
0	drug_0050f02cadea	2021-09-01	33.0	0.0	33.0	NaN	NaN	NaN	NaN	NaN	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	0.000000	0.0
1	drug_0050f02cadea	2021-10-01	52.0	0.0	52.0	33.0	NaN	NaN	NaN	NaN	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	0.000000	0.0
2	drug_0050f02cadea	2021-11-01	44.0	0.0	44.0	52.0	33.0	NaN	NaN	NaN	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	0.000000	0.0
3	drug_0050f02cadea	2021-12-01	64.0	0.0	64.0	44.0	52.0	33.0	NaN	NaN	...	43.000000	9.539392	33.0	52.0	NaN	NaN	NaN	NaN	0.000000	0.0
4	drug_0050f02cadea	2022-01-01	40.0	0.0	40.0	64.0	44.0	52.0	33.0	NaN	...	53.333333	10.066446	44.0	64.0	NaN	NaN	NaN	NaN	0.000000	0.0
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
2803	drug_f80641253af7	2025-10-01	6.0	0.0	6.0	0.0	6.0	6.0	0.0	0.0	...	4.000000	3.464102	0.0	6.0	3.0	3.286335	0.0	6.0	1.333333	0.0
2804	drug_f80641253af7	2025-11-01	0.0	0.0	0.0	6.0	0.0	6.0	6.0	0.0	...	4.000000	3.464102	0.0	6.0	3.0	3.286335	0.0	6.0	1.333333	2.0
2805	drug_f80641253af7	2025-12-01	0.0	0.0	0.0	0.0	6.0	0.0	6.0	6.0	...	2.000000	3.464102	0.0	6.0	3.0	3.286335	0.0	6.0	0.666667	0.0
2806	drug_f80641253af7	2026-01-01	6.0	0.0	6.0	0.0	0.0	6.0	0.0	6.0	...	2.000000	3.464102	0.0	6.0	3.0	3.286335	0.0	6.0	0.666667	0.0
2807	drug_f80641253af7	2026-02-01	0.0	0.0	0.0	6.0	0.0	0.0	6.0	0.0	...	2.000000	3.464102	0.0	6.0	3.0	3.286335	0.0	6.0	0.666667	2.0

2808 rows × 26 columns

# Create homecare and hospital share features

monthly_units["homecare_share"] = (
    monthly_units["homecare_units"] /
    monthly_units["total_units"].replace(0, np.nan)
)

monthly_units["hospital_share"] = (
    monthly_units["hospital_units"] /
    monthly_units["total_units"].replace(0, np.nan)
)

monthly_units["homecare_share"] = monthly_units["homecare_share"].fillna(0)
monthly_units["hospital_share"] = monthly_units["hospital_share"].fillna(0)

# Month on month change features

monthly_units["change_1m"] = (
    monthly_units
    .groupby(["drugfull"])["total_units"]
    .diff(1)
)

monthly_units["change_3m"] = (
    monthly_units
    .groupby(["drugfull"])["total_units"]
    .diff(3)
)

monthly_units["change_6m"] = (
    monthly_units
    .groupby(["drugfull"])["total_units"]
    .diff(6)
)

# Additional peak and momentum features

monthly_units["lag_1_to_roll_max_6_ratio"] = (
    monthly_units["lag_1"] /
    monthly_units["roll_max_6"].replace(0, np.nan)
)

monthly_units["lag_1_to_roll_max_6_ratio"] = (
    monthly_units["lag_1_to_roll_max_6_ratio"]
    .replace([np.inf, -np.inf], np.nan)
    .fillna(0)
)

monthly_units["lag_1_to_roll_mean_3_ratio"] = (
    monthly_units["lag_1"] /
    monthly_units["roll_mean_3"].replace(0, np.nan)
)

monthly_units["lag_1_to_roll_mean_3_ratio"] = (
    monthly_units["lag_1_to_roll_mean_3_ratio"]
    .replace([np.inf, -np.inf], np.nan)
    .fillna(0)
)

monthly_units["recent_peak_flag"] = (
    monthly_units["lag_1"] >= monthly_units["roll_max_6"]
).astype(int)

monthly_units["positive_growth_1m_flag"] = (
    monthly_units["change_1m"] > 0
).astype(int)

monthly_units["positive_growth_3m_flag"] = (
    monthly_units["change_3m"] > 0
).astype(int)

# Drug-level trend and acceleration features

monthly_units["growth_1m"] = (
    monthly_units["lag_1"] - monthly_units["lag_2"]
)

monthly_units["growth_3m"] = (
    monthly_units["roll_mean_3"] - monthly_units["roll_mean_6"]
)

monthly_units["growth_1m_pct"] = (
    monthly_units["growth_1m"] /
    monthly_units["lag_2"].replace(0, np.nan)
)

monthly_units["growth_3m_pct"] = (
    monthly_units["growth_3m"] /
    monthly_units["roll_mean_6"].replace(0, np.nan)
)

monthly_units["growth_1m_pct"] = (
    monthly_units["growth_1m_pct"]
    .replace([np.inf, -np.inf], np.nan)
    .fillna(0)
)

monthly_units["growth_3m_pct"] = (
    monthly_units["growth_3m_pct"]
    .replace([np.inf, -np.inf], np.nan)
    .fillna(0)
)

monthly_units["usage_acceleration"] = (
    monthly_units["change_1m"] - monthly_units["change_3m"] / 3
)

monthly_units["above_6m_average_flag"] = (
    monthly_units["lag_1"] > monthly_units["roll_mean_6"]
).astype(int)

monthly_units["above_6m_max_flag"] = (
    monthly_units["lag_1"] >= monthly_units["roll_max_6"]
).astype(int)

# Flag drugs showing recent sustained growth

monthly_units["recent_growth_flag"] = (
    (monthly_units["roll_mean_3"] > monthly_units["roll_mean_6"] * 1.2)
    & (monthly_units["lag_1"] > monthly_units["roll_mean_6"])
).astype(int)

# Add calendar-based features

monthly_units["month_number"] = (
    monthly_units["month_period"].dt.month
)

monthly_units["year"] = (
    monthly_units["month_period"].dt.year
)

# Encode categorical fields

monthly_units["drug_code"] = (
    monthly_units["drugfull"]
    .astype("category")
    .cat.codes
)

6. Target variable creation

# Create target variables
# These predict total units 1, 3 and 6 months ahead for each drug and cost centre
group_cols = ["drugfull"]

monthly_units["target_1m"] = (
    monthly_units
    .groupby(group_cols)["total_units"]
    .shift(-1)
)

monthly_units["target_3m"] = (
    monthly_units.groupby(group_cols)["total_units"].shift(-1)
    + monthly_units.groupby(group_cols)["total_units"].shift(-2)
    + monthly_units.groupby(group_cols)["total_units"].shift(-3)
)

monthly_units["target_6m"] = (
    monthly_units.groupby(group_cols)["total_units"].shift(-1)
    + monthly_units.groupby(group_cols)["total_units"].shift(-2)
    + monthly_units.groupby(group_cols)["total_units"].shift(-3)
    + monthly_units.groupby(group_cols)["total_units"].shift(-4)
    + monthly_units.groupby(group_cols)["total_units"].shift(-5)
    + monthly_units.groupby(group_cols)["total_units"].shift(-6)
)

7. Feature list and modelling dataset

# Set features and target

features = [
    # Lag features
    "lag_1",
    "lag_2",
    "lag_3",
    "lag_4",
    "lag_5",
    "lag_6",
    "lag_12",

    # Recent usage features
    "used_last_month",
    "used_last_3m_count",
    "used_last_6m_count",
    "nonzero_rate_6m",

    # Rolling mean features
    "roll_mean_3",
    "roll_mean_6",

    # Rolling variation features
    "roll_std_3",
    "roll_std_6",

    # Rolling minimum / maximum features
    "roll_min_3",
    "roll_min_6",
    "roll_max_3",
    "roll_max_6",
    "recent_growth_3m",
    "lag_1_to_mean_6_ratio",
    "lag_1_to_roll_max_6_ratio",
    "lag_1_to_roll_mean_3_ratio",
    "recent_peak_flag",
    "positive_growth_1m_flag",
    "positive_growth_3m_flag",
    "growth_1m",
    "growth_3m",
    "growth_1m_pct",
    "growth_3m_pct",
    "usage_acceleration",
    "above_6m_average_flag",
    "above_6m_max_flag",
    "recent_growth_flag",

    # Homecare / hospital split
    "homecare_share",
    "hospital_share",

    # Change features
    "change_1m",
    "change_3m",
    "change_6m",

    # Calendar features
    "month_number",
    "year",

    # Encoded categorical fields
    "drug_code",
    
]

targets = {
    "1_month": "target_1m",
    "3_month": "target_3m",
    "6_month": "target_6m"
}

# Remove rows with missing feature values only
# Target-specific missing values are handled separately for each forecast horizon

base_df = monthly_units.dropna(
    subset=features
).reset_index(drop=True)

display(base_df.head(25))

	drugfull	month_period	homecare_units	hospital_units	total_units	lag_1	lag_2	lag_3	lag_4	lag_5	...	usage_acceleration	above_6m_average_flag	above_6m_max_flag	recent_growth_flag	month_number	year	drug_code	target_1m	target_3m	target_6m
0	drug_0050f02cadea	2022-09-01	84.0	0.0	84.0	56.0	40.0	30.0	48.0	48.0	...	10.000000	1	1	0	9	2022	0	4.0	156.0	312.0
1	drug_0050f02cadea	2022-10-01	4.0	0.0	4.0	84.0	56.0	40.0	30.0	48.0	...	-68.000000	1	1	0	10	2022	0	66.0	192.0	370.0
2	drug_0050f02cadea	2022-11-01	66.0	0.0	66.0	4.0	84.0	56.0	40.0	30.0	...	58.666667	0	0	0	11	2022	0	86.0	170.0	387.0
3	drug_0050f02cadea	2022-12-01	86.0	0.0	86.0	66.0	4.0	84.0	56.0	40.0	...	19.333333	1	0	0	12	2022	0	40.0	156.0	393.0
4	drug_0050f02cadea	2023-01-01	40.0	0.0	40.0	86.0	66.0	4.0	84.0	56.0	...	-58.000000	1	1	0	1	2023	0	44.0	178.0	410.0
5	drug_0050f02cadea	2023-02-01	44.0	0.0	44.0	40.0	86.0	66.0	4.0	84.0	...	11.333333	0	0	0	2	2023	0	72.0	217.0	463.0
6	drug_0050f02cadea	2023-03-01	72.0	0.0	72.0	44.0	40.0	86.0	66.0	4.0	...	32.666667	0	0	0	3	2023	0	62.0	237.0	495.0
7	drug_0050f02cadea	2023-04-01	62.0	0.0	62.0	72.0	44.0	40.0	86.0	66.0	...	-17.333333	1	0	0	4	2023	0	83.0	232.0	507.0
8	drug_0050f02cadea	2023-05-01	83.0	0.0	83.0	62.0	72.0	44.0	40.0	86.0	...	8.000000	1	0	0	5	2023	0	92.0	246.0	511.0
9	drug_0050f02cadea	2023-06-01	92.0	0.0	92.0	83.0	62.0	72.0	44.0	40.0	...	2.333333	1	0	0	6	2023	0	57.0	258.0	495.0
10	drug_0050f02cadea	2023-07-01	57.0	0.0	57.0	92.0	83.0	62.0	72.0	44.0	...	-33.333333	1	1	1	7	2023	0	97.0	275.0	512.0
11	drug_0050f02cadea	2023-08-01	97.0	0.0	97.0	57.0	92.0	83.0	62.0	72.0	...	35.333333	0	0	0	8	2023	0	104.0	265.0	546.0
12	drug_0050f02cadea	2023-09-01	104.0	0.0	104.0	97.0	57.0	92.0	83.0	62.0	...	3.000000	1	1	0	9	2023	0	74.0	237.0	512.0
13	drug_0050f02cadea	2023-10-01	74.0	0.0	74.0	104.0	97.0	57.0	92.0	83.0	...	-35.666667	1	1	0	10	2023	0	87.0	237.0	505.0
14	drug_0050f02cadea	2023-11-01	87.0	0.0	87.0	74.0	104.0	97.0	57.0	92.0	...	16.333333	0	0	0	11	2023	0	76.0	281.0	510.0
15	drug_0050f02cadea	2023-12-01	76.0	0.0	76.0	87.0	74.0	104.0	97.0	57.0	...	-1.666667	1	0	0	12	2023	0	74.0	275.0	522.0
16	drug_0050f02cadea	2024-01-01	74.0	0.0	74.0	76.0	87.0	74.0	104.0	97.0	...	-2.000000	0	0	0	1	2024	0	131.0	268.0	632.0
17	drug_0050f02cadea	2024-02-01	131.0	0.0	131.0	74.0	76.0	87.0	74.0	104.0	...	42.333333	0	0	0	2	2024	0	70.0	229.0	535.0
18	drug_0050f02cadea	2024-03-01	70.0	0.0	70.0	131.0	74.0	76.0	87.0	74.0	...	-59.000000	1	1	0	3	2024	0	67.0	247.0	671.0
19	drug_0050f02cadea	2024-04-01	67.0	0.0	67.0	70.0	131.0	74.0	76.0	87.0	...	-0.666667	0	0	0	4	2024	0	92.0	364.0	696.0
20	drug_0050f02cadea	2024-05-01	92.0	0.0	92.0	67.0	70.0	131.0	74.0	76.0	...	38.000000	0	0	0	5	2024	0	88.0	306.0	670.0
21	drug_0050f02cadea	2024-06-01	88.0	0.0	88.0	92.0	67.0	70.0	131.0	74.0	...	-10.000000	1	0	0	6	2024	0	184.0	424.0	693.0
22	drug_0050f02cadea	2024-07-01	184.0	0.0	184.0	88.0	92.0	67.0	70.0	131.0	...	57.000000	1	0	0	7	2024	0	34.0	332.0	615.0
23	drug_0050f02cadea	2024-08-01	32.0	2.0	34.0	184.0	88.0	92.0	67.0	70.0	...	-130.666667	1	1	0	8	2024	0	206.0	364.0	661.0
24	drug_0050f02cadea	2024-09-01	206.0	0.0	206.0	34.0	184.0	88.0	92.0	67.0	...	132.666667	0	0	0	9	2024	0	92.0	269.0	497.0

25 rows × 50 columns

8. XGBoost configuration and hyperparameter tuning

# Define XGBoost model

base_model = XGBRegressor(
    objective="reg:squarederror",
    random_state=42,
    n_jobs=-1
)

# Parameter grid for tuning

param_grid = {
    "n_estimators": [500, 700, 900, 1100],
    "max_depth": [4, 5, 6],
    "learning_rate": [0.02, 0.03, 0.05],
    "min_child_weight": [1, 2, 3],
    "subsample": [0.8, 0.9, 1.0],
    "colsample_bytree": [0.8, 0.9, 1.0],
    "reg_alpha": [0.0, 0.1, 0.25],
    "reg_lambda": [0.5, 1.0, 1.5],
}

# Time-series cross-validation

tscv = TimeSeriesSplit(n_splits=5)

rmse_scorer = make_scorer(
    lambda y_true, y_pred: np.sqrt(mean_squared_error(y_true, y_pred)),
    greater_is_better=False
)

# Run grid search separately for each forecast horizon

grid_searches = {}
best_models = {}
best_params = {}
best_cv_scores = {}
horizon_data = {}

weight_config = {
    "1_month": {"q85": 4, "q90": 8, "q95": 12},
    "3_month": {"q85": 3, "q90": 6, "q95": 10},
    "6_month": {"q85": 2, "q90": 5, "q95": 8},
}

for horizon_name, target_col in targets.items():

    print("\n==============================")
    print(f"Preparing data for {horizon_name}")
    print("==============================")

    # Keep only rows where this horizon's target exists
    horizon_df = base_df.dropna(
        subset=[target_col]
    ).reset_index(drop=True)

    # Split this horizon's data by time
    split_date = horizon_df["month_period"].quantile(0.8)

    train = horizon_df[horizon_df["month_period"] <= split_date].copy()
    valid = horizon_df[horizon_df["month_period"] > split_date].copy()

    X_train = train[features]
    y_train_h = train[target_col]

    X_valid = valid[features]
    y_valid_h = valid[target_col]

    high_volume_scale = y_train_h.quantile(0.90)

    cfg = weight_config[horizon_name]

    q85 = y_train_h.quantile(0.85)
    q90 = y_train_h.quantile(0.90)
    q95 = y_train_h.quantile(0.95)

    sample_weights = np.select(
        [
            y_train_h >= q95,
            y_train_h >= q90,
            y_train_h >= q85
        ],
        [
            cfg["q95"],
            cfg["q90"],
            cfg["q85"]
        ],
        default=1
    )

    print("Sample weight summary:")
    print(pd.Series(sample_weights).describe())

    horizon_data[horizon_name] = {
        "train": train,
        "valid": valid,
        "X_train": X_train,
        "X_valid": X_valid,
        "y_train": y_train_h,
        "y_valid": y_valid_h,
        "sample_weights": sample_weights,
        "target_col": target_col,
        "split_date": split_date
    }

    print("Training rows:", len(train))
    print("Validation rows:", len(valid))
    print("Split date:", split_date)

    grid_search = GridSearchCV(
        estimator=base_model,
        param_grid=param_grid,
        scoring=rmse_scorer,
        cv=tscv,
        verbose=2,
        n_jobs=-1
    )

    grid_search.fit(
        X_train,
        y_train_h,
        sample_weight=sample_weights
    )

    grid_searches[horizon_name] = grid_search
    best_models[horizon_name] = grid_search.best_estimator_
    best_params[horizon_name] = grid_search.best_params_
    best_cv_scores[horizon_name] = -grid_search.best_score_

print("Grid search complete for all horizons.")

==============================
Preparing data for 1_month
==============================
Sample weight summary:
count    1716.000000
mean        2.078089
std         2.836614
min         1.000000
25%         1.000000
50%         1.000000
75%         1.000000
max        12.000000
dtype: float64
Training rows: 1716
Validation rows: 416
Split date: 2025-05-01 00:00:00

c:\Users\th252839\Desktop\Python Projects\Python\Final Project XGBoost\final_project\.venv\Lib\site-packages\sklearn\model_selection\_search.py:883: UserWarning: The scoring make_scorer(<lambda>, greater_is_better=False, response_method='predict') does not support sample_weight, which may lead to statistically incorrect results when fitting GridSearchCV(cv=TimeSeriesSplit(gap=0, max_train_size=None, n_splits=5, test_size=None),
             estimator=XGBRegressor(base_score=None, booster=None,
                                    callbacks=None, colsample_bylevel=None,
                                    colsample_bynode=None,
                                    colsample_bytree=None, device=None,
                                    early_stopping_rounds=None,
                                    enable_categorical=False, eval_metric=None,
                                    feature_types=None, feature_weights=None,
                                    gamma=None...
                                    n_jobs=-1, num_parallel_tree=None, ...),
             n_jobs=-1,
             param_grid={'colsample_bytree': [0.8, 0.9, 1.0],
                         'learning_rate': [0.02, 0.03, 0.05],
                         'max_depth': [4, 5, 6], 'min_child_weight': [1, 2, 3],
                         'n_estimators': [500, 700, 900, 1100],
                         'reg_alpha': [0.0, 0.1, 0.25],
                         'reg_lambda': [0.5, 1.0, 1.5],
                         'subsample': [0.8, 0.9, 1.0]},
             scoring=make_scorer(<lambda>, greater_is_better=False, response_method='predict'),
             verbose=2) with sample_weight. 
  warnings.warn(

Fitting 5 folds for each of 8748 candidates, totalling 43740 fits

==============================
Preparing data for 3_month
==============================
Sample weight summary:
count    1664.000000
mean        1.803486
std         2.211969
min         1.000000
25%         1.000000
50%         1.000000
75%         1.000000
max        10.000000
dtype: float64
Training rows: 1664
Validation rows: 364
Split date: 2025-04-01 00:00:00
Fitting 5 folds for each of 8748 candidates, totalling 43740 fits

c:\Users\th252839\Desktop\Python Projects\Python\Final Project XGBoost\final_project\.venv\Lib\site-packages\sklearn\model_selection\_search.py:883: UserWarning: The scoring make_scorer(<lambda>, greater_is_better=False, response_method='predict') does not support sample_weight, which may lead to statistically incorrect results when fitting GridSearchCV(cv=TimeSeriesSplit(gap=0, max_train_size=None, n_splits=5, test_size=None),
             estimator=XGBRegressor(base_score=None, booster=None,
                                    callbacks=None, colsample_bylevel=None,
                                    colsample_bynode=None,
                                    colsample_bytree=None, device=None,
                                    early_stopping_rounds=None,
                                    enable_categorical=False, eval_metric=None,
                                    feature_types=None, feature_weights=None,
                                    gamma=None...
                                    n_jobs=-1, num_parallel_tree=None, ...),
             n_jobs=-1,
             param_grid={'colsample_bytree': [0.8, 0.9, 1.0],
                         'learning_rate': [0.02, 0.03, 0.05],
                         'max_depth': [4, 5, 6], 'min_child_weight': [1, 2, 3],
                         'n_estimators': [500, 700, 900, 1100],
                         'reg_alpha': [0.0, 0.1, 0.25],
                         'reg_lambda': [0.5, 1.0, 1.5],
                         'subsample': [0.8, 0.9, 1.0]},
             scoring=make_scorer(<lambda>, greater_is_better=False, response_method='predict'),
             verbose=2) with sample_weight. 
  warnings.warn(

==============================
Preparing data for 6_month
==============================
Sample weight summary:
count    1508.000000
mean        1.604111
std         1.721224
min         1.000000
25%         1.000000
50%         1.000000
75%         1.000000
max         8.000000
dtype: float64
Training rows: 1508
Validation rows: 364
Split date: 2025-01-01 00:00:00
Fitting 5 folds for each of 8748 candidates, totalling 43740 fits

c:\Users\th252839\Desktop\Python Projects\Python\Final Project XGBoost\final_project\.venv\Lib\site-packages\sklearn\model_selection\_search.py:883: UserWarning: The scoring make_scorer(<lambda>, greater_is_better=False, response_method='predict') does not support sample_weight, which may lead to statistically incorrect results when fitting GridSearchCV(cv=TimeSeriesSplit(gap=0, max_train_size=None, n_splits=5, test_size=None),
             estimator=XGBRegressor(base_score=None, booster=None,
                                    callbacks=None, colsample_bylevel=None,
                                    colsample_bynode=None,
                                    colsample_bytree=None, device=None,
                                    early_stopping_rounds=None,
                                    enable_categorical=False, eval_metric=None,
                                    feature_types=None, feature_weights=None,
                                    gamma=None...
                                    n_jobs=-1, num_parallel_tree=None, ...),
             n_jobs=-1,
             param_grid={'colsample_bytree': [0.8, 0.9, 1.0],
                         'learning_rate': [0.02, 0.03, 0.05],
                         'max_depth': [4, 5, 6], 'min_child_weight': [1, 2, 3],
                         'n_estimators': [500, 700, 900, 1100],
                         'reg_alpha': [0.0, 0.1, 0.25],
                         'reg_lambda': [0.5, 1.0, 1.5],
                         'subsample': [0.8, 0.9, 1.0]},
             scoring=make_scorer(<lambda>, greater_is_better=False, response_method='predict'),
             verbose=2) with sample_weight. 
  warnings.warn(

Grid search complete for all horizons.

9. Best model parameters

# Display best parameters and CV scores for each forecast horizon

best_summary = pd.DataFrame({
    "Best_CV_RMSE": best_cv_scores,
    "Best_Params": best_params
})

display(best_summary)

for horizon_name in targets.keys():
    print("\n==============================")
    print(f"Best model for {horizon_name}")
    print("==============================")
    print("Best CV RMSE:", best_cv_scores[horizon_name])
    print("Best parameters:")
    print(best_params[horizon_name])

	Best_CV_RMSE	Best_Params
1_month	20.759878	{'colsample_bytree': 0.8, 'learning_rate': 0.0...
3_month	52.252289	{'colsample_bytree': 0.9, 'learning_rate': 0.0...
6_month	103.979847	{'colsample_bytree': 0.8, 'learning_rate': 0.0...

==============================
Best model for 1_month
==============================
Best CV RMSE: 20.759877536912324
Best parameters:
{'colsample_bytree': 0.8, 'learning_rate': 0.03, 'max_depth': 4, 'min_child_weight': 2, 'n_estimators': 900, 'reg_alpha': 0.1, 'reg_lambda': 0.5, 'subsample': 0.8}

==============================
Best model for 3_month
==============================
Best CV RMSE: 52.252289479007246
Best parameters:
{'colsample_bytree': 0.9, 'learning_rate': 0.03, 'max_depth': 4, 'min_child_weight': 3, 'n_estimators': 500, 'reg_alpha': 0.25, 'reg_lambda': 0.5, 'subsample': 0.8}

==============================
Best model for 6_month
==============================
Best CV RMSE: 103.97984653468572
Best parameters:
{'colsample_bytree': 0.8, 'learning_rate': 0.03, 'max_depth': 4, 'min_child_weight': 1, 'n_estimators': 500, 'reg_alpha': 0.25, 'reg_lambda': 1.5, 'subsample': 0.8}

10. Naïve baseline and validation predictions

# Naive forecast baseline for each forecast horizon
# This assumes future demand is the same as the most recent observed month

naive_results = {}

for horizon_name, target_col in targets.items():

    valid = horizon_data[horizon_name]["valid"]

    y_valid_h = valid[target_col]

    if horizon_name == "1_month":
        naive_forecast = valid["total_units"]

    elif horizon_name == "3_month":
        naive_forecast = valid["total_units"] * 3

    elif horizon_name == "6_month":
        naive_forecast = valid["total_units"] * 6

    naive_mae = mean_absolute_error(y_valid_h, naive_forecast)
    naive_rmse = np.sqrt(mean_squared_error(y_valid_h, naive_forecast))

    naive_results[horizon_name] = {
        "Naive_MAE": naive_mae,
        "Naive_RMSE": naive_rmse
    }

display(pd.DataFrame(naive_results).T)

	Naive_MAE	Naive_RMSE
1_month	13.091346	35.492077
3_month	21.343407	43.149147
6_month	45.791209	119.063978

# Predict validation data for each forecast horizon

predictions = {}

for horizon_name, target_col in targets.items():

    model = best_models[horizon_name]

    X_valid = horizon_data[horizon_name]["X_valid"]
    y_valid_h = horizon_data[horizon_name]["y_valid"]

    y_pred_h = model.predict(X_valid)
    y_pred_h = np.clip(y_pred_h, 0, None)

    predictions[horizon_name] = {
        "actual": y_valid_h,
        "predicted": y_pred_h
    }

11. Model evaluation metrics

# Evaluation metrics for each forecast horizon

model_results = {}

for horizon_name, target_col in targets.items():

    y_valid_h = predictions[horizon_name]["actual"]
    y_pred_h = predictions[horizon_name]["predicted"]

    mae = mean_absolute_error(y_valid_h, y_pred_h)
    rmse = np.sqrt(mean_squared_error(y_valid_h, y_pred_h))
    r2 = r2_score(y_valid_h, y_pred_h)

    # MAPE - avoids division by zero
    mape = np.nanmean(
        np.abs((y_valid_h - y_pred_h) / np.where(y_valid_h == 0, np.nan, y_valid_h))
    ) * 100

    # WMAPE - more stable than MAPE for low-volume / zero-heavy demand
    wmape = (
        np.sum(np.abs(y_valid_h - y_pred_h)) /
        np.sum(np.abs(y_valid_h))
    ) * 100

    # SMAPE - symmetric percentage error
    smape = (
        np.nanmean(
            2 * np.abs(y_pred_h - y_valid_h) /
            np.where(
                (np.abs(y_valid_h) + np.abs(y_pred_h)) == 0,
                np.nan,
                (np.abs(y_valid_h) + np.abs(y_pred_h))
            )
        )
    ) * 100

    model_results[horizon_name] = {
        "MAE": mae,
        "RMSE": rmse,
        "MAPE": mape,
        "WMAPE": wmape,
        "SMAPE": smape,
        "R2": r2,
        "Naive_MAE": naive_results[horizon_name]["Naive_MAE"],
        "Naive_RMSE": naive_results[horizon_name]["Naive_RMSE"],
        "Best_CV_RMSE": best_cv_scores[horizon_name]
    }

results_df = pd.DataFrame(model_results).T
results_df["MAE_Improvement_%"] = (
    (results_df["Naive_MAE"] - results_df["MAE"]) /
    results_df["Naive_MAE"]
) * 100

results_df["RMSE_Improvement_%"] = (
    (results_df["Naive_RMSE"] - results_df["RMSE"]) /
    results_df["Naive_RMSE"]
) * 100

display(results_df)

	MAE	RMSE	MAPE	WMAPE	SMAPE	R2	Naive_MAE	Naive_RMSE	Best_CV_RMSE	MAE_Improvement_%	RMSE_Improvement_%
1_month	11.263416	29.872100	68.103716	39.043255	120.363141	0.871214	13.091346	35.492077	20.759878	13.962889	15.834455
3_month	20.190313	51.039271	58.469541	22.710281	98.674976	0.957831	21.343407	43.149147	52.252289	5.402575	-18.285701
6_month	41.336035	87.083867	57.803144	23.829332	94.218482	0.969298	45.791209	119.063978	103.979847	9.729322	26.859602

12. Forecast result tables

# Create results dataframes for each forecast horizon

results_by_horizon = {}

for horizon_name, target_col in targets.items():

    valid = horizon_data[horizon_name]["valid"]

    results_h = valid[[
        "drugfull",
        "month_period",
        "total_units",
        target_col,
        "lag_1"
    ]].copy()

    results_h = results_h.rename(columns={
        "total_units": "current_month_units",
        target_col: "future_actual_usage"
    })

    results_h["xgboost_prediction"] = predictions[horizon_name]["predicted"]
    # Correct naïve forecast depending on forecast horizon
    if horizon_name == "1_month":
        results_h["naive_forecast"] = results_h["current_month_units"]
    elif horizon_name == "3_month":
        results_h["naive_forecast"] = results_h["current_month_units"] * 3
    elif horizon_name == "6_month":
        results_h["naive_forecast"] = results_h["current_month_units"] * 6

    results_h["absolute_error"] = np.abs(
        results_h["future_actual_usage"] - results_h["xgboost_prediction"]
    )

    results_h["squared_error"] = (
        results_h["future_actual_usage"] - results_h["xgboost_prediction"]
    ) ** 2

    results_h["horizon"] = horizon_name

    results_by_horizon[horizon_name] = results_h

# Combine all horizon results into one dataframe
all_forecast_results = pd.concat(
    results_by_horizon.values(),
    ignore_index=True
)

display(all_forecast_results.head(50))

	drugfull	month_period	current_month_units	future_actual_usage	lag_1	xgboost_prediction	naive_forecast	absolute_error	squared_error	horizon
0	drug_0050f02cadea	2025-06-01	79.0	100.0	102.0	80.544350	79.0	19.455650	378.522330	1_month
1	drug_0050f02cadea	2025-07-01	100.0	114.0	79.0	72.962692	100.0	41.037308	1684.060626	1_month
2	drug_0050f02cadea	2025-08-01	114.0	58.0	100.0	98.449707	114.0	40.449707	1636.178799	1_month
3	drug_0050f02cadea	2025-09-01	58.0	120.0	114.0	74.587189	58.0	45.412811	2062.323428	1_month
4	drug_0050f02cadea	2025-10-01	120.0	96.0	58.0	71.184875	120.0	24.815125	615.790405	1_month
5	drug_0050f02cadea	2025-11-01	96.0	99.0	120.0	80.370819	96.0	18.629181	347.046381	1_month
6	drug_0050f02cadea	2025-12-01	99.0	125.0	96.0	91.588509	99.0	33.411491	1116.327757	1_month
7	drug_0050f02cadea	2026-01-01	125.0	48.0	99.0	92.073563	125.0	44.073563	1942.478922	1_month
8	drug_0d3c27ea4824	2025-06-01	0.0	0.0	0.0	1.127236	0.0	1.127236	1.270662	1_month
9	drug_0d3c27ea4824	2025-07-01	0.0	0.0	0.0	0.703001	0.0	0.703001	0.494210	1_month
10	drug_0d3c27ea4824	2025-08-01	0.0	0.0	0.0	1.398705	0.0	1.398705	1.956374	1_month
11	drug_0d3c27ea4824	2025-09-01	0.0	0.0	0.0	0.391127	0.0	0.391127	0.152980	1_month
12	drug_0d3c27ea4824	2025-10-01	0.0	0.0	0.0	0.000000	0.0	0.000000	0.000000	1_month
13	drug_0d3c27ea4824	2025-11-01	0.0	0.0	0.0	0.000000	0.0	0.000000	0.000000	1_month
14	drug_0d3c27ea4824	2025-12-01	0.0	0.0	0.0	2.040641	0.0	2.040641	4.164217	1_month
15	drug_0d3c27ea4824	2026-01-01	0.0	0.0	0.0	0.000000	0.0	0.000000	0.000000	1_month
16	drug_119bcb15803e	2025-06-01	14.0	10.0	16.0	17.645451	14.0	7.645451	58.452915	1_month
17	drug_119bcb15803e	2025-07-01	10.0	20.0	14.0	21.232262	10.0	1.232262	1.518469	1_month
18	drug_119bcb15803e	2025-08-01	20.0	20.0	10.0	24.187265	20.0	4.187265	17.533191	1_month
19	drug_119bcb15803e	2025-09-01	20.0	12.0	20.0	24.831306	20.0	12.831306	164.642425	1_month
20	drug_119bcb15803e	2025-10-01	12.0	10.0	20.0	20.132589	12.0	10.132589	102.669367	1_month
21	drug_119bcb15803e	2025-11-01	10.0	18.0	12.0	16.877401	10.0	1.122599	1.260228	1_month
22	drug_119bcb15803e	2025-12-01	18.0	16.0	10.0	11.292971	18.0	4.707029	22.156125	1_month
23	drug_119bcb15803e	2026-01-01	16.0	6.0	18.0	11.183419	16.0	5.183419	26.867835	1_month
24	drug_139f6537068d	2025-06-01	0.0	0.0	0.0	0.374211	0.0	0.374211	0.140034	1_month
25	drug_139f6537068d	2025-07-01	0.0	0.0	0.0	0.092667	0.0	0.092667	0.008587	1_month
26	drug_139f6537068d	2025-08-01	0.0	0.0	0.0	1.412827	0.0	1.412827	1.996079	1_month
27	drug_139f6537068d	2025-09-01	0.0	0.0	0.0	0.447313	0.0	0.447313	0.200089	1_month
28	drug_139f6537068d	2025-10-01	0.0	0.0	0.0	0.051859	0.0	0.051859	0.002689	1_month
29	drug_139f6537068d	2025-11-01	0.0	3.0	0.0	0.000000	0.0	3.000000	9.000000	1_month
30	drug_139f6537068d	2025-12-01	3.0	0.0	0.0	3.548332	3.0	3.548332	12.590663	1_month
31	drug_139f6537068d	2026-01-01	0.0	0.0	3.0	0.000000	0.0	0.000000	0.000000	1_month
32	drug_1b7500ddafc8	2025-06-01	0.0	20.0	12.0	16.085110	0.0	3.914890	15.326366	1_month
33	drug_1b7500ddafc8	2025-07-01	20.0	0.0	0.0	15.472316	20.0	15.472316	239.392556	1_month
34	drug_1b7500ddafc8	2025-08-01	0.0	40.0	20.0	12.106170	0.0	27.893830	778.065769	1_month
35	drug_1b7500ddafc8	2025-09-01	40.0	16.0	0.0	24.145588	40.0	8.145588	66.350603	1_month
36	drug_1b7500ddafc8	2025-10-01	16.0	16.0	40.0	13.388783	16.0	2.611217	6.818452	1_month
37	drug_1b7500ddafc8	2025-11-01	16.0	0.0	16.0	16.922821	16.0	16.922821	286.381872	1_month
38	drug_1b7500ddafc8	2025-12-01	0.0	8.0	16.0	6.909886	0.0	1.090114	1.188349	1_month
39	drug_1b7500ddafc8	2026-01-01	8.0	16.0	0.0	10.810439	8.0	5.189561	26.931542	1_month
40	drug_1f2ebad2b6dd	2025-06-01	7.0	7.0	14.0	8.987826	7.0	1.987826	3.951454	1_month
41	drug_1f2ebad2b6dd	2025-07-01	7.0	14.0	7.0	13.426182	7.0	0.573818	0.329267	1_month
42	drug_1f2ebad2b6dd	2025-08-01	14.0	0.0	7.0	14.898687	14.0	14.898687	221.970885	1_month
43	drug_1f2ebad2b6dd	2025-09-01	0.0	14.0	14.0	3.088803	0.0	10.911197	119.054224	1_month
44	drug_1f2ebad2b6dd	2025-10-01	14.0	7.0	0.0	9.019532	14.0	2.019532	4.078510	1_month
45	drug_1f2ebad2b6dd	2025-11-01	7.0	0.0	14.0	6.578100	7.0	6.578100	43.271402	1_month
46	drug_1f2ebad2b6dd	2025-12-01	0.0	12.0	7.0	9.641685	0.0	2.358315	5.561652	1_month
47	drug_1f2ebad2b6dd	2026-01-01	12.0	7.0	0.0	6.419474	12.0	0.580526	0.337011	1_month
48	drug_25654b50582c	2025-06-01	6.0	0.0	0.0	5.087248	6.0	5.087248	25.880091	1_month
49	drug_25654b50582c	2025-07-01	0.0	6.0	6.0	4.126983	0.0	1.873017	3.508194	1_month

13. Actual vs predicted visualisations

# Plot actual vs predicted values for each forecast horizon

for horizon_name in targets.keys():

    results_h = results_by_horizon[horizon_name]

    plt.figure(figsize=(10, 6))
    plt.scatter(
        results_h["future_actual_usage"],
        results_h["xgboost_prediction"],
        alpha=0.6
    )

    plt.xlabel("Actual Units")
    plt.ylabel("Predicted Units")
    plt.title(f"XGBoost Actual vs Predicted Monthly Units: {horizon_name} Forecast")

    plt.show()

14. Feature importance and SHAP explainability

# Feature importance for each forecast horizon

importance_by_horizon = {}

for horizon_name in targets.keys():

    model = best_models[horizon_name]

    importance_h = pd.DataFrame({
        "feature": features,
        "importance": model.feature_importances_
    }).sort_values("importance", ascending=False)

    importance_h["horizon"] = horizon_name

    importance_by_horizon[horizon_name] = importance_h

    print(f"\nFeature importance: {horizon_name}")
    display(importance_h)

Feature importance: 1_month

	feature	importance	horizon
11	roll_mean_3	0.359544	1_month
17	roll_max_3	0.201819	1_month
18	roll_max_6	0.116647	1_month
12	roll_mean_6	0.101061	1_month
16	roll_min_6	0.039074	1_month
15	roll_min_3	0.035824	1_month
2	lag_3	0.020443	1_month
34	homecare_share	0.019146	1_month
0	lag_1	0.011246	1_month
40	year	0.010332	1_month
3	lag_4	0.009575	1_month
20	lag_1_to_mean_6_ratio	0.008506	1_month
35	hospital_share	0.007411	1_month
13	roll_std_3	0.006746	1_month
5	lag_6	0.006346	1_month
1	lag_2	0.005133	1_month
39	month_number	0.003394	1_month
14	roll_std_6	0.003380	1_month
37	change_3m	0.003364	1_month
22	lag_1_to_roll_mean_3_ratio	0.003094	1_month
26	growth_1m	0.003025	1_month
41	drug_code	0.002918	1_month
38	change_6m	0.002407	1_month
4	lag_5	0.002317	1_month
31	above_6m_average_flag	0.002300	1_month
30	usage_acceleration	0.002177	1_month
6	lag_12	0.001777	1_month
29	growth_3m_pct	0.001726	1_month
36	change_1m	0.001661	1_month
19	recent_growth_3m	0.001546	1_month
9	used_last_6m_count	0.001036	1_month
27	growth_3m	0.000953	1_month
21	lag_1_to_roll_max_6_ratio	0.000842	1_month
28	growth_1m_pct	0.000840	1_month
24	positive_growth_1m_flag	0.000667	1_month
32	above_6m_max_flag	0.000439	1_month
23	recent_peak_flag	0.000363	1_month
7	used_last_month	0.000323	1_month
8	used_last_3m_count	0.000280	1_month
25	positive_growth_3m_flag	0.000131	1_month
33	recent_growth_flag	0.000112	1_month
10	nonzero_rate_6m	0.000074	1_month

Feature importance: 3_month

	feature	importance	horizon
11	roll_mean_3	0.534766	3_month
17	roll_max_3	0.207711	3_month
15	roll_min_3	0.075591	3_month
12	roll_mean_6	0.060015	3_month
18	roll_max_6	0.055770	3_month
6	lag_12	0.010465	3_month
1	lag_2	0.008763	3_month
40	year	0.004625	3_month
24	positive_growth_1m_flag	0.003715	3_month
16	roll_min_6	0.003681	3_month
0	lag_1	0.003043	3_month
3	lag_4	0.002828	3_month
41	drug_code	0.002674	3_month
34	homecare_share	0.002557	3_month
39	month_number	0.002556	3_month
38	change_6m	0.002172	3_month
5	lag_6	0.001763	3_month
2	lag_3	0.001643	3_month
37	change_3m	0.001553	3_month
13	roll_std_3	0.001467	3_month
27	growth_3m	0.001216	3_month
8	used_last_3m_count	0.001096	3_month
19	recent_growth_3m	0.001065	3_month
14	roll_std_6	0.000979	3_month
30	usage_acceleration	0.000893	3_month
20	lag_1_to_mean_6_ratio	0.000891	3_month
36	change_1m	0.000810	3_month
21	lag_1_to_roll_max_6_ratio	0.000741	3_month
9	used_last_6m_count	0.000723	3_month
10	nonzero_rate_6m	0.000619	3_month
25	positive_growth_3m_flag	0.000574	3_month
4	lag_5	0.000566	3_month
23	recent_peak_flag	0.000491	3_month
26	growth_1m	0.000449	3_month
35	hospital_share	0.000288	3_month
7	used_last_month	0.000268	3_month
28	growth_1m_pct	0.000242	3_month
22	lag_1_to_roll_mean_3_ratio	0.000227	3_month
32	above_6m_max_flag	0.000206	3_month
29	growth_3m_pct	0.000164	3_month
33	recent_growth_flag	0.000134	3_month
31	above_6m_average_flag	0.000000	3_month

Feature importance: 6_month

	feature	importance	horizon
11	roll_mean_3	0.500459	6_month
17	roll_max_3	0.212025	6_month
12	roll_mean_6	0.087932	6_month
15	roll_min_3	0.085026	6_month
18	roll_max_6	0.051429	6_month
1	lag_2	0.008159	6_month
16	roll_min_6	0.007762	6_month
6	lag_12	0.007096	6_month
40	year	0.005422	6_month
4	lag_5	0.004748	6_month
0	lag_1	0.004153	6_month
14	roll_std_6	0.003000	6_month
41	drug_code	0.002956	6_month
3	lag_4	0.002315	6_month
34	homecare_share	0.001945	6_month
38	change_6m	0.001855	6_month
37	change_3m	0.001608	6_month
7	used_last_month	0.001194	6_month
2	lag_3	0.001137	6_month
39	month_number	0.001129	6_month
13	roll_std_3	0.001051	6_month
36	change_1m	0.000796	6_month
5	lag_6	0.000777	6_month
27	growth_3m	0.000748	6_month
30	usage_acceleration	0.000740	6_month
19	recent_growth_3m	0.000542	6_month
26	growth_1m	0.000511	6_month
10	nonzero_rate_6m	0.000445	6_month
23	recent_peak_flag	0.000442	6_month
9	used_last_6m_count	0.000426	6_month
8	used_last_3m_count	0.000376	6_month
22	lag_1_to_roll_mean_3_ratio	0.000333	6_month
25	positive_growth_3m_flag	0.000316	6_month
20	lag_1_to_mean_6_ratio	0.000308	6_month
21	lag_1_to_roll_max_6_ratio	0.000173	6_month
28	growth_1m_pct	0.000146	6_month
35	hospital_share	0.000138	6_month
29	growth_3m_pct	0.000131	6_month
32	above_6m_max_flag	0.000080	6_month
24	positive_growth_1m_flag	0.000080	6_month
31	above_6m_average_flag	0.000070	6_month
33	recent_growth_flag	0.000021	6_month

# Plot feature importance for each forecast horizon

for horizon_name in targets.keys():

    importance_h = importance_by_horizon[horizon_name]

    plt.figure(figsize=(10, 6))
    sns.barplot(
        data=importance_h,
        x="importance",
        y="feature"
    )

    plt.title(f"XGBoost Feature Importance: {horizon_name} Forecast")
    plt.xlabel("Importance")
    plt.ylabel("Feature")
    plt.show()

# SHAP analysis for all forecast horizons

shap_results = {}

for horizon_name in targets.keys():

    model = best_models[horizon_name]
    X_valid_h = horizon_data[horizon_name]["X_valid"]

    explainer = shap.Explainer(model)
    shap_values = explainer(X_valid_h)

    # Retain the results for later local explanations
    shap_results[horizon_name] = {
        "explainer": explainer,
        "shap_values": shap_values,
        "X_valid": X_valid_h
     }

    plt.title(f"SHAP Summary Plot: {horizon_name} Forecast")

    shap.summary_plot(
        shap_values,
        X_valid_h,
        show=True
    )

# Plot a local SHAP waterfall explanation for one drug and forecast origin

def plot_shap_waterfall(
    horizon_name,
    selected_drug,
    selected_month,
    max_display=15
):

    if horizon_name not in shap_results:
        raise ValueError(
            f"Unknown horizon '{horizon_name}'. "
            f"Choose from: {list(shap_results.keys())}"
        )

    # These lines must be outside the if block
    result = shap_results[horizon_name]

    # SHAP objects retained in shap_results
    X_valid_h = (
        result["X_valid"]
        .reset_index(drop=True)
    )

    shap_values_h = result["shap_values"]

    # Validation metadata retained in horizon_data
    valid_h = (
        horizon_data[horizon_name]["valid"]
        .reset_index(drop=True)
    )

    target_col = horizon_data[horizon_name]["target_col"]

    # Convert dates to monthly periods for reliable comparison
    requested_month = (
        pd.Timestamp(selected_month)
        .to_period("M")
    )

    validation_months = (
        pd.to_datetime(valid_h["month_period"])
        .dt.to_period("M")
    )

    selection_mask = (
        valid_h["drugfull"]
        .astype(str)
        .str.strip()
        .str.casefold()
        .eq(str(selected_drug).strip().casefold())
        &
        validation_months.eq(requested_month)
    )

    matching_positions = np.flatnonzero(
        selection_mask.to_numpy()
    )

    if len(matching_positions) == 0:

        available_months = (
            valid_h.loc[
                valid_h["drugfull"]
                .astype(str)
                .str.strip()
                .str.casefold()
                .eq(str(selected_drug).strip().casefold()),
                "month_period"
            ]
            .dt.strftime("%Y-%m")
            .unique()
        )

        raise ValueError(
            f"No validation row found for '{selected_drug}' "
            f"in {requested_month}.\n"
            f"Available validation months for this drug: "
            f"{list(available_months)}"
        )

    if len(matching_positions) > 1:
        print(
            f"Warning: {len(matching_positions)} matching rows found. "
            "The first row will be displayed."
        )

    row_position = int(matching_positions[0])

    # Obtain the model prediction and actual target
    model = best_models[horizon_name]

    model_prediction = float(
        model.predict(
            X_valid_h.iloc[[row_position]]
        )[0]
    )

    actual_usage = float(
        valid_h[target_col].iloc[row_position]
    )

    current_usage = float(
        valid_h["total_units"].iloc[row_position]
    )

    # Determine target period
    horizon_months = {
        "1_month": 1,
        "3_month": 3,
        "6_month": 6
    }

    number_of_months = horizon_months[horizon_name]

    forecast_start = requested_month + 1
    forecast_end = requested_month + number_of_months

    if number_of_months == 1:
        forecast_period_label = forecast_start.strftime(
            "%B %Y"
        )
    else:
        forecast_period_label = (
            f"{forecast_start.strftime('%B %Y')}–"
            f"{forecast_end.strftime('%B %Y')}"
        )

    # Local SHAP waterfall chart
    shap.plots.waterfall(
        shap_values_h[row_position],
        max_display=max_display,
        show=False
    )

    plt.title(
        f"SHAP Waterfall Explanation\n"
        f"{selected_drug}\n"
        f"Forecast origin: "
        f"{requested_month.strftime('%B %Y')} | "
        f"Target period: {forecast_period_label}\n"
        f"Predicted: {model_prediction:,.1f} units | "
        f"Actual: {actual_usage:,.1f} units",
        pad=20
    )

    plt.tight_layout()
    plt.show()

    explanation = shap_values_h[row_position]

    base_value = float(
        np.asarray(
            explanation.base_values
        ).reshape(-1)[0]
    )

    reconstructed_prediction = (
        base_value
        + explanation.values.sum()
    )

    result_summary = {
        "horizon": horizon_name,
        "drug": selected_drug,
        "forecast_origin": str(requested_month),
        "forecast_period": forecast_period_label,
        "current_month_usage": current_usage,
        "model_prediction": model_prediction,
        "actual_usage": actual_usage,
        "base_value": base_value,
        "total_shap_contribution": float(
            explanation.values.sum()
        ),
        "reconstructed_prediction": float(
            reconstructed_prediction
        ),
        "validation_row_position": row_position
    }

    display(
        pd.DataFrame(
            [result_summary]
        ).T.rename(columns={0: "Value"})
    )

    return result_summary

waterfall_result = plot_shap_waterfall(
    horizon_name="3_month",
    selected_drug="drug_2b4752a1e609",
    selected_month="2025-11",
    max_display=15
)

	Value
horizon	3_month
drug	drug_2b4752a1e609
forecast_origin	2025-11
forecast_period	December 2025–February 2026
current_month_usage	504.0
model_prediction	1291.934082
actual_usage	1479.0
base_value	271.367279
total_shap_contribution	1020.567505
reconstructed_prediction	1291.934814
validation_row_position	62

SHAP explainability

SHAP is used to interpret how the fitted XGBoost models use the engineered features when producing forecasts. The summary plots help identify whether recent usage, rolling statistics, growth indicators, homecare/hospital split, and calendar features are contributing meaningfully to predictions.

15. Diebold-Mariano tests

# Function to run Diebold-Mariano test between the forecast horizons

def run_dm_test(results_a, results_b, h):

    merge_cols = ["drugfull", "month_period"]

    aligned = results_a.merge(
        results_b,
        on=merge_cols,
        suffixes=("_a", "_b")
    )

    actual = aligned["future_actual_usage_a"].values
    forecast_a = aligned["xgboost_prediction_a"].values
    forecast_b = aligned["xgboost_prediction_b"].values

    stat, p_value = dm_test(
        actual,
        forecast_a,
        forecast_b,
        h=h,
        loss=lambda actual, predicted: (actual - predicted) ** 2,
        one_sided=False
    )

    return {
        "DM Statistic": stat,
        "p-value": p_value,
        "Compared Rows": len(aligned)
    }

# Run Diebold-Mariano comparisons for the horizons

dm_results = {
    "1m_vs_3m": run_dm_test(
        results_by_horizon["1_month"],
        results_by_horizon["3_month"],
        h=3
    ),
    "1m_vs_6m": run_dm_test(
        results_by_horizon["1_month"],
        results_by_horizon["6_month"],
        h=6
    ),
    "3m_vs_6m": run_dm_test(
        results_by_horizon["3_month"],
        results_by_horizon["6_month"],
        h=6
    )

}

dm_results_df = pd.DataFrame(dm_results).T

display(dm_results_df)

	DM Statistic	p-value	Compared Rows
1m_vs_3m	-1.619193	0.106419	312.0
1m_vs_6m	-1.210738	0.227839	156.0
3m_vs_6m	-1.142006	0.254771	208.0

# Diebold-Mariano test: XGBoost vs naive baseline for each horizon

def run_dm_test_against_naive(results_h, h):

    actual = results_h["future_actual_usage"].values
    xgb_forecast = results_h["xgboost_prediction"].values
    naive_forecast = results_h["naive_forecast"].values

    stat, p_value = dm_test(
        actual,
        xgb_forecast,
        naive_forecast,
        h=h,
        loss=lambda actual, predicted: (actual - predicted) ** 2,
        one_sided=False
    )

    return {
        "DM Statistic": stat,
        "p-value": p_value,
        "Compared Rows": len(results_h)
    }

# Run Diebold-Mariano tests against naive baseline

dm_naive_results = {
    "1m_xgboost_vs_naive": run_dm_test_against_naive(
        results_by_horizon["1_month"],
        h=1
    ),
    "3m_xgboost_vs_naive": run_dm_test_against_naive(
        results_by_horizon["3_month"],
        h=3
    ),
    "6m_xgboost_vs_naive": run_dm_test_against_naive(
        results_by_horizon["6_month"],
        h=6
    )
}

dm_naive_results_df = pd.DataFrame(dm_naive_results).T

display(dm_naive_results_df)

	DM Statistic	p-value	Compared Rows
1m_xgboost_vs_naive	-1.095287	0.274026	416.0
3m_xgboost_vs_naive	0.499815	0.617508	364.0
6m_xgboost_vs_naive	-1.105235	0.269790	364.0

16. Export main forecast results

all_forecast_results.to_csv("xgboost_forecast_results.csv", index=False)

17. Drug-level, monthly, and high-volume residual diagnostics

drug_error = (
    all_forecast_results
    .groupby(["horizon", "drugfull"], as_index=False)
    .agg(
        actual_total=("future_actual_usage", "sum"),
        predicted_total=("xgboost_prediction", "sum"),
        mae=("absolute_error", "mean"),
        rmse=("squared_error", lambda x: np.sqrt(np.mean(x))),
        rows=("future_actual_usage", "size")
    )
)

drug_error["absolute_total_error"] = (
    drug_error["actual_total"] - drug_error["predicted_total"]
).abs()

display(
    drug_error
    .sort_values("absolute_total_error", ascending=False)
    .head(25)
)

	horizon	drugfull	actual_total	predicted_total	mae	rmse	rows	absolute_total_error
112	6_month	drug_2b4752a1e609	23934.0	21295.802734	376.885393	396.784973	7	2638.197266
139	6_month	drug_9e65f61d9995	6364.0	4484.419922	268.511475	291.086574	7	1879.580078
60	3_month	drug_2b4752a1e609	11842.0	10081.479492	254.625959	303.516003	7	1760.520508
140	6_month	drug_a3ffd5cb15c9	5264.0	3690.171387	224.832663	228.362405	7	1573.828613
87	3_month	drug_9e65f61d9995	3488.0	2602.306396	126.527653	133.294231	7	885.693604
116	6_month	drug_3eb69a66366e	2360.0	1563.537842	113.780299	118.078846	7	796.462158
143	6_month	drug_b3fdd633217f	2432.0	1767.658325	94.905960	101.155260	7	664.341675
108	6_month	drug_1b7500ddafc8	628.0	1234.799438	103.273078	137.371361	7	606.799438
148	6_month	drug_d033588a88cb	2128.0	1524.357300	86.234674	99.521943	7	603.642700
8	1_month	drug_2b4752a1e609	4349.0	3762.479492	146.191967	171.174236	8	586.520508
88	3_month	drug_a3ffd5cb15c9	2797.0	2289.917480	72.440369	80.336142	7	507.082520
111	6_month	drug_29f00a08b8da	1523.0	1121.322144	86.317089	94.775766	7	401.677856
128	6_month	drug_6cc0e0edce7d	411.0	43.275055	53.772938	73.382107	7	367.724945
113	6_month	drug_2fa7ea18cba1	1475.0	1176.429077	56.967905	68.694825	7	298.570923
64	3_month	drug_3eb69a66366e	1310.0	1024.887451	40.730366	46.678960	7	285.112549
122	6_month	drug_57f8fe0a5f8a	1528.0	1245.492310	52.986300	65.024093	7	282.507690
123	6_month	drug_58e790c405a1	4462.0	4180.391602	65.258702	75.121236	7	281.608398
154	6_month	drug_edea72b252ae	370.0	125.252769	34.963890	38.094367	7	244.747231
104	6_month	drug_0050f02cadea	3975.0	3740.293213	56.669695	68.591568	7	234.706787
35	1_month	drug_9e65f61d9995	1348.0	1117.303467	58.621080	72.302505	8	230.696533
76	3_month	drug_6cc0e0edce7d	335.0	118.115036	31.850308	40.012004	7	216.884964
131	6_month	drug_7eae036aa176	9.0	213.996277	29.285183	34.729068	7	204.996277
130	6_month	drug_7cab3c1e3f92	660.0	494.698090	23.614559	26.993348	7	165.301910
117	6_month	drug_3edead72ce43	720.0	555.340637	44.130114	46.844532	7	164.659363
97	3_month	drug_d174b29edcbb	546.0	697.972717	28.078498	29.937360	7	151.972717

plt.figure(figsize=(8, 6))
plt.scatter(all_forecast_results["future_actual_usage"], all_forecast_results["xgboost_prediction"], alpha=0.5)
plt.xlabel("Actual usage")
plt.ylabel("Predicted usage")
plt.title("Actual vs Predicted Usage")

# Add perfect prediction line
max_val = max(all_forecast_results["future_actual_usage"].max(), all_forecast_results["xgboost_prediction"].max())
plt.plot([0, max_val], [0, max_val], linestyle="--")

plt.show()

all_forecast_results["residual"] = all_forecast_results["future_actual_usage"] - all_forecast_results["xgboost_prediction"]

high_volume = all_forecast_results[all_forecast_results["future_actual_usage"] >= all_forecast_results["future_actual_usage"].quantile(0.95)]

high_volume_summary = high_volume.agg(
    actual_total=("future_actual_usage", "sum"),
    predicted_total=("xgboost_prediction", "sum"),
    mean_residual=("residual", "mean"),
    total_residual=("residual", "sum"),
    mean_absolute_error=("absolute_error", "mean")
)

high_volume_summary

	future_actual_usage	xgboost_prediction	residual	absolute_error
actual_total	64311.0	NaN	NaN	NaN
predicted_total	NaN	54004.949219	NaN	NaN
mean_residual	NaN	NaN	177.690510	NaN
total_residual	NaN	NaN	10306.049591	NaN
mean_absolute_error	NaN	NaN	NaN	183.880944

monthly_error = (
    all_forecast_results
    .groupby("month_period")
    .agg(
        actual_total=("future_actual_usage", "sum"),
        predicted_total=("xgboost_prediction", "sum"),
        total_absolute_error=("absolute_error", "sum"),
        mean_absolute_error=("absolute_error", "mean"),
        residual_total=("residual", "sum")
    )
    .reset_index()
)

monthly_error["wmape"] = (
    monthly_error["total_absolute_error"] / monthly_error["actual_total"]
) * 100

monthly_error["bias_percent"] = (
    monthly_error["residual_total"] / monthly_error["actual_total"]
) * 100

monthly_error.sort_values("residual_total", ascending=False).head(10)

	month_period	actual_total	predicted_total	total_absolute_error	mean_absolute_error	residual_total	wmape	bias_percent
5	2025-07-01	15977.0	13279.587891	3824.995919	24.519205	2697.411812	23.940639	16.883093
4	2025-06-01	15521.0	13887.346680	3471.614327	22.253938	1633.653401	22.367208	10.525439
3	2025-05-01	13501.0	11902.104492	2992.671320	28.775686	1598.895580	22.166294	11.842794
7	2025-09-01	6676.0	5234.981934	2059.819594	19.805958	1441.017991	30.854098	21.585051
1	2025-03-01	8731.0	7306.142578	2205.801053	42.419251	1424.857323	25.264014	16.319520
6	2025-08-01	15134.0	13712.192383	3603.326337	23.098246	1421.807493	23.809478	9.394790
2	2025-04-01	9166.0	7783.924805	2536.449369	48.777872	1382.075148	27.672369	15.078280
0	2025-02-01	8193.0	7461.247559	1849.636615	35.569935	731.752589	22.575816	8.931436
8	2025-10-01	6109.0	5621.580078	1680.541461	16.159053	487.419934	27.509273	7.978719
10	2025-12-01	1973.0	1516.913452	652.209100	12.542483	456.086582	33.056721	23.116400

drug_error = (
    all_forecast_results
    .groupby("drugfull")
    .agg(
        actual_total=("future_actual_usage", "sum"),
        predicted_total=("xgboost_prediction", "sum"),
        total_absolute_error=("absolute_error", "sum"),
        residual_total=("residual", "sum"),
        count=("future_actual_usage", "count")
    )
    .reset_index()
)

drug_error["wmape"] = (
    drug_error["total_absolute_error"] / drug_error["actual_total"]
) * 100

drug_error["bias_percent"] = (
    drug_error["residual_total"] / drug_error["actual_total"]
) * 100

drug_error.sort_values("total_absolute_error", ascending=False).head(20)

	drugfull	actual_total	predicted_total	total_absolute_error	residual_total	count	wmape	bias_percent
8	drug_2b4752a1e609	40125.0	35139.761719	5590.115204	4985.238495	22	13.931751	12.424270
35	drug_9e65f61d9995	11200.0	8204.029297	3234.242538	2995.970490	22	28.877166	26.749737
36	drug_a3ffd5cb15c9	9083.0	7004.543945	2499.540840	2078.456245	22	27.518891	22.882927
12	drug_3eb69a66366e	4174.0	3032.086182	1230.029751	1141.913891	22	29.468849	27.357784
7	drug_29f00a08b8da	2562.0	2299.563477	1045.829578	262.436588	22	40.820827	10.243427
19	drug_58e790c405a1	7498.0	7049.458984	1022.614098	448.540764	22	13.638492	5.982139
4	drug_1b7500ddafc8	1032.0	1696.499390	1014.207725	-664.499403	22	98.275942	-64.389477
44	drug_d033588a88cb	3618.0	2980.554932	976.724354	637.445141	22	26.996251	17.618716
0	drug_0050f02cadea	6773.0	6359.144043	882.003128	413.856003	22	13.022341	6.110379
39	drug_b3fdd633217f	4168.0	3522.036377	829.066372	645.963627	22	19.891228	15.498168
9	drug_2fa7ea18cba1	2487.0	2303.697754	788.003723	183.302231	22	31.684910	7.370415
13	drug_3edead72ce43	1234.0	1230.294678	763.819198	3.705279	22	61.897828	0.300266
24	drug_6cc0e0edce7d	880.0	251.531891	677.052040	628.468112	22	76.937732	71.416831
18	drug_57f8fe0a5f8a	2588.0	2401.154541	562.198442	186.845377	22	21.723278	7.219682
50	drug_edea72b252ae	624.0	356.279114	438.809605	267.720867	22	70.322052	42.903985
45	drug_d174b29edcbb	1865.0	2066.634521	416.763077	-201.634548	22	22.346546	-10.811504
32	drug_91d77d9a9e30	382.0	580.453369	409.363501	-198.453351	22	107.163220	-51.951139
26	drug_7cab3c1e3f92	1136.0	921.774170	294.445118	214.225828	22	25.919465	18.857907
29	drug_86055c3c4933	401.0	237.198715	265.095566	163.801283	22	66.108620	40.848200
5	drug_1f2ebad2b6dd	562.0	726.501099	257.810661	-164.501130	22	45.873783	-29.270664

for horizon_name, results_h in results_by_horizon.items():

    plt.figure(figsize=(8, 6))

    plt.scatter(
        results_h["future_actual_usage"],
        results_h["xgboost_prediction"],
        alpha=0.6
    )

    max_val = max(
        results_h["future_actual_usage"].max(),
        results_h["xgboost_prediction"].max()
    )

    plt.plot([0, max_val], [0, max_val], linestyle="--")

    plt.xlabel("Actual usage")
    plt.ylabel("Predicted usage")
    plt.title(f"Actual vs Predicted Usage: {horizon_name}")

    plt.show()

horizon_bias = {}

for horizon_name, results_h in results_by_horizon.items():

    results_h = results_h.copy()

    results_h["residual"] = (
        results_h["future_actual_usage"]
        - results_h["xgboost_prediction"]
    )

    high_volume = results_h[
        results_h["future_actual_usage"]
        >= results_h["future_actual_usage"].quantile(0.95)
    ]

    actual_total = high_volume["future_actual_usage"].sum()
    predicted_total = high_volume["xgboost_prediction"].sum()
    residual_total = high_volume["residual"].sum()

    horizon_bias[horizon_name] = {
        "actual_total": actual_total,
        "predicted_total": predicted_total,
        "residual_total": residual_total,
        "bias_percent": (residual_total / actual_total) * 100
    }

horizon_bias_df = pd.DataFrame(horizon_bias).T
display(horizon_bias_df)

	actual_total	predicted_total	residual_total	bias_percent
1_month	6586.0	5364.646484	1221.353157	18.544688
3_month	17435.0	14340.492188	3094.506470	17.748818
6_month	34272.0	28587.982422	5684.017456	16.585018

18. Diagnosing conservative high-volume predictions

train_valid_range_check = {}

for horizon_name, data in horizon_data.items():

    train = data["train"]
    valid = data["valid"]
    target_col = data["target_col"]

    train_valid_range_check[horizon_name] = {
        "train_target_max": train[target_col].max(),
        "valid_target_max": valid[target_col].max(),
        "train_target_95th": train[target_col].quantile(0.95),
        "valid_target_95th": valid[target_col].quantile(0.95),
        "train_target_mean": train[target_col].mean(),
        "valid_target_mean": valid[target_col].mean()
    }

range_check_df = pd.DataFrame(train_valid_range_check).T
display(range_check_df)

	train_target_max	valid_target_max	train_target_95th	valid_target_95th	train_target_mean	valid_target_mean
1_month	736.0	702.0	88.00	122.75	22.440559	28.848558
3_month	1850.0	1780.0	269.85	395.70	68.561899	88.903846
6_month	3604.0	3529.0	511.65	692.80	137.541777	173.467033

prediction_range_check = {}

for horizon_name, results_h in results_by_horizon.items():

    data = horizon_data[horizon_name]
    train = data["train"]
    target_col = data["target_col"]

    prediction_range_check[horizon_name] = {
        "train_target_max": train[target_col].max(),
        "valid_actual_max": results_h["future_actual_usage"].max(),
        "prediction_max": results_h["xgboost_prediction"].max(),
        "valid_actual_95th": results_h["future_actual_usage"].quantile(0.95),
        "prediction_95th": results_h["xgboost_prediction"].quantile(0.95)
    }

prediction_range_df = pd.DataFrame(prediction_range_check).T
display(prediction_range_df)

	train_target_max	valid_actual_max	prediction_max	valid_actual_95th	prediction_95th
1_month	736.0	702.0	598.765747	122.75	108.368975
3_month	1850.0	1780.0	1719.930786	395.70	317.343959
6_month	3604.0	3529.0	3178.005127	692.80	579.210150

drug_bias_by_horizon = []

for horizon_name, results_h in results_by_horizon.items():

    temp = results_h.copy()

    temp["residual"] = (
        temp["future_actual_usage"]
        - temp["xgboost_prediction"]
    )

    drug_bias = (
        temp
        .groupby("drugfull", as_index=False)
        .agg(
            actual_total=("future_actual_usage", "sum"),
            predicted_total=("xgboost_prediction", "sum"),
            residual_total=("residual", "sum"),
            absolute_error_total=("absolute_error", "sum"),
            rows=("future_actual_usage", "size")
        )
    )

    drug_bias["bias_percent"] = (
        drug_bias["residual_total"] / drug_bias["actual_total"]
    ) * 100

    drug_bias["horizon"] = horizon_name

    drug_bias_by_horizon.append(drug_bias)

drug_bias_by_horizon = pd.concat(
    drug_bias_by_horizon,
    ignore_index=True
)

display(
    drug_bias_by_horizon
    .sort_values("residual_total", ascending=False)
    
)

	drugfull	actual_total	predicted_total	residual_total	absolute_error_total	rows	bias_percent	horizon
112	drug_2b4752a1e609	23934.0	21295.802734	2638.197754	2638.197754	7	11.022803	6_month
139	drug_9e65f61d9995	6364.0	4484.419922	1879.580322	1879.580322	7	29.534575	6_month
60	drug_2b4752a1e609	11842.0	10081.479492	1760.520142	1782.381714	7	14.866747	3_month
140	drug_a3ffd5cb15c9	5264.0	3690.171387	1573.828644	1573.828644	7	29.897961	6_month
87	drug_9e65f61d9995	3488.0	2602.306396	885.693573	885.693573	7	25.392591	3_month
...	...	...	...	...	...	...	...	...
138	drug_9a1aa198dda4	4.0	108.684013	-104.684013	104.684013	7	-2617.100334	6_month
109	drug_1f2ebad2b6dd	344.0	457.388214	-113.388218	165.430668	7	-32.961691	6_month
97	drug_d174b29edcbb	546.0	697.972717	-151.972729	196.549488	7	-27.833833	3_month
131	drug_7eae036aa176	9.0	213.996277	-204.996283	204.996283	7	-2277.736473	6_month
108	drug_1b7500ddafc8	628.0	1234.799438	-606.799438	722.911545	7	-96.624114	6_month

156 rows × 8 columns

problem_drugs = [
    "drug_2b4752a1e609",
    "drug_9e65f61d9995",
    "drug_a3ffd5cb15c9",
    "drug_3eb69a66366e",
    "drug_b3fdd633217f",
    "drug_d033588a88cb"
]

problem_drug_shift = []

for horizon_name, data in horizon_data.items():

    train = data["train"]
    valid = data["valid"]
    target_col = data["target_col"]

    for drug in problem_drugs:

        train_drug = train[train["drugfull"] == drug]
        valid_drug = valid[valid["drugfull"] == drug]

        problem_drug_shift.append({
            "horizon": horizon_name,
            "drugfull": drug,
            "train_mean_target": train_drug[target_col].mean(),
            "valid_mean_target": valid_drug[target_col].mean(),
            "train_max_target": train_drug[target_col].max(),
            "valid_max_target": valid_drug[target_col].max(),
            "train_total_target": train_drug[target_col].sum(),
            "valid_total_target": valid_drug[target_col].sum(),
            "train_rows": len(train_drug),
            "valid_rows": len(valid_drug)
        })

problem_drug_shift_df = pd.DataFrame(problem_drug_shift)

problem_drug_shift_df["mean_change_percent"] = (
    (
        problem_drug_shift_df["valid_mean_target"]
        - problem_drug_shift_df["train_mean_target"]
    )
    / problem_drug_shift_df["train_mean_target"]
) * 100

display(problem_drug_shift_df)

	horizon	drugfull	train_mean_target	valid_mean_target	train_max_target	valid_max_target	train_total_target	valid_total_target	train_rows	valid_rows	mean_change_percent
0	1_month	drug_2b4752a1e609	436.454545	543.625000	736.0	702.0	14403.0	4349.0	33	8	24.554780
1	1_month	drug_9e65f61d9995	32.454545	168.500000	108.0	272.0	1071.0	1348.0	33	8	419.187675
2	1_month	drug_a3ffd5cb15c9	59.606061	127.750000	115.0	197.0	1967.0	1022.0	33	8	114.323843
3	1_month	drug_3eb69a66366e	9.272727	63.000000	50.0	86.0	306.0	504.0	33	8	579.411765
4	1_month	drug_b3fdd633217f	44.424242	58.500000	89.0	87.0	1466.0	468.0	33	8	31.684857
5	1_month	drug_d033588a88cb	45.727273	51.375000	88.0	80.0	1509.0	411.0	33	8	12.350895
6	3_month	drug_2b4752a1e609	1329.500000	1691.714286	1850.0	1780.0	42544.0	11842.0	32	7	27.244399
7	3_month	drug_9e65f61d9995	104.281250	498.285714	404.0	552.0	3337.0	3488.0	32	7	377.828674
8	3_month	drug_a3ffd5cb15c9	181.437500	399.571429	303.0	445.0	5806.0	2797.0	32	7	120.225383
9	3_month	drug_3eb69a66366e	29.875000	187.142857	122.0	214.0	956.0	1310.0	32	7	526.419605
10	3_month	drug_b3fdd633217f	134.968750	181.142857	223.0	211.0	4319.0	1268.0	32	7	34.210962
11	3_month	drug_d033588a88cb	141.093750	154.142857	201.0	178.0	4515.0	1079.0	32	7	9.248537
12	6_month	drug_2b4752a1e609	2653.689655	3419.142857	3604.0	3529.0	76957.0	23934.0	29	7	28.844865
13	6_month	drug_9e65f61d9995	196.758621	909.142857	638.0	1044.0	5706.0	6364.0	29	7	362.059987
14	6_month	drug_a3ffd5cb15c9	361.448276	752.000000	578.0	844.0	10482.0	5264.0	29	7	108.051898
15	6_month	drug_3eb69a66366e	54.000000	337.142857	210.0	408.0	1566.0	2360.0	29	7	524.338624
16	6_month	drug_b3fdd633217f	271.620690	347.428571	379.0	389.0	7877.0	2432.0	29	7	27.909465
17	6_month	drug_d033588a88cb	285.379310	304.000000	367.0	317.0	8276.0	2128.0	29	7	6.524891

19. Persistent underprediction and problem-drug analysis

persistent_underpredicted_drugs = (
    drug_bias_by_horizon
    .query("bias_percent > 10")
    .groupby("drugfull")
    .agg(
        horizons_underpredicted=("horizon", "nunique"),
        total_actual=("actual_total", "sum"),
        total_predicted=("predicted_total", "sum"),
        total_residual=("residual_total", "sum")
    )
    .reset_index()
)

persistent_underpredicted_drugs["overall_bias_percent"] = (
    persistent_underpredicted_drugs["total_residual"]
    / persistent_underpredicted_drugs["total_actual"]
) * 100

persistent_underpredicted_drugs = (
    persistent_underpredicted_drugs
    .sort_values("total_residual", ascending=False)
)

display(persistent_underpredicted_drugs)

	drugfull	horizons_underpredicted	total_actual	total_predicted	total_residual	overall_bias_percent
3	drug_2b4752a1e609	3	40125.0	35139.761719	4985.238495	12.424270
14	drug_9e65f61d9995	3	11200.0	8204.029297	2995.970490	26.749737
15	drug_a3ffd5cb15c9	2	8061.0	5980.088867	2080.911224	25.814554
5	drug_3eb69a66366e	3	4174.0	3032.085938	1141.913891	27.357784
18	drug_b3fdd633217f	1	2432.0	1767.658325	664.341721	27.316683
9	drug_6cc0e0edce7d	3	880.0	251.531891	628.468112	71.416831
19	drug_d033588a88cb	1	2128.0	1524.357300	603.642715	28.366669
2	drug_29f00a08b8da	1	1523.0	1121.322144	401.677872	26.374122
4	drug_2fa7ea18cba1	1	1475.0	1176.429077	298.570869	20.242093
8	drug_57f8fe0a5f8a	1	1528.0	1245.492310	282.507736	18.488726
24	drug_edea72b252ae	2	557.0	286.234192	270.765821	48.611458
11	drug_7cab3c1e3f92	2	1006.0	800.963379	205.036652	20.381377
20	drug_d34593a50053	3	312.0	140.981552	171.018446	54.813605
6	drug_3edead72ce43	1	720.0	555.340637	164.659330	22.869351
12	drug_86055c3c4933	3	401.0	237.198715	163.801283	40.848200
0	drug_0050f02cadea	1	760.0	661.761719	98.238297	12.926092
13	drug_959f0efe359e	1	603.0	516.578552	86.421432	14.331913
16	drug_a8f17e018e69	2	116.0	42.025211	73.974790	63.771370
21	drug_d6440d12de4f	3	112.0	52.897522	59.102479	52.770070
25	drug_f80641253af7	3	208.0	150.550674	57.449317	27.619864
23	drug_e2391f09ae75	3	180.0	139.665741	40.334249	22.407916
10	drug_6f84645e7b70	2	138.0	102.340256	35.659736	25.840389
1	drug_25654b50582c	1	102.0	76.971970	25.028033	24.537287
17	drug_af2e33228425	1	38.0	31.508583	6.491417	17.082676
7	drug_46ae92e1fec0	1	10.0	6.338623	3.661378	36.613776
22	drug_def417c600be	2	16.0	13.557957	2.442043	15.262771

February cut-off forecasting demonstration

This section simulates an operational forecasting scenario. The model is trained using data available before February 2025, then February is used as the forecast origin to predict cumulative March–May usage. The analysis compares the predicted three-month usage against actual observed usage at total and drug level.

20. February cut-off practical forecast

# Practical cut-off analysis
# Forecast 3-month usage from February actual usage
# This simulates using data available at the end of February to forecast March-May usage

forecast_cutoff = pd.Timestamp("2025-02-01")

horizon_name = "3_month"
target_col = "target_3m"

# Keep rows where 3-month target exists
cutoff_df = base_df.dropna(
    subset=[target_col]
).reset_index(drop=True)

# Train on data before February
train_cutoff = cutoff_df[
    cutoff_df["month_period"] < forecast_cutoff
].copy()

# Forecast from February
test_cutoff = cutoff_df[
    cutoff_df["month_period"] == forecast_cutoff
].copy()

X_train_cutoff = train_cutoff[features]
y_train_cutoff = train_cutoff[target_col]

X_test_cutoff = test_cutoff[features]
y_test_cutoff = test_cutoff[target_col]

print("Training data up to:", train_cutoff["month_period"].max())
print("Forecast origin:", forecast_cutoff)
print("Forecast period: March-May 2025")
print("Training rows:", len(train_cutoff))
print("Test rows:", len(test_cutoff))

Training data up to: 2025-01-01 00:00:00
Forecast origin: 2025-02-01 00:00:00
Forecast period: March-May 2025
Training rows: 1508
Test rows: 52

# Train cut-off model using best 3-month parameters

best_3m_params = best_params["3_month"]

cutoff_model = XGBRegressor(
    objective="reg:squarederror",
    random_state=42,
    n_jobs=-1,
    **best_3m_params
)

# Use the same tiered weighting approach as the main model
q85 = y_train_cutoff.quantile(0.85)
q90 = y_train_cutoff.quantile(0.90)
q95 = y_train_cutoff.quantile(0.95)

sample_weights_cutoff = np.select(
    [
        y_train_cutoff >= q95,
        y_train_cutoff >= q90,
        y_train_cutoff >= q85
    ],
    [
        12,
        8,
        4
    ],
    default=1
)

cutoff_model.fit(
    X_train_cutoff,
    y_train_cutoff,
    sample_weight=sample_weights_cutoff
)

XGBRegressor(base_score=None, booster=None, callbacks=None,
             colsample_bylevel=None, colsample_bynode=None,
             colsample_bytree=0.9, device=None, early_stopping_rounds=None,
             enable_categorical=False, eval_metric=None, feature_types=None,
             feature_weights=None, gamma=None, grow_policy=None,
             importance_type=None, interaction_constraints=None,
             learning_rate=0.03, max_bin=None, max_cat_threshold=None,
             max_cat_to_onehot=None, max_delta_step=None, max_depth=4,
             max_leaves=None, min_child_weight=3, missing=nan,
             monotone_constraints=None, multi_strategy=None, n_estimators=500,
             n_jobs=-1, num_parallel_tree=None, ...)

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

# Predict March-May cumulative usage from February

test_cutoff["predicted_3m_usage"] = cutoff_model.predict(X_test_cutoff)

test_cutoff["predicted_3m_usage"] = np.clip(
    test_cutoff["predicted_3m_usage"],
    0,
    None
)

test_cutoff["actual_3m_usage"] = y_test_cutoff

test_cutoff["absolute_error"] = np.abs(
    test_cutoff["actual_3m_usage"]
    - test_cutoff["predicted_3m_usage"]
)

test_cutoff["squared_error"] = (
    test_cutoff["actual_3m_usage"]
    - test_cutoff["predicted_3m_usage"]
) ** 2

display(
    test_cutoff[
        [
            "drugfull",
            "month_period",
            "total_units",
            "actual_3m_usage",
            "predicted_3m_usage",
            "absolute_error"
        ]
    ].sort_values("absolute_error", ascending=False)
)

	drugfull	month_period	total_units	actual_3m_usage	predicted_3m_usage	absolute_error
341	drug_2b4752a1e609	2025-02-01	381.0	1568.0	1388.362305	179.637695
185	drug_1b7500ddafc8	2025-02-01	16.0	68.0	174.238800	106.238800
1433	drug_a3ffd5cb15c9	2025-02-01	72.0	289.0	204.864685	84.135315
1394	drug_9e65f61d9995	2025-02-01	62.0	280.0	209.498734	70.501266
770	drug_58e790c405a1	2025-02-01	96.0	260.0	323.167389	63.167389
497	drug_3eb69a66366e	2025-02-01	22.0	122.0	62.211151	59.788849
380	drug_2fa7ea18cba1	2025-02-01	30.0	93.0	50.212574	42.787426
107	drug_119bcb15803e	2025-02-01	4.0	42.0	16.031942	25.968058
1784	drug_d174b29edcbb	2025-02-01	22.0	79.0	102.852882	23.852882
1082	drug_7eae036aa176	2025-02-01	14.0	9.0	32.525101	23.525101
1628	drug_b8556e04d039	2025-02-01	0.0	31.0	11.264912	19.735088
1043	drug_7cab3c1e3f92	2025-02-01	7.0	54.0	36.406353	17.593647
1550	drug_b3fdd633217f	2025-02-01	20.0	143.0	126.440796	16.559204
536	drug_3edead72ce43	2025-02-01	0.0	28.0	13.484136	14.515864
731	drug_57f8fe0a5f8a	2025-02-01	30.0	96.0	82.630074	13.369926
848	drug_5df05892de2b	2025-02-01	6.0	0.0	11.664601	11.664601
1667	drug_c6283c502b77	2025-02-01	0.0	12.0	0.341338	11.658662
224	drug_1f2ebad2b6dd	2025-02-01	0.0	33.0	23.279350	9.720650
2018	drug_f80641253af7	2025-02-01	0.0	10.0	0.349325	9.650675
419	drug_37fd409a7a07	2025-02-01	8.0	30.0	20.754797	9.245203
1277	drug_91d77d9a9e30	2025-02-01	4.0	22.0	12.784729	9.215271
68	drug_0d3c27ea4824	2025-02-01	0.0	0.0	8.599671	8.599671
263	drug_25654b50582c	2025-02-01	0.0	6.0	14.167142	8.167142
809	drug_5be70debf5b5	2025-02-01	1.0	2.0	10.033858	8.033858
1316	drug_959f0efe359e	2025-02-01	4.0	43.0	37.395466	5.604534
1940	drug_e2391f09ae75	2025-02-01	0.0	5.0	0.000000	5.000000
575	drug_3fb3e5685f75	2025-02-01	0.0	0.0	4.968575	4.968575
1160	drug_86055c3c4933	2025-02-01	0.0	5.0	0.401656	4.598344
29	drug_0050f02cadea	2025-02-01	80.0	244.0	248.556046	4.556046
302	drug_29f00a08b8da	2025-02-01	26.0	76.0	71.608765	4.391235
1511	drug_af2e33228425	2025-02-01	0.0	4.0	8.078116	4.078116
146	drug_139f6537068d	2025-02-01	0.0	0.0	2.627234	2.627234
653	drug_4b3d76e7691e	2025-02-01	0.0	0.0	2.281160	2.281160
1121	drug_825fe42cacc9	2025-02-01	0.0	0.0	2.175764	2.175764
1979	drug_edea72b252ae	2025-02-01	3.0	4.0	2.328003	1.671997
692	drug_4c8fbeacf83d	2025-02-01	0.0	0.0	1.593337	1.593337
1472	drug_a8f17e018e69	2025-02-01	0.0	0.0	1.562813	1.562813
1355	drug_9a1aa198dda4	2025-02-01	0.0	4.0	2.574070	1.425930
458	drug_3c3d57da1b08	2025-02-01	0.0	0.0	0.969707	0.969707
1745	drug_d033588a88cb	2025-02-01	25.0	121.0	120.360420	0.639580
614	drug_46ae92e1fec0	2025-02-01	0.0	3.0	2.568740	0.431260
887	drug_5e71b18ec9d4	2025-02-01	0.0	0.0	0.401656	0.401656
926	drug_6c5601ca5cde	2025-02-01	0.0	0.0	0.401656	0.401656
965	drug_6cc0e0edce7d	2025-02-01	0.0	0.0	0.401656	0.401656
1004	drug_6f84645e7b70	2025-02-01	0.0	0.0	0.401656	0.401656
1199	drug_8b8863408f28	2025-02-01	0.0	0.0	0.401656	0.401656
1589	drug_b6b85ff6a7c2	2025-02-01	0.0	0.0	0.392804	0.392804
1706	drug_c75065c2f034	2025-02-01	0.0	0.0	0.341338	0.341338
1238	drug_907adc10f7c8	2025-02-01	0.0	2.0	1.923756	0.076244
1823	drug_d34593a50053	2025-02-01	0.0	0.0	0.000000	0.000000
1862	drug_d6440d12de4f	2025-02-01	0.0	0.0	0.000000	0.000000
1901	drug_def417c600be	2025-02-01	0.0	0.0	0.000000	0.000000

# Plot total actual vs predicted March-May usage

total_actual_3m = test_cutoff["actual_3m_usage"].sum()
total_predicted_3m = test_cutoff["predicted_3m_usage"].sum()

plt.figure(figsize=(6, 5))

plt.bar(
    ["Actual March-May", "Predicted March-May"],
    [total_actual_3m, total_predicted_3m]
)

plt.ylabel("Total units")
plt.title("3-Month Biologics Usage Forecast from February Cut-Off")

plt.show()

# Drug-level February cut-off forecast comparison

drug_cutoff_results = test_cutoff[
    [
        "drugfull",
        "total_units",
        "actual_3m_usage",
        "predicted_3m_usage",
        "absolute_error"
    ]
].copy()

drug_cutoff_results["bias"] = (
    drug_cutoff_results["actual_3m_usage"]
    - drug_cutoff_results["predicted_3m_usage"]
)

drug_cutoff_results["bias_percent"] = (
    drug_cutoff_results["bias"]
    / drug_cutoff_results["actual_3m_usage"].replace(0, np.nan)
) * 100

drug_cutoff_results = drug_cutoff_results.sort_values(
    "absolute_error",
    ascending=False
)

display(drug_cutoff_results.head(20))

	drugfull	total_units	actual_3m_usage	predicted_3m_usage	absolute_error	bias	bias_percent
341	drug_2b4752a1e609	381.0	1568.0	1388.362305	179.637695	179.637695	11.456486
185	drug_1b7500ddafc8	16.0	68.0	174.238800	106.238800	-106.238800	-156.233529
1433	drug_a3ffd5cb15c9	72.0	289.0	204.864685	84.135315	84.135315	29.112566
1394	drug_9e65f61d9995	62.0	280.0	209.498734	70.501266	70.501266	25.179024
770	drug_58e790c405a1	96.0	260.0	323.167389	63.167389	-63.167389	-24.295150
497	drug_3eb69a66366e	22.0	122.0	62.211151	59.788849	59.788849	49.007253
380	drug_2fa7ea18cba1	30.0	93.0	50.212574	42.787426	42.787426	46.007985
107	drug_119bcb15803e	4.0	42.0	16.031942	25.968058	25.968058	61.828709
1784	drug_d174b29edcbb	22.0	79.0	102.852882	23.852882	-23.852882	-30.193522
1082	drug_7eae036aa176	14.0	9.0	32.525101	23.525101	-23.525101	-261.390008
1628	drug_b8556e04d039	0.0	31.0	11.264912	19.735088	19.735088	63.661575
1043	drug_7cab3c1e3f92	7.0	54.0	36.406353	17.593647	17.593647	32.580828
1550	drug_b3fdd633217f	20.0	143.0	126.440796	16.559204	16.559204	11.579863
536	drug_3edead72ce43	0.0	28.0	13.484136	14.515864	14.515864	51.842373
731	drug_57f8fe0a5f8a	30.0	96.0	82.630074	13.369926	13.369926	13.927007
848	drug_5df05892de2b	6.0	0.0	11.664601	11.664601	-11.664601	NaN
1667	drug_c6283c502b77	0.0	12.0	0.341338	11.658662	11.658662	97.155514
224	drug_1f2ebad2b6dd	0.0	33.0	23.279350	9.720650	9.720650	29.456514
2018	drug_f80641253af7	0.0	10.0	0.349325	9.650675	9.650675	96.506751
419	drug_37fd409a7a07	8.0	30.0	20.754797	9.245203	9.245203	30.817343

# Plot top 10 drug-level actual vs predicted 3-month usage

top_10_drugs = drug_cutoff_results.head(10)

x = np.arange(len(top_10_drugs))

plt.figure(figsize=(12, 6))

plt.bar(
    x - 0.2,
    top_10_drugs["actual_3m_usage"],
    width=0.4,
    label="Actual March-May usage"
)

plt.bar(
    x + 0.2,
    top_10_drugs["predicted_3m_usage"],
    width=0.4,
    label="Predicted March-May usage"
)

plt.xticks(
    x,
    top_10_drugs["drugfull"],
    rotation=45,
    ha="right"
)

plt.ylabel("3-month usage")
plt.title("Actual vs Predicted 3-Month Usage by Drug")
plt.legend()
plt.tight_layout()
plt.show()

24. Rolling three-month forecast-origin comparison

# Compare 3-month model prediction against actual 3-month usage

quarter_forecast_results = results_by_horizon["3_month"].copy()

quarter_forecast_results["forecast_period"] = (
    quarter_forecast_results["month_period"].dt.strftime("%b %Y")
    + " forecast for next 3 months"
)

quarter_summary = (
    quarter_forecast_results
    .groupby("month_period", as_index=False)
    .agg(
        actual_3m_usage=("future_actual_usage", "sum"),
        predicted_3m_usage=("xgboost_prediction", "sum")
    )
)

quarter_summary["absolute_error"] = (
    quarter_summary["actual_3m_usage"]
    - quarter_summary["predicted_3m_usage"]
).abs()

quarter_summary["bias"] = (
    quarter_summary["actual_3m_usage"]
    - quarter_summary["predicted_3m_usage"]
)

quarter_summary["bias_percent"] = (
    quarter_summary["bias"]
    / quarter_summary["actual_3m_usage"]
) * 100

quarter_summary["wmape"] = (
    quarter_summary["absolute_error"]
    / quarter_summary["actual_3m_usage"]
) * 100

display(quarter_summary)

	month_period	actual_3m_usage	predicted_3m_usage	absolute_error	bias	bias_percent	wmape
0	2025-05-01	4405.0	4049.140625	355.859375	355.859375	8.078533	8.078533
1	2025-06-01	4566.0	4344.816406	221.183594	221.183594	4.844144	4.844144
2	2025-07-01	4853.0	4379.093262	473.906738	473.906738	9.765233	9.765233
3	2025-08-01	4691.0	4290.887207	400.112793	400.112793	8.529371	8.529371
4	2025-09-01	4767.0	3886.316406	880.683594	880.683594	18.474588	18.474588
5	2025-10-01	4831.0	4183.995117	647.004883	647.004883	13.392773	13.392773
6	2025-11-01	4248.0	4161.677734	86.322266	86.322266	2.032068	2.032068

plt.figure(figsize=(10, 6))

plt.plot(
    quarter_summary["month_period"],
    quarter_summary["actual_3m_usage"],
    marker="o",
    label="Actual 3-month usage"
)

plt.plot(
    quarter_summary["month_period"],
    quarter_summary["predicted_3m_usage"],
    marker="o",
    linestyle="--",
    label="Predicted 3-month usage"
)

plt.ylim(0, quarter_summary[["actual_3m_usage", "predicted_3m_usage"]].values.max() * 1.1)
plt.xlim(quarter_summary["month_period"].min() - pd.Timedelta(days=15), quarter_summary["month_period"].max() + pd.Timedelta(days=15))
plt.xlabel("Forecast origin")
plt.ylabel("Total units")
plt.title("Actual vs Predicted 3-Month Usage by Forecast Origin")
plt.legend()
plt.xticks(rotation=45)
plt.tight_layout()
plt.show()