{ "cells": [ { "cell_type": "markdown", "id": "85d28757", "metadata": { "id": "85d28757" }, "source": [ "\n", "# Optuna‑tuned XGBoost + SHAP — **2018 data, 80/20 split**\n", "\n", "**What this notebook does** \n", "- Uses sample data monitored in 2018 and performs a **chronological 80/20** train/test split. \n", "- Tunes hyperparameters per band using Optuna. \n", "- Trains a full-band model with tuned params. \n", "- Produces SHAP summary/bar/dependence/waterfall plots.\n", "\n", "**Inputs**\n", "- CSV path (default): `Sample_data.csv` — update `DATA_PATH` if your filename differs.\n", "\n", "**Outputs**\n", "- `optuna_results_table2.csv` \n", "- `optuna_rmse_by_band.csv`\n" ] }, { "cell_type": "code", "execution_count": null, "id": "4602e465", "metadata": { "id": "4602e465" }, "outputs": [], "source": [ "\n", "# --- (optional) installs ---\n", "# If packages are missing, uncomment the line below:\n", "# !pip install xgboost shap optuna scikit-learn pandas matplotlib seaborn\n" ] }, { "cell_type": "code", "execution_count": null, "id": "ecca1641", "metadata": { "id": "ecca1641" }, "outputs": [], "source": [ "\n", "# ======================================================================\n", "# Optuna-tuned XGBoost + SHAP (80/20 split)\n", "# ======================================================================\n", "import warnings\n", "warnings.filterwarnings(\"ignore\")\n", "\n", "import pandas as pd\n", "import numpy as np\n", "import optuna\n", "from xgboost import XGBRegressor\n", "from sklearn.metrics import mean_squared_error\n", "import shap\n", "import matplotlib.pyplot as plt\n", "import seaborn as sns\n", "from IPython.display import display\n", "\n", "pd.set_option(\"display.max_columns\", 200)\n", "sns.set(context=\"notebook\", style=\"whitegrid\")\n" ] }, { "cell_type": "code", "execution_count": null, "id": "8b87d19e", "metadata": { "id": "8b87d19e" }, "outputs": [], "source": [ "\n", "# ------------------ CONFIG ------------------\n", "# Default to the uploaded 2018-only CSV\n", "DATA_PATH = \"Sample_data.csv\"\n", "\n", "# Targets and band columns (Refer to Eq 4, Section 2.4)\n", "COLS = {\n", " \"targets\": {\n", " \"SPL10-300 Hz\": \"SPL_10_300\",\n", " \"SPL300-3000 Hz\": \"SPL_300_3000\",\n", " \"SPL3000-24000 Hz\": \"SPL_3000_24000\",\n", " \"SPL10-24000 Hz\": \"SPL_full\",\n", " },\n", " \"bands\": {\n", " \"low\": \"SPL_10_300\",\n", " \"mid\": \"SPL_300_3000\",\n", " \"high\": \"SPL_3000_24000\",\n", " \"full\": \"SPL_full\",\n", " },\n", " # Base predictors consistent with Methods (one-hot for categoricals; others native units)\n", " \"base\": [\"temperature\",\"season\",\"tide\",\"moon_phase\",\"month\",\"hour\",\"H\",\"SE\"],\n", " \"categoricals\": [\"season\",\"tide\",\"moon_phase\",\"month\",\"hour\"],\n", "\n", " # Kept for reference; **not used** for splitting in this notebook\n", " \"year\": \"Year\"\n", "}\n", "\n", "# Broadband vs off-band control for band models\n", "INCLUDE_BROADBAND_CONTROL = True\n", "\n", "# Optuna settings\n", "N_TRIALS = 60\n", "RANDOM_STATE = 42\n", "np.random.seed(RANDOM_STATE)\n", "# --------------------------------------------\n" ] }, { "cell_type": "code", "execution_count": null, "id": "5d82168a", "metadata": { "id": "5d82168a" }, "outputs": [], "source": [ "\n", "# ---------------- LOAD DATA -----------------\n", "df = pd.read_csv(DATA_PATH)\n", "print(\"Data shape:\", df.shape)\n", "display(df.head())\n" ] }, { "cell_type": "code", "execution_count": null, "id": "b36a1885", "metadata": { "id": "b36a1885" }, "outputs": [], "source": [ "\n", "def time_split_80_20(frame: pd.DataFrame):\n", " \"\"\"Chronological 80/20 split (no shuffling).\"\"\"\n", " n = len(frame)\n", " cut = int(n * 0.8)\n", " train = frame.iloc[:cut].copy()\n", " test = frame.iloc[cut:].copy()\n", " return train, test\n", "\n", "def predictors_for_target(label: str):\n", " base = COLS[\"base\"].copy()\n", " b = COLS[\"bands\"]\n", " if label == \"SPL10-24000 Hz\": # full band uses sub-bands + base\n", " return [b[\"low\"], b[\"mid\"], b[\"high\"]] + base\n", " # band models: broadband control or off-band SPLs\n", " if label == \"SPL10-300 Hz\":\n", " return ([b[\"full\"]] if INCLUDE_BROADBAND_CONTROL else [b[\"mid\"], b[\"high\"]]) + base\n", " if label == \"SPL300-3000 Hz\":\n", " return ([b[\"full\"]] if INCLUDE_BROADBAND_CONTROL else [b[\"low\"], b[\"high\"]]) + base\n", " if label == \"SPL3000-24000 Hz\":\n", " return ([b[\"full\"]] if INCLUDE_BROADBAND_CONTROL else [b[\"low\"], b[\"mid\"]]) + base\n", " return base\n", "\n", "def prepare_xy(frame: pd.DataFrame, target_col: str, pred_cols: list, fit_columns=None, return_columns=False):\n", " \"\"\"\n", " Builds X, y. If fit_columns is provided, aligns dummy columns to that schema.\n", " Set return_columns=True to also return the column list used for X.\n", " \"\"\"\n", " X = frame[pred_cols].copy()\n", " y = frame[target_col].copy()\n", " cats = [c for c in COLS[\"categoricals\"] if c in X.columns]\n", " X = pd.get_dummies(X, columns=cats, drop_first=True) # one-hot categoricals\n", "\n", " if fit_columns is not None:\n", " # Add any missing columns as zeros\n", " missing = [c for c in fit_columns if c not in X.columns]\n", " for c in missing:\n", " X[c] = 0\n", " # Drop any unexpected columns\n", " extra = [c for c in X.columns if c not in fit_columns]\n", " if extra:\n", " X = X.drop(columns=extra)\n", " # Reorder to match the fitted schema\n", " X = X[fit_columns]\n", "\n", " if return_columns:\n", " return X, y, list(X.columns)\n", " return X, y\n", "\n", "def objective_factory(X_tr, y_tr, X_va, y_va):\n", " def objective(trial):\n", " params = {\n", " \"n_estimators\": trial.suggest_int(\"n_estimators\", 80, 300),\n", " \"max_depth\": trial.suggest_int(\"max_depth\", 3, 10),\n", " \"learning_rate\": trial.suggest_float(\"learning_rate\", 1e-3, 0.3, log=True),\n", " \"min_child_weight\": trial.suggest_int(\"min_child_weight\", 1, 10),\n", " \"subsample\": trial.suggest_float(\"subsample\", 0.5, 1.0),\n", " \"colsample_bytree\": 1.0,\n", " \"objective\": \"reg:squarederror\",\n", " \"n_jobs\": -1,\n", " \"tree_method\": \"hist\",\n", " \"random_state\": RANDOM_STATE\n", " }\n", " mdl = XGBRegressor(**params)\n", " mdl.fit(\n", " X_tr, y_tr,\n", " eval_set=[(X_va, y_va)],\n", " eval_metric=\"rmse\",\n", " verbose=False,\n", " early_stopping_rounds=50\n", " )\n", " pred = mdl.predict(X_va)\n", " rmse = mean_squared_error(y_va, pred, squared=False)\n", " return rmse\n", " return objective\n", "\n", "def run_band(tgt_label: str, data: pd.DataFrame):\n", " target_col = COLS[\"targets\"][tgt_label]\n", " pred_cols = predictors_for_target(tgt_label)\n", " needed = [target_col] + pred_cols\n", " frame = data.dropna(subset=[c for c in needed if c in data.columns]).copy()\n", "\n", " # 80/20 chronological split\n", " train_df, test_df = time_split_80_20(frame)\n", "\n", " # validation slice from the end of the TRAIN period for early stopping\n", " cut = max(int(len(train_df) * 0.9), len(train_df) - 1000)\n", " inner_tr = train_df.iloc[:cut].copy()\n", " inner_va = train_df.iloc[cut:].copy()\n", "\n", " # Prepare aligned design matrices\n", " X_tr, y_tr, cols = prepare_xy(inner_tr, target_col, pred_cols, return_columns=True)\n", " X_va, y_va = prepare_xy(inner_va, target_col, pred_cols, fit_columns=cols)\n", " X_te, y_te = prepare_xy(test_df, target_col, pred_cols, fit_columns=cols)\n", "\n", " study = optuna.create_study(direction=\"minimize\")\n", " study.optimize(objective_factory(X_tr, y_tr, X_va, y_va), n_trials=N_TRIALS)\n", " best = study.best_params\n", "\n", " # Refit on full training (aligned to train columns), evaluate on 20% test\n", " X_full, y_full = prepare_xy(train_df, target_col, pred_cols, fit_columns=cols)\n", " mdl = XGBRegressor(objective=\"reg:squarederror\",\n", " n_jobs=-1, tree_method=\"hist\",\n", " random_state=RANDOM_STATE, **best)\n", " mdl.fit(X_full, y_full, verbose=False)\n", " rmse_test = mean_squared_error(y_te, mdl.predict(X_te), squared=False)\n", "\n", " return {\n", " \"n_estimators\": int(best[\"n_estimators\"]),\n", " \"max_depth\": int(best[\"max_depth\"]),\n", " \"learning_rate\": float(round(best[\"learning_rate\"], 3)),\n", " \"min_child_weight\": int(best[\"min_child_weight\"]),\n", " \"subsample\": float(round(best[\"subsample\"], 2)),\n", " \"RMSE_test\": float(round(rmse_test, 2))\n", " }\n" ] }, { "cell_type": "code", "execution_count": null, "id": "ee64c998", "metadata": { "id": "ee64c998" }, "outputs": [], "source": [ "\n", "# ---------------- RUN OPTUNA ----------------\n", "bands = [\"SPL10-300 Hz\",\"SPL300-3000 Hz\",\"SPL3000-24000 Hz\",\"SPL10-24000 Hz\"]\n", "\n", "header = {\n", " \"Hyperparameter\": [\"n_estimators\",\"max_depth\",\"learning_rate\",\"min_child_weight\",\"subsample\"],\n", " \"Meaning\": [\"Number of trees\",\"Tree max depth\",\"Step size at each iteration\",\n", " \"Minimum sum of instance weight\",\"Subsample ratio of training instances\"],\n", " \"Search space\": [\"80–300\",\"3–10\",\"0.001–0.3\",\"1–10\",\"0.5–1.0\"]\n", "}\n", "table2 = pd.DataFrame(header)\n", "\n", "rmse_rows = []\n", "results = {}\n", "for label in bands:\n", " print(f\"Tuning for {label} ...\")\n", " res = run_band(label, df)\n", " results[label] = res\n", " table2[label] = [\n", " res[\"n_estimators\"], res[\"max_depth\"], res[\"learning_rate\"],\n", " res[\"min_child_weight\"], res[\"subsample\"]\n", " ]\n", " rmse_rows.append({\"Band\": label, \"RMSE_test\": res[\"RMSE_test\"]})\n", "\n", "display(table2)\n", "rmse_df = pd.DataFrame(rmse_rows)\n", "display(rmse_df)\n", "\n", "# Save outputs for manuscript Table 2\n", "table2.to_csv(\"optuna_results_table2.csv\", index=False)\n", "rmse_df.to_csv(\"optuna_rmse_by_band.csv\", index=False)\n", "print(\"Saved: optuna_results_table2.csv, optuna_rmse_by_band.csv\")\n" ] }, { "cell_type": "code", "execution_count": null, "id": "794681b3", "metadata": { "id": "794681b3" }, "outputs": [], "source": [ "\n", "# ------- Train full-band model with tuned params -------\n", "\n", "try:\n", " opt_df = pd.read_csv(\"optuna_results_table2.csv\")\n", " full_row = opt_df[[\"Hyperparameter\",\"SPL10-24000 Hz\"]].set_index(\"Hyperparameter\")[\"SPL10-24000 Hz\"]\n", " best_full = {\n", " \"n_estimators\": int(full_row[\"n_estimators\"]),\n", " \"max_depth\": int(full_row[\"max_depth\"]),\n", " \"learning_rate\": float(full_row[\"learning_rate\"]),\n", " \"min_child_weight\": int(full_row[\"min_child_weight\"]),\n", " \"subsample\": float(full_row[\"subsample\"]),\n", " }\n", "except Exception:\n", " # Fallback (manuscript Table 2)\n", " best_full = {\"n_estimators\": 157, \"max_depth\": 5, \"learning_rate\": 0.09,\n", " \"min_child_weight\": 3, \"subsample\": 0.60}\n", "\n", "pred_cols_full = predictors_for_target(\"SPL10-24000 Hz\")\n", "# Build full design matrix ON ALL ROWS for SHAP and final model\n", "X_all, y_all = prepare_xy(df, COLS[\"targets\"][\"SPL10-24000 Hz\"], pred_cols_full)\n", "\n", "full_mdl = XGBRegressor(objective=\"reg:squarederror\",\n", " n_estimators=best_full[\"n_estimators\"],\n", " max_depth=best_full[\"max_depth\"],\n", " learning_rate=best_full[\"learning_rate\"],\n", " min_child_weight=best_full[\"min_child_weight\"],\n", " subsample=best_full[\"subsample\"],\n", " colsample_bytree=1.0,\n", " random_state=RANDOM_STATE,\n", " tree_method=\"hist\")\n", "full_mdl.fit(X_all, y_all, verbose=False)\n", "print(\"Full-band model trained with params:\", best_full)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "c400eacc", "metadata": { "id": "c400eacc" }, "outputs": [], "source": [ "\n", "# ---------------------- SHAP --------------------------\n", "explainer = shap.TreeExplainer(full_mdl)\n", "shap_values = explainer.shap_values(X_all)\n", "\n", "# Rank features by mean |SHAP|\n", "mean_abs = np.mean(np.abs(shap_values), axis=0)\n", "ranked_idx = np.argsort(mean_abs)[::-1]\n", "top_k = 9\n", "top_feats = X_all.columns[ranked_idx][:top_k]\n", "print(\"Top features:\", list(top_feats))\n", "\n", "# SHAP summary beeswarm\n", "plt.figure(figsize=(10,6))\n", "shap.summary_plot(shap_values, X_all, plot_type=\"dot\", max_display=top_k, show=False)\n", "plt.title(\"SHAP summary (full band)\")\n", "plt.tight_layout()\n", "plt.show()\n", "\n", "# SHAP bar plot\n", "plt.figure(figsize=(8,5))\n", "shap.summary_plot(shap_values, X_all, plot_type=\"bar\", max_display=top_k, show=False)\n", "plt.title(\"Mean |SHAP| (top features)\")\n", "plt.tight_layout()\n", "plt.show()\n", "\n", "# SHAP dependence plots for top features (adjust how many to render)\n", "for feat in list(top_feats)[:6]:\n", " plt.figure(figsize=(6,4))\n", " shap.dependence_plot(feat, shap_values, X_all, show=False)\n", " plt.title(f\"SHAP dependence: {feat}\")\n", " plt.tight_layout()\n", " plt.show()\n", "\n", "# Waterfall for one observation (row 0); fallback if SHAP version differs\n", "try:\n", " expl = shap.Explanation(values=shap_values[0,:],\n", " base_values=explainer.expected_value,\n", " data=X_all.iloc[0,:].values,\n", " feature_names=list(X_all.columns))\n", " shap.plots.waterfall(expl, max_display=15)\n", "except Exception:\n", " print(\"Waterfall not supported by this SHAP version; showing top contributions for row 0:\")\n", " contrib = pd.Series(shap_values[0,:], index=X_all.columns)\\\n", " .sort_values(key=np.abs, ascending=False)[:15]\n", " display(contrib.to_frame(\"SHAP\").style.background_gradient(axis=None, cmap=\"coolwarm\"))\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "name": "python", "pygments_lexer": "ipython3" }, "colab": { "provenance": [] } }, "nbformat": 4, "nbformat_minor": 5 }