Imports
import pandas as pd
import numpy as np
import time
from sklearn.model_selection import train_test_split, GridSearchCV, StratifiedKFold
from sklearn.preprocessing import StandardScaler
from sklearn.impute import SimpleImputer
from sklearn.pipeline import Pipeline
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import (
roc_auc_score, recall_score, precision_score, accuracy_score,
confusion_matrix, fbeta_score, precision_recall_curve, roc_curve
)
from sklearn.dummy import DummyClassifier
from sklearn.utils import resample
import matplotlib.pyplot as plt
from sklearn.model_selection import cross_validate
from sklearn.metrics import make_scorer
import os
from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import OneHotEncoder
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
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])
categorical_features = X.select_dtypes(
include=["object", "category"]
).columns.tolist()
numeric_features = X.select_dtypes(
include=["int64", "float64"]
).columns.tolist()
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)
npv = tn / (tn + fn) if (tn + fn) > 0 else 0
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),
}
CV metric helpers
def specificity_cv(y_true, y_pred):
tn, fp, fn, tp = confusion_matrix(y_true, y_pred).ravel()
return tn / (tn + fp) if (tn + fp) > 0 else 0
def npv_cv(y_true, y_pred):
tn, fp, fn, tp = confusion_matrix(y_true, y_pred).ravel()
return tn / (tn + fn) if (tn + fn) > 0 else 0
CV scoring dictionary
cv_scoring = {
"accuracy": make_scorer(accuracy_score),
"sensitivity": make_scorer(recall_score),
"specificity": make_scorer(specificity_cv),
"ppv": make_scorer(precision_score, zero_division=0),
"npv": make_scorer(npv_cv),
"f1": make_scorer(fbeta_score, beta=1),
"f2": make_scorer(fbeta_score, beta=2),
"auc": "roc_auc"
}
Logistic Regression Pipeline w/ Weighted Fitting + GridSearchCV
param_grid = {
"clf__C": [0.01, 0.1, 1, 10],
"clf__class_weight": [None, "balanced"],
}
cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import OneHotEncoder
def build_model(numeric_features, categorical_features):
numeric_transformer = Pipeline([
("imputer", SimpleImputer(strategy="median")),
("scaler", StandardScaler())
])
categorical_transformer = Pipeline([
("imputer", SimpleImputer(strategy="most_frequent")),
("onehot", OneHotEncoder(
handle_unknown="ignore",
sparse_output=False
))
])
preprocessor = ColumnTransformer(
transformers=[
("num", numeric_transformer, numeric_features),
("cat", categorical_transformer, categorical_features),
]
)
pipe = Pipeline([
("preprocessor", preprocessor),
("clf", LogisticRegression(
max_iter=3000,
penalty="l2",
solver="lbfgs"
))
])
model = GridSearchCV(
estimator=pipe,
param_grid=param_grid,
scoring="recall",
cv=cv,
n_jobs=-1
)
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)
c:\Users\agila\Documents\STM\Predicting-Past-Year-Major-Depressive-Episode\venv\Lib\site-packages\sklearn\metrics\_classification.py:1731: UndefinedMetricWarning: Precision is ill-defined and being set to 0.0 due to no predicted samples. Use `zero_division` parameter to control this behavior.
_warn_prf(average, modifier, f"{metric.capitalize()} is", result.shape[0])
Threshold tuning helper function
def evaluate_at_threshold(y_true, probs, threshold):
preds = (probs >= threshold).astype(int)
return evaluate_model(y_true, preds, probs)
Train + Evaluate Logistic Regression on Each tier
results = []
for tier_name, features in tiers.items():
print(f"\n----- Training Logistic Regression on {tier_name} -----")
Xtr = X_train[features]
Xte = X_test[features]
# Identify feature types for THIS tier only
categorical_features = Xtr.select_dtypes(
include=["object", "category"]
).columns.tolist()
numeric_features = Xtr.select_dtypes(
include=["int64", "float64"]
).columns.tolist()
model = build_model(
numeric_features=numeric_features,
categorical_features=categorical_features
)
# Fit WITH survey weights
model.fit(Xtr, y_train, clf__sample_weight=w_train)
cv_mean_recall = model.best_score_
cv_std_recall = model.cv_results_["std_test_score"][model.best_index_]
cv_scores = {k: [] for k in cv_scoring.keys()}
for train_idx, val_idx in cv.split(Xtr, y_train):
X_cv_tr, X_cv_val = Xtr.iloc[train_idx], Xtr.iloc[val_idx]
y_cv_tr, y_cv_val = y_train.iloc[train_idx], y_train.iloc[val_idx]
w_cv_tr = w_train.iloc[train_idx]
est = model.best_estimator_
est.fit(X_cv_tr, y_cv_tr, clf__sample_weight=w_cv_tr)
preds = est.predict(X_cv_val)
probs = est.predict_proba(X_cv_val)[:, 1]
# compute metrics fold-by-fold
fold_metrics = {
"accuracy": accuracy_score(y_cv_val, preds),
"sensitivity": recall_score(y_cv_val, preds),
"specificity": specificity_cv(y_cv_val, preds),
"ppv": precision_score(y_cv_val, preds, zero_division=0),
"npv": npv_cv(y_cv_val, preds),
"f1": fbeta_score(y_cv_val, preds, beta=1),
"f2": fbeta_score(y_cv_val, preds, beta=2),
"auc": roc_auc_score(y_cv_val, probs),
}
for k, v in fold_metrics.items():
cv_scores[k].append(v)
# summarize CV metrics
cv_summary = {}
for metric, values in cv_scores.items():
cv_summary[f"CV_{metric}_mean"] = np.mean(values)
cv_summary[f"CV_{metric}_std"] = np.std(values)
metrics_order = [
"accuracy", "sensitivity", "specificity",
"ppv", "npv", "auc", "f1", "f2"
]
fold_df = pd.DataFrame({
f"Fold {i+1}": [cv_scores[m][i] for m in metrics_order]
for i in range(len(next(iter(cv_scores.values()))))
}, index=metrics_order)
fold_df["Mean"] = fold_df.mean(axis=1)
fold_df["Std Dev"] = fold_df.std(axis=1)
# Add metadata
fold_df.insert(0, "Tier", tier_name)
fold_df.insert(1, "Model", "LogisticRegression")
# Reset index so Metric is a column
fold_df = fold_df.reset_index().rename(columns={"index": "Metric"})
# (optional but recommended) save per-tier fold table here
output_csv = "../results/logreg_5fold_threshold0.50_by_tier.csv"
if os.path.exists(output_csv):
fold_df.to_csv(output_csv, mode="a", header=False, index=False)
else:
fold_df.to_csv(output_csv, index=False)
print(f"Saved fold table to: {output_csv}")
# Probabilities (computed once)
probs = model.predict_proba(Xte)[:, 1]
# ---- Threshold sweep ----
thresholds = np.arange(0.1, 0.51, 0.05)
rows = []
for t in thresholds:
m = evaluate_at_threshold(y_test, probs, threshold=t)
m["Threshold"] = t
rows.append(m)
threshold_df = pd.DataFrame(rows)
'''# Select screening threshold (Recall ≥ 0.90)
screening_candidates = threshold_df[
(threshold_df["Sensitivity"] >= 0.90) &
(threshold_df["NPV"] >= 0.95)
]
screening_row = screening_candidates.sort_values(
by=["Specificity", "F2"],
ascending=[False, False]
).iloc[0]
# Sanity check print
print(
f"Selected screening threshold for {tier_name}: "
f"{screening_row['Threshold']:.2f} "
f"(Sens={screening_row['Sensitivity']:.2f}, "
f"Spec={screening_row['Specificity']:.2f}, "
f"NPV={screening_row['NPV']:.2f})"
)'''
# ---- Bootstrapped CI (AUC is threshold-independent) ----
ci_low, ci_high = bootstrap_auc_ci(model, Xte, y_test)
# ---- Default threshold metrics (0.5) ----
default_metrics = evaluate_at_threshold(y_test, probs, threshold=0.5)
default_metrics.update({
"Tier": tier_name,
"Threshold": 0.5,
"Best_Params": model.best_params_,
**cv_summary,
"AUC_CI_Low": ci_low,
"AUC_CI_High": ci_high,
"Mode": "Default"
})
'''# ---- Screening threshold metrics ----
screening_metrics = screening_row.to_dict()
screening_metrics.update({
"Tier": tier_name,
"Best_Params": model.best_params_,
**cv_summary, # ← ALL CV mean/std metrics
"AUC_CI_Low": ci_low,
"AUC_CI_High": ci_high,
"Mode": "Screening"
})'''
results.append(default_metrics)
# results.append(screening_metrics)
----- Training Logistic Regression on Tier_1_Basic -----
Saved fold table to: ../results/logreg_5fold_threshold0.50_by_tier.csv
----- Training Logistic Regression on Tier_2_Clinical -----
Saved fold table to: ../results/logreg_5fold_threshold0.50_by_tier.csv
----- Training Logistic Regression on Tier_3_Personalized -----
Saved fold table to: ../results/logreg_5fold_threshold0.50_by_tier.csv
Convert to DataFrame + Export
# Convert to DataFrame + Export
results_df = pd.DataFrame(results)
#results_df.to_csv("../results/logreg_tier_performance.csv", index=False)
# 🔍 Sanity check view (THIS LINE)
print(results_df[["Tier", "Mode", "Threshold", "Sensitivity", "NPV", "AUC"]])
Tier Mode Threshold Sensitivity NPV AUC
0 Tier_1_Basic Default 0.5 0.764765 0.943414 0.701264
1 Tier_2_Clinical Default 0.5 0.893894 0.983286 0.922726
2 Tier_3_Personalized Default 0.5 0.905906 0.985301 0.930458