predictive_model_xgboost.ipynb
XGBoost Predictive Model
Imports
import pandas as pd
import numpy as np
import time
import os
from sklearn.model_selection import train_test_split, GridSearchCV, StratifiedKFold, RandomizedSearchCV
from sklearn.preprocessing import StandardScaler
from sklearn.impute import SimpleImputer
from sklearn.pipeline import Pipeline
from sklearn.base import clone
from xgboost import XGBClassifier
from sklearn.metrics import (
roc_auc_score, recall_score, precision_score, accuracy_score,
confusion_matrix, fbeta_score
)
from sklearn.dummy import DummyClassifier
from sklearn.utils import resample
import matplotlib.pyplot as plt
Load tiers
from tiering_preparation import TIERS
tiers = {
"Tier_1_Basic": TIERS["Tier_1_Basic"],
"Tier_2_Clinical": TIERS["Tier_2_Clinical"],
"Tier_3_Personalized": TIERS["Tier_3_Personalized"],
}
from pathlib import Path
# 🔒 HARD-ANCHOR PROJECT ROOT (EDIT THIS PATH ONCE)
PROJECT_ROOT = Path(
r"C:\Users\agila\Documents\STM\Predicting-Past-Year-Major-Depressive-Episode"
)
RESULTS_DIR = PROJECT_ROOT / "results"
FIGURES_DIR = RESULTS_DIR / "figures"
RESULTS_DIR.mkdir(exist_ok=True)
FIGURES_DIR.mkdir(exist_ok=True)
print("Saving results to:", RESULTS_DIR.resolve())
Saving results to: C:\Users\agila\Documents\STM\Predicting-Past-Year-Major-Depressive-Episode\results
Load dataset, drop weight from X, but SAVE w
df = pd.read_pickle("../data/prepared/NSDUH_2023_prepared.pkl")
TARGET = "IRAMDEYR"
WEIGHT = "ANALWT2_C"
y = df[TARGET]
w = df[WEIGHT]
X = df.drop(columns=[TARGET, WEIGHT])
Split train/test (consistent with feature selection)
X_train, X_test, y_train, y_test, w_train, w_test = train_test_split(
X, y, w,
test_size=0.2,
stratify=y,
random_state=42
)
Compute fixed class imbalance
pos = y_train.sum()
neg = len(y_train) - pos
SCALE_POS_WEIGHT = neg / pos
Define metric helper
def evaluate_model(y_true, preds, probs):
tn, fp, fn, tp = confusion_matrix(y_true, preds).ravel()
specificity = tn / (tn + fp)
sensitivity = recall_score(y_true, preds) # same as recall
ppv = precision_score(y_true, preds)
npv = tn / (tn + fn)
return {
"Accuracy": accuracy_score(y_true, preds),
"Sensitivity": sensitivity,
"Specificity": specificity,
"PPV": ppv,
"NPV": npv,
"AUC": roc_auc_score(y_true, probs),
"F1": fbeta_score(y_true, preds, beta=1),
"F2": fbeta_score(y_true, preds, beta=2),
}
XGBoost Pipeline w/ Weighted Fitting + RandomizedSearchCV
param_grid = {
"clf__n_estimators": [50, 100, 150, 200],
"clf__max_depth": [3, 5, 7],
"clf__learning_rate": [0.01, 0.05, 0.1],
"clf__subsample": [0.8, 1.0],
"clf__scale_pos_weight": [1, 9],
}
cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
# RandomizedSearchCV
def build_model():
pipe = Pipeline([
("clf", XGBClassifier(
objective="binary:logistic",
eval_metric="logloss",
enable_categorical=True,
tree_method="hist",
random_state=42,
n_jobs=-1
))
])
model = RandomizedSearchCV(
estimator=pipe,
param_distributions=param_grid,
n_iter=20,
scoring="recall",
cv=cv,
n_jobs=-1,
random_state=42,
verbose=2
)
return model
Bootstrapped Confidence Intervals (ROC + PR)
def bootstrap_auc_ci(model, X, y, n_boot=1000):
aucs = []
for _ in range(n_boot):
X_bs, y_bs = resample(X, y, replace=True)
probs = model.predict_proba(X_bs)[:, 1]
aucs.append(roc_auc_score(y_bs, probs))
return np.percentile(aucs, [2.5, 97.5])
Baseline Model (Dummy Classifier)
# dummy = DummyClassifier(strategy="most_frequent")
# dummy.fit(X_train, y_train)
# dummy_preds = dummy.predict(X_test)
# dummy_probs = np.zeros_like(dummy_preds)
# dummy_metrics = evaluate_model(y_test, dummy_preds, dummy_probs)
Train + Evaluate XGBoost on Each tier
results = []
best_models = {}
for tier_name, features in tiers.items():
print(f"\n----- Training XGBoost on {tier_name} ({len(features)} features) -----")
start_time = time.time()
Xtr = X_train[features]
Xte = X_test[features]
model = build_model()
# Fit WITH survey weights
model.fit(Xtr, y_train, clf__sample_weight=w_train)
# Save best model for fold-wise CV later
best_models[tier_name] = model.best_estimator_
preds = model.predict(Xte)
probs = model.predict_proba(Xte)[:, 1]
metrics = evaluate_model(y_test, preds, probs)
# Bootstrapped CI
ci_low, ci_high = bootstrap_auc_ci(model, Xte, y_test)
metrics["tier"] = tier_name
metrics["threshold"] = 0.50
metrics["n_features"] = len(features)
metrics["AUC_CI_Low"] = ci_low
metrics["AUC_CI_High"] = ci_high
metrics["Best_Params"] = str(model.best_params_)
results.append(metrics)
print(f"Best params: {model.best_params_}")
print(f"Completed in {time.time() - start_time:.1f} seconds")
----- Training XGBoost on Tier_1_Basic (6 features) -----
Fitting 5 folds for each of 20 candidates, totalling 100 fits
Best params: {'clf__subsample': 1.0, 'clf__scale_pos_weight': 9, 'clf__n_estimators': 50, 'clf__max_depth': 3, 'clf__learning_rate': 0.05}
Completed in 77.4 seconds
----- Training XGBoost on Tier_2_Clinical (10 features) -----
Fitting 5 folds for each of 20 candidates, totalling 100 fits
Best params: {'clf__subsample': 1.0, 'clf__scale_pos_weight': 9, 'clf__n_estimators': 50, 'clf__max_depth': 3, 'clf__learning_rate': 0.05}
Completed in 71.7 seconds
----- Training XGBoost on Tier_3_Personalized (19 features) -----
Fitting 5 folds for each of 20 candidates, totalling 100 fits
Best params: {'clf__subsample': 1.0, 'clf__scale_pos_weight': 9, 'clf__n_estimators': 50, 'clf__max_depth': 3, 'clf__learning_rate': 0.05}
Completed in 93.4 seconds
Convert to DataFrame + Export
# 5-Fold Cross-Validation on TRAIN set (no test leakage)
skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
metrics_order = ["Accuracy", "Sensitivity", "Specificity", "PPV", "NPV", "AUC", "F1", "F2"]
all_fold_dfs = []
for tier_name, features in tiers.items():
print(f"\n===== {tier_name} =====")
fold_rows = []
best_model = best_models[tier_name]
X_tier = X_train[features]
w_tier = w_train
for fold_i, (tr_idx, val_idx) in enumerate(skf.split(X_tier, y_train), start=1):
Xtr_fold, Xval_fold = X_tier.iloc[tr_idx], X_tier.iloc[val_idx]
ytr_fold, yval_fold = y_train.iloc[tr_idx], y_train.iloc[val_idx]
wtr_fold = w_tier.iloc[tr_idx]
model_clone = clone(best_model)
model_clone.fit(Xtr_fold, ytr_fold, clf__sample_weight=wtr_fold)
preds = model_clone.predict(Xval_fold)
probs = model_clone.predict_proba(Xval_fold)[:, 1]
fold_metrics = evaluate_model(yval_fold, preds, probs)
fold_metrics["Fold"] = f"Fold {fold_i}"
fold_rows.append(fold_metrics)
# Table: rows=Metric, cols=Fold1..Fold5 (TRANSPOSE like LightGBM)
fold_df = pd.DataFrame(fold_rows).set_index("Fold")[metrics_order].T
fold_df["Mean"] = fold_df.mean(axis=1)
fold_df["Std Dev"] = fold_df.std(axis=1)
fold_df.insert(0, "Tier", tier_name)
fold_df = fold_df.reset_index().rename(columns={"index": "Metric"})
all_fold_dfs.append(fold_df)
print(fold_df)
# Combine and export
cv_results_df = pd.concat(all_fold_dfs, ignore_index=True)
print("\nSaved: xgb_5fold_threshold0.50_by_tier.csv")
===== Tier_1_Basic =====
Fold Metric Tier Fold 1 Fold 2 Fold 3 Fold 4 \
0 Accuracy Tier_1_Basic 0.558941 0.532761 0.601717 0.568526
1 Sensitivity Tier_1_Basic 0.724656 0.783750 0.707500 0.736250
2 Specificity Tier_1_Basic 0.537555 0.500323 0.588045 0.546850
3 PPV Tier_1_Basic 0.168216 0.168548 0.181643 0.173542
4 NPV Tier_1_Basic 0.937993 0.947095 0.939597 0.941324
5 AUC Tier_1_Basic 0.693584 0.709003 0.714747 0.701434
6 F1 Tier_1_Basic 0.273049 0.277434 0.289070 0.280877
7 F2 Tier_1_Basic 0.436125 0.453035 0.448068 0.446618
Fold Fold 5 Mean Std Dev
0 0.624696 0.577328 0.032374
1 0.685857 0.727603 0.032811
2 0.616801 0.557915 0.040569
3 0.187671 0.175924 0.007618
4 0.938314 0.940865 0.003328
5 0.716441 0.707042 0.008533
6 0.294703 0.283027 0.007853
7 0.448005 0.446370 0.005568
===== Tier_2_Clinical =====
Fold Metric Tier Fold 1 Fold 2 Fold 3 Fold 4 \
0 Accuracy Tier_2_Clinical 0.831903 0.830758 0.827754 0.826753
1 Sensitivity Tier_2_Clinical 0.903630 0.895000 0.908750 0.903750
2 Specificity Tier_2_Clinical 0.822646 0.822456 0.817286 0.816801
3 PPV Tier_2_Clinical 0.396703 0.394490 0.391281 0.389338
4 NPV Tier_2_Clinical 0.985106 0.983768 0.985776 0.984999
5 AUC Tier_2_Clinical 0.931716 0.919203 0.929653 0.925441
6 F1 Tier_2_Clinical 0.551355 0.547610 0.547028 0.544223
7 F2 Tier_2_Clinical 0.719697 0.713858 0.718664 0.714851
Fold Fold 5 Mean Std Dev
0 0.836314 0.830696 0.003383
1 0.917397 0.905705 0.007331
2 0.825848 0.821007 0.003457
3 0.404749 0.395312 0.005359
4 0.987254 0.985381 0.001139
5 0.933265 0.927855 0.005061
6 0.561686 0.550380 0.006093
7 0.731975 0.719809 0.006470
===== Tier_3_Personalized =====
Fold Metric Tier Fold 1 Fold 2 Fold 3 \
0 Accuracy Tier_3_Personalized 0.836481 0.824320 0.825608
1 Sensitivity Tier_3_Personalized 0.909887 0.881250 0.923750
2 Specificity Tier_3_Personalized 0.827007 0.816963 0.812924
3 PPV Tier_3_Personalized 0.404338 0.383569 0.389562
4 NPV Tier_3_Personalized 0.986133 0.981561 0.988023
5 AUC Tier_3_Personalized 0.937897 0.925039 0.933029
6 F1 Tier_3_Personalized 0.559877 0.534496 0.548016
7 F2 Tier_3_Personalized 0.727873 0.699682 0.724936
Fold Fold 4 Fold 5 Mean Std Dev
0 0.828755 0.840607 0.831154 0.006339
1 0.901250 0.918648 0.906957 0.014971
2 0.819386 0.830533 0.821363 0.006487
3 0.392061 0.411666 0.396239 0.010258
4 0.984663 0.987514 0.985579 0.002324
5 0.931067 0.939726 0.933352 0.005210
6 0.546419 0.568552 0.551472 0.011730
7 0.715420 0.737096 0.721002 0.012707
Saved: xgb_5fold_threshold0.50_by_tier.csv
import os
cv_results_df.to_csv(
RESULTS_DIR / "xgb_5fold_threshold0.50_by_tier.csv",
index=False
)
print("✅ Saved: ./results/xgb_5fold_threshold0.50_by_tier.csv")
print("Rows:", len(cv_results_df))
print(cv_results_df.head())
✅ Saved: ./results/xgb_5fold_threshold0.50_by_tier.csv
Rows: 24
Fold Metric Tier Fold 1 Fold 2 Fold 3 Fold 4 \
0 Accuracy Tier_1_Basic 0.558941 0.532761 0.601717 0.568526
1 Sensitivity Tier_1_Basic 0.724656 0.783750 0.707500 0.736250
2 Specificity Tier_1_Basic 0.537555 0.500323 0.588045 0.546850
3 PPV Tier_1_Basic 0.168216 0.168548 0.181643 0.173542
4 NPV Tier_1_Basic 0.937993 0.947095 0.939597 0.941324
Fold Fold 5 Mean Std Dev
0 0.624696 0.577328 0.032374
1 0.685857 0.727603 0.032811
2 0.616801 0.557915 0.040569
3 0.187671 0.175924 0.007618
4 0.938314 0.940865 0.003328
results_df = pd.DataFrame(results)
# Reorder columns to match LightGBM format
front_cols = ["tier", "threshold", "n_features"]
other_cols = [c for c in results_df.columns if c not in front_cols]
results_df = results_df[front_cols + other_cols]
results_df.to_csv("./results/xgb_two_rows_per_tier.csv", index=False)
print(results_df)
---------------------------------------------------------------------------
OSError Traceback (most recent call last)
Cell In[67], line 8
5 other_cols = [c for c in results_df.columns if c not in front_cols]
6 results_df = results_df[front_cols + other_cols]
----> 8 results_df.to_csv("./results/xgb_two_rows_per_tier.csv", index=False)
10 print(results_df)
File c:\Users\agila\Documents\STM\Predicting-Past-Year-Major-Depressive-Episode\venv\Lib\site-packages\pandas\util\_decorators.py:333, in deprecate_nonkeyword_arguments.<locals>.decorate.<locals>.wrapper(*args, **kwargs)
327 if len(args) > num_allow_args:
328 warnings.warn(
329 msg.format(arguments=_format_argument_list(allow_args)),
330 FutureWarning,
331 stacklevel=find_stack_level(),
332 )
--> 333 return func(*args, **kwargs)
File c:\Users\agila\Documents\STM\Predicting-Past-Year-Major-Depressive-Episode\venv\Lib\site-packages\pandas\core\generic.py:3989, in NDFrame.to_csv(self, path_or_buf, sep, na_rep, float_format, columns, header, index, index_label, mode, encoding, compression, quoting, quotechar, lineterminator, chunksize, date_format, doublequote, escapechar, decimal, errors, storage_options)
3978 df = self if isinstance(self, ABCDataFrame) else self.to_frame()
3980 formatter = DataFrameFormatter(
3981 frame=df,
3982 header=header,
(...) 3986 decimal=decimal,
3987 )
-> 3989 return DataFrameRenderer(formatter).to_csv(
3990 path_or_buf,
3991 lineterminator=lineterminator,
3992 sep=sep,
3993 encoding=encoding,
3994 errors=errors,
3995 compression=compression,
3996 quoting=quoting,
3997 columns=columns,
3998 index_label=index_label,
3999 mode=mode,
4000 chunksize=chunksize,
4001 quotechar=quotechar,
4002 date_format=date_format,
4003 doublequote=doublequote,
4004 escapechar=escapechar,
4005 storage_options=storage_options,
4006 )
File c:\Users\agila\Documents\STM\Predicting-Past-Year-Major-Depressive-Episode\venv\Lib\site-packages\pandas\io\formats\format.py:1014, in DataFrameRenderer.to_csv(self, path_or_buf, encoding, sep, columns, index_label, mode, compression, quoting, quotechar, lineterminator, chunksize, date_format, doublequote, escapechar, errors, storage_options)
993 created_buffer = False
995 csv_formatter = CSVFormatter(
996 path_or_buf=path_or_buf,
997 lineterminator=lineterminator,
(...) 1012 formatter=self.fmt,
1013 )
-> 1014 csv_formatter.save()
1016 if created_buffer:
1017 assert isinstance(path_or_buf, StringIO)
File c:\Users\agila\Documents\STM\Predicting-Past-Year-Major-Depressive-Episode\venv\Lib\site-packages\pandas\io\formats\csvs.py:251, in CSVFormatter.save(self)
247 """
248 Create the writer & save.
249 """
250 # apply compression and byte/text conversion
--> 251 with get_handle(
252 self.filepath_or_buffer,
253 self.mode,
254 encoding=self.encoding,
255 errors=self.errors,
256 compression=self.compression,
257 storage_options=self.storage_options,
258 ) as handles:
259 # Note: self.encoding is irrelevant here
260 self.writer = csvlib.writer(
261 handles.handle,
262 lineterminator=self.lineterminator,
(...) 267 quotechar=self.quotechar,
268 )
270 self._save()
File c:\Users\agila\Documents\STM\Predicting-Past-Year-Major-Depressive-Episode\venv\Lib\site-packages\pandas\io\common.py:749, in get_handle(path_or_buf, mode, encoding, compression, memory_map, is_text, errors, storage_options)
747 # Only for write methods
748 if "r" not in mode and is_path:
--> 749 check_parent_directory(str(handle))
751 if compression:
752 if compression != "zstd":
753 # compression libraries do not like an explicit text-mode
File c:\Users\agila\Documents\STM\Predicting-Past-Year-Major-Depressive-Episode\venv\Lib\site-packages\pandas\io\common.py:616, in check_parent_directory(path)
614 parent = Path(path).parent
615 if not parent.is_dir():
--> 616 raise OSError(rf"Cannot save file into a non-existent directory: '{parent}'")
OSError: Cannot save file into a non-existent directory: 'results'
SHAP Analysis
# SHAP Analysis for XGBoost
import shap
import os
os.makedirs("./results/figures", exist_ok=True)
shap_importance_all = []
for tier_name, features in tiers.items():
print(f"\n===== {tier_name} =====")
best_model = best_models[tier_name]
clf = best_model.named_steps["clf"] # only step now
# Use raw test data (NO IMPUTATION)
X_tier_test = X_test[features]
# TreeExplainer for XGBoost
explainer = shap.TreeExplainer(clf)
shap_values = explainer.shap_values(X_tier_test)
# --- Global importance (bar plot) ---
plt.figure(figsize=(10, 6))
shap.summary_plot(
shap_values,
X_tier_test,
plot_type="bar",
show=False
)
plt.title(f"SHAP Feature Importance - {tier_name}")
plt.tight_layout()
plt.savefig(
f"./results/figures/xgb_shap_importance_{tier_name}.png",
dpi=300,
bbox_inches="tight"
)
plt.show()
# --- Beeswarm plot ---
plt.figure(figsize=(10, 6))
shap.summary_plot(
shap_values,
X_tier_test,
show=False
)
plt.title(f"SHAP Summary - {tier_name}")
plt.tight_layout()
plt.savefig(
f"./results/figures/xgb_shap_beeswarm_{tier_name}.png",
dpi=300,
bbox_inches="tight"
)
plt.show()
# --- Feature importance table ---
mean_abs_shap = np.abs(shap_values).mean(axis=0)
shap_importance = pd.DataFrame({
"feature": features,
"mean_abs_shap": mean_abs_shap
}).sort_values("mean_abs_shap", ascending=False)
shap_importance["tier"] = tier_name
shap_importance["rank"] = range(1, len(features) + 1)
shap_importance_all.append(shap_importance)
print("\nTop 10 Features by SHAP:")
print(shap_importance.head(10).to_string(index=False))
# --- Dependence plots (top 3 features) ---
top_features = shap_importance["feature"].head(3).tolist()
for i, feat in enumerate(top_features, start=1):
plt.figure(figsize=(8, 5))
shap.dependence_plot(
feat,
shap_values,
X_tier_test,
show=False
)
plt.title(f"SHAP Dependence: {feat} - {tier_name}")
plt.tight_layout()
plt.savefig(
f"./results/figures/xgb_shap_dependence_{tier_name}_top{i}_{feat}.png",
dpi=300,
bbox_inches="tight"
)
plt.show()
# --- Save SHAP values ---
shap_df = pd.DataFrame(shap_values, columns=features, index=X_tier_test.index)
shap_df.to_csv(
f"./results/xgb_shap_values_{tier_name}.csv"
)
# --- Combined importance table ---
shap_importance_combined = pd.concat(shap_importance_all, ignore_index=True)
shap_importance_combined = shap_importance_combined[
["tier", "rank", "feature", "mean_abs_shap"]
]
shap_importance_combined.to_csv(
"./results/xgb_shap_importance_by_tier.csv",
index=False
)
print("\nSHAP analysis complete")
===== Tier_1_Basic =====
Top 10 Features by SHAP:
feature mean_abs_shap tier rank
CATAG3 0.341068 Tier_1_Basic 1
SVYRDUDANY 0.262221 Tier_1_Basic 2
IRMARIT_ABS 0.247915 Tier_1_Basic 3
IRSEX 0.227461 Tier_1_Basic 4
IRPINC3_ABS 0.071741 Tier_1_Basic 5
WRKDRGHLP_ABS 0.010195 Tier_1_Basic 6
<Figure size 800x500 with 0 Axes>
<Figure size 800x500 with 0 Axes>
<Figure size 800x500 with 0 Axes>
===== Tier_2_Clinical =====
Top 10 Features by SHAP:
feature mean_abs_shap tier rank
KSSLR6MAX_ABS 1.625626 Tier_2_Clinical 1
IRIMPGOUT_ABS 0.334070 Tier_2_Clinical 2
RCVYMHPRB 0.173065 Tier_2_Clinical 3
IRSUICTHNK 0.093117 Tier_2_Clinical 4
IRSEX 0.015487 Tier_2_Clinical 5
IRMARIT_ABS 0.014273 Tier_2_Clinical 6
WRKDRGHLP_ABS 0.004962 Tier_2_Clinical 7
CATAG3 0.003675 Tier_2_Clinical 8
SVYRDUDANY 0.000000 Tier_2_Clinical 9
IRPINC3_ABS 0.000000 Tier_2_Clinical 10
<Figure size 800x500 with 0 Axes>
<Figure size 800x500 with 0 Axes>
<Figure size 800x500 with 0 Axes>
===== Tier_3_Personalized =====
Top 10 Features by SHAP:
feature mean_abs_shap tier rank
KSSLR6MAX_ABS 1.357722 Tier_3_Personalized 1
WHODASDASC_ABS 0.653802 Tier_3_Personalized 2
RCVYMHPRB 0.158053 Tier_3_Personalized 3
IRIMPGOUT_ABS 0.090075 Tier_3_Personalized 4
IRSUICTHNK 0.089366 Tier_3_Personalized 5
IRSEX 0.018302 Tier_3_Personalized 6
COCLALCUSE_ABS 0.013641 Tier_3_Personalized 7
IRMARIT_ABS 0.009763 Tier_3_Personalized 8
IRPINC3_ABS 0.006501 Tier_3_Personalized 9
DIFOBTCRK 0.006431 Tier_3_Personalized 10
<Figure size 800x500 with 0 Axes>
<Figure size 800x500 with 0 Axes>
<Figure size 800x500 with 0 Axes>
SHAP analysis complete