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.ensemble import RandomForestClassifier
from sklearn.base import clone
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
import os
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"],
}
Saved: C:\Users\agila\Documents\STM\Predicting-Past-Year-Major-Depressive-Episode\tiers\tier_1_basic_features.csv
Saved: C:\Users\agila\Documents\STM\Predicting-Past-Year-Major-Depressive-Episode\tiers\tier_2_clinical_features.csv
Saved: C:\Users\agila\Documents\STM\Predicting-Past-Year-Major-Depressive-Episode\tiers\tier_3_personalized_features.csv
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"
# Encode ABS categorical variables
ABS_CATEGORICAL = [
"IRMARIT_ABS",
"IRPINC3_ABS",
"WRKDRGHLP_ABS",
"KSSLR6MAX_ABS",
"IRIMPGOUT_ABS",
"WHODASDASC_ABS",
"COCLALCUSE_ABS",
]
ABS_CATEGORY_MAPS = {}
for col in ABS_CATEGORICAL:
if col in df.columns:
df[col] = df[col].astype("category")
ABS_CATEGORY_MAPS[col] = dict(enumerate(df[col].cat.categories))
df[col] = df[col].cat.codes
# 🔥 REBUILD X AFTER ENCODING
y = df[TARGET]
w = df[WEIGHT]
X = df.drop(columns=[TARGET, WEIGHT])
print(f'Total samples: {X.shape[0]}')
Total samples: 43687
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
)
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),
}
Random Forest Pipeline w/ Weighted Fitting + RandomizedSearchCV
param_grid = {
"clf__n_estimators": [50, 100, 150, 200],
"clf__max_depth": [5, 10, 15],
"clf__min_samples_split": [2, 5],
"clf__class_weight": [None, "balanced"],
}
cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
# RandomizedSearchCV
def build_model():
pipe = Pipeline([
("clf", RandomForestClassifier(
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 Random Forest on Each tier
results = []
best_models = {}
for tier_name, features in tiers.items():
print(f"\n----- Training Random Forest on {tier_name} ({len(features)} features) -----")
start_time = time.time()
Xtr = X_train[features]
Xte = X_test[features]
model = build_model()
bad_cols = Xtr.select_dtypes(include=["object", "string", "category"]).columns.tolist()
if bad_cols:
raise ValueError(
f"[{tier_name}] Non-numeric columns found: {bad_cols}"
)
# 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 Random Forest on Tier_1_Basic (6 features) -----
Fitting 5 folds for each of 20 candidates, totalling 100 fits
Best params: {'clf__n_estimators': 150, 'clf__min_samples_split': 2, 'clf__max_depth': 5, 'clf__class_weight': 'balanced'}
Completed in 184.7 seconds
----- Training Random Forest on Tier_2_Clinical (10 features) -----
Fitting 5 folds for each of 20 candidates, totalling 100 fits
Best params: {'clf__n_estimators': 50, 'clf__min_samples_split': 2, 'clf__max_depth': 5, 'clf__class_weight': 'balanced'}
Completed in 171.8 seconds
----- Training Random Forest on Tier_3_Personalized (19 features) -----
Fitting 5 folds for each of 20 candidates, totalling 100 fits
Best params: {'clf__n_estimators': 100, 'clf__min_samples_split': 2, 'clf__max_depth': 5, 'clf__class_weight': 'balanced'}
Completed in 179.4 seconds
# 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)
assert not X_tier.select_dtypes(include=["object", "category"]).any().any(), \
f"{tier_name} still contains categorical dtypes"
if hasattr(model_clone, "named_steps"):
model_clone.fit(Xtr_fold, ytr_fold, clf__sample_weight=wtr_fold)
else:
model_clone.fit(Xtr_fold, ytr_fold, 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)
cv_results_df.to_csv(
RESULTS_DIR / "rf_5fold_threshold0.50_by_tier.csv",
index=False
)
print("\nSaved: rf_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.625036 0.661087 0.675536 0.655365
1 Sensitivity Tier_1_Basic 0.628285 0.636250 0.620000 0.622500
2 Specificity Tier_1_Basic 0.624616 0.664297 0.682714 0.659612
3 PPV Tier_1_Basic 0.177636 0.196753 0.201626 0.191171
4 NPV Tier_1_Basic 0.928674 0.933909 0.932892 0.931129
5 AUC Tier_1_Basic 0.692042 0.711625 0.718023 0.701979
6 F1 Tier_1_Basic 0.276966 0.300561 0.304294 0.292511
7 F2 Tier_1_Basic 0.416805 0.439779 0.438163 0.428941
Fold Fold 5 Mean Std Dev
0 0.682787 0.659962 0.020029
1 0.613267 0.624060 0.007770
2 0.691761 0.664600 0.023193
3 0.204337 0.194305 0.009465
4 0.932694 0.931860 0.001824
5 0.715546 0.707843 0.009606
6 0.306537 0.296174 0.010724
7 0.437969 0.432331 0.008647
===== Tier_2_Clinical =====
Fold Metric Tier Fold 1 Fold 2 Fold 3 Fold 4 \
0 Accuracy Tier_2_Clinical 0.833619 0.826753 0.834621 0.821888
1 Sensitivity Tier_2_Clinical 0.901126 0.893750 0.898750 0.902500
2 Specificity Tier_2_Clinical 0.824907 0.818094 0.826333 0.811470
3 PPV Tier_2_Clinical 0.399113 0.388376 0.400780 0.382213
4 NPV Tier_2_Clinical 0.984767 0.983492 0.984411 0.984709
5 AUC Tier_2_Clinical 0.931201 0.917699 0.928792 0.924019
6 F1 Tier_2_Clinical 0.553208 0.541462 0.554356 0.537003
7 F2 Tier_2_Clinical 0.720000 0.709185 0.719864 0.709373
Fold Fold 5 Mean Std Dev
0 0.840607 0.831498 0.006513
1 0.908636 0.900952 0.004862
2 0.831826 0.822526 0.007050
3 0.410866 0.396270 0.010014
4 0.986021 0.984680 0.000811
5 0.931139 0.926570 0.005147
6 0.565861 0.550378 0.010217
7 0.731412 0.717967 0.008240
===== Tier_3_Personalized =====
Fold Metric Tier Fold 1 Fold 2 Fold 3 \
0 Accuracy Tier_3_Personalized 0.841917 0.836624 0.833047
1 Sensitivity Tier_3_Personalized 0.904881 0.866250 0.913750
2 Specificity Tier_3_Personalized 0.833791 0.832795 0.822617
3 PPV Tier_3_Personalized 0.412671 0.401042 0.399672
4 NPV Tier_3_Personalized 0.985491 0.979666 0.986630
5 AUC Tier_3_Personalized 0.934691 0.922078 0.931718
6 F1 Tier_3_Personalized 0.566837 0.548259 0.556105
7 F2 Tier_3_Personalized 0.730598 0.703125 0.726785
Fold Fold 4 Fold 5 Mean Std Dev
0 0.834907 0.845185 0.838336 0.004525
1 0.903750 0.906133 0.898953 0.016722
2 0.826010 0.837318 0.830506 0.005383
3 0.401667 0.418255 0.406661 0.007429
4 0.985164 0.985736 0.984537 0.002484
5 0.932145 0.939194 0.931965 0.005613
6 0.556154 0.572332 0.559937 0.008563
7 0.723000 0.734727 0.723647 0.010977
Saved: rf_5fold_threshold0.50_by_tier.csv
Convert to DataFrame + Export
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_DIR / "rf_two_rows_per_tier.csv",
index=False
)
print(results_df)
tier threshold n_features Accuracy Sensitivity \
0 Tier_1_Basic 0.5 6 0.682879 0.602603
1 Tier_2_Clinical 0.5 10 0.835203 0.874875
2 Tier_3_Personalized 0.5 19 0.834173 0.885886
Specificity PPV NPV AUC F1 F2 AUC_CI_Low \
0 0.693242 0.202285 0.931100 0.711132 0.302893 0.431727 0.694228
1 0.830081 0.399269 0.980913 0.923757 0.548306 0.706548 0.915998
2 0.827497 0.398649 0.982510 0.931227 0.549860 0.711873 0.924279
AUC_CI_High Best_Params
0 0.726792 {'clf__n_estimators': 150, 'clf__min_samples_s...
1 0.930662 {'clf__n_estimators': 50, 'clf__min_samples_sp...
2 0.938016 {'clf__n_estimators': 100, 'clf__min_samples_s...