{ "cells": [ { "cell_type": "markdown", "id": "a8e44e75", "metadata": {}, "source": [ "# Random Forest Predictive Model" ] }, { "cell_type": "markdown", "id": "c467f8df", "metadata": {}, "source": [ "Imports" ] }, { "cell_type": "code", "execution_count": 1, "id": "030f5163", "metadata": {}, "outputs": [], "source": [ "import pandas as pd\n", "import numpy as np\n", "import time\n", "import os\n", "from sklearn.model_selection import train_test_split, GridSearchCV, StratifiedKFold, RandomizedSearchCV\n", "from sklearn.preprocessing import StandardScaler\n", "from sklearn.impute import SimpleImputer\n", "from sklearn.pipeline import Pipeline\n", "from sklearn.ensemble import RandomForestClassifier\n", "from sklearn.base import clone\n", "from sklearn.metrics import (\n", " roc_auc_score, recall_score, precision_score, accuracy_score, \n", " confusion_matrix, fbeta_score\n", ")\n", "from sklearn.dummy import DummyClassifier\n", "from sklearn.utils import resample\n", "import matplotlib.pyplot as plt\n", "import os" ] }, { "cell_type": "markdown", "id": "e5ba3077", "metadata": {}, "source": [ "Load tiers" ] }, { "cell_type": "code", "execution_count": 2, "id": "babee804", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Saved: C:\\Users\\agila\\Documents\\STM\\Predicting-Past-Year-Major-Depressive-Episode\\tiers\\tier_1_basic_features.csv\n", "Saved: C:\\Users\\agila\\Documents\\STM\\Predicting-Past-Year-Major-Depressive-Episode\\tiers\\tier_2_clinical_features.csv\n", "Saved: C:\\Users\\agila\\Documents\\STM\\Predicting-Past-Year-Major-Depressive-Episode\\tiers\\tier_3_personalized_features.csv\n" ] } ], "source": [ "from tiering_preparation import TIERS\n", "\n", "tiers = {\n", " \"Tier_1_Basic\": TIERS[\"Tier_1_Basic\"],\n", " \"Tier_2_Clinical\": TIERS[\"Tier_2_Clinical\"],\n", " \"Tier_3_Personalized\": TIERS[\"Tier_3_Personalized\"],\n", "}" ] }, { "cell_type": "code", "execution_count": 3, "id": "5679cc87", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Saving results to: C:\\Users\\agila\\Documents\\STM\\Predicting-Past-Year-Major-Depressive-Episode\\results\n" ] } ], "source": [ "from pathlib import Path\n", "\n", "# 🔒 HARD-ANCHOR PROJECT ROOT (EDIT THIS PATH ONCE)\n", "PROJECT_ROOT = Path(\n", " r\"C:\\Users\\agila\\Documents\\STM\\Predicting-Past-Year-Major-Depressive-Episode\"\n", ")\n", "\n", "RESULTS_DIR = PROJECT_ROOT / \"results\"\n", "FIGURES_DIR = RESULTS_DIR / \"figures\"\n", "\n", "RESULTS_DIR.mkdir(exist_ok=True)\n", "FIGURES_DIR.mkdir(exist_ok=True)\n", "\n", "print(\"Saving results to:\", RESULTS_DIR.resolve())\n" ] }, { "cell_type": "markdown", "id": "0fa8e944", "metadata": {}, "source": [ "Load dataset, drop weight from X, but SAVE w" ] }, { "cell_type": "code", "execution_count": 4, "id": "6d728bf6", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Total samples: 43687\n" ] } ], "source": [ "df = pd.read_pickle(\"../data/prepared/NSDUH_2023_prepared.pkl\")\n", "\n", "TARGET = \"IRAMDEYR\"\n", "WEIGHT = \"ANALWT2_C\"\n", "\n", "# Encode ABS categorical variables\n", "ABS_CATEGORICAL = [\n", " \"IRMARIT_ABS\",\n", " \"IRPINC3_ABS\",\n", " \"WRKDRGHLP_ABS\",\n", " \"KSSLR6MAX_ABS\",\n", " \"IRIMPGOUT_ABS\",\n", " \"WHODASDASC_ABS\",\n", " \"COCLALCUSE_ABS\",\n", "]\n", "\n", "ABS_CATEGORY_MAPS = {}\n", "\n", "for col in ABS_CATEGORICAL:\n", " if col in df.columns:\n", " df[col] = df[col].astype(\"category\")\n", " ABS_CATEGORY_MAPS[col] = dict(enumerate(df[col].cat.categories))\n", " df[col] = df[col].cat.codes\n", "\n", "\n", "# 🔥 REBUILD X AFTER ENCODING\n", "y = df[TARGET]\n", "w = df[WEIGHT]\n", "X = df.drop(columns=[TARGET, WEIGHT])\n", "\n", "print(f'Total samples: {X.shape[0]}')" ] }, { "cell_type": "markdown", "id": "9a84949a", "metadata": {}, "source": [ "Split train/test (consistent with feature selection)" ] }, { "cell_type": "code", "execution_count": 5, "id": "739bca76", "metadata": {}, "outputs": [], "source": [ "X_train, X_test, y_train, y_test, w_train, w_test = train_test_split(\n", " X, y, w,\n", " test_size=0.2,\n", " stratify=y,\n", " random_state=42\n", ")\n" ] }, { "cell_type": "markdown", "id": "60d4bb0e", "metadata": {}, "source": [ "Define metric helper" ] }, { "cell_type": "code", "execution_count": 6, "id": "75f72c77", "metadata": {}, "outputs": [], "source": [ "def evaluate_model(y_true, preds, probs):\n", " tn, fp, fn, tp = confusion_matrix(y_true, preds).ravel()\n", "\n", " specificity = tn / (tn + fp)\n", " sensitivity = recall_score(y_true, preds) # same as recall\n", " ppv = precision_score(y_true, preds)\n", " npv = tn / (tn + fn)\n", "\n", " return {\n", " \"Accuracy\": accuracy_score(y_true, preds),\n", " \"Sensitivity\": sensitivity,\n", " \"Specificity\": specificity,\n", " \"PPV\": ppv,\n", " \"NPV\": npv,\n", " \"AUC\": roc_auc_score(y_true, probs),\n", " \"F1\": fbeta_score(y_true, preds, beta=1),\n", " \"F2\": fbeta_score(y_true, preds, beta=2),\n", " }" ] }, { "cell_type": "markdown", "id": "391c4136", "metadata": {}, "source": [ "Random Forest Pipeline w/ Weighted Fitting + RandomizedSearchCV" ] }, { "cell_type": "code", "execution_count": 7, "id": "966b088c", "metadata": {}, "outputs": [], "source": [ "param_grid = {\n", " \"clf__n_estimators\": [50, 100, 150, 200],\n", " \"clf__max_depth\": [5, 10, 15],\n", " \"clf__min_samples_split\": [2, 5],\n", " \"clf__class_weight\": [None, \"balanced\"],\n", "}\n", "\n", "cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)\n", "\n", "# RandomizedSearchCV\n", "def build_model():\n", " pipe = Pipeline([\n", " (\"clf\", RandomForestClassifier(\n", " random_state=42,\n", " n_jobs=-1\n", " ))\n", " ])\n", "\n", " model = RandomizedSearchCV(\n", " estimator=pipe,\n", " param_distributions=param_grid,\n", " n_iter=20, \n", " scoring=\"recall\",\n", " cv=cv,\n", " n_jobs=-1,\n", " random_state=42,\n", " verbose=2\n", " )\n", "\n", " return model" ] }, { "cell_type": "markdown", "id": "ff87b9f7", "metadata": {}, "source": [ "Bootstrapped Confidence Intervals (ROC + PR)" ] }, { "cell_type": "code", "execution_count": 8, "id": "0517b4e1", "metadata": {}, "outputs": [], "source": [ "def bootstrap_auc_ci(model, X, y, n_boot=1000):\n", " aucs = []\n", "\n", " for _ in range(n_boot):\n", " X_bs, y_bs = resample(X, y, replace=True)\n", " probs = model.predict_proba(X_bs)[:, 1]\n", " aucs.append(roc_auc_score(y_bs, probs))\n", "\n", " return np.percentile(aucs, [2.5, 97.5])" ] }, { "cell_type": "markdown", "id": "d84eeafd", "metadata": {}, "source": [ "Baseline Model (Dummy Classifier)" ] }, { "cell_type": "code", "execution_count": 9, "id": "6a8839f0", "metadata": {}, "outputs": [], "source": [ "# dummy = DummyClassifier(strategy=\"most_frequent\")\n", "# dummy.fit(X_train, y_train)\n", "# dummy_preds = dummy.predict(X_test)\n", "# dummy_probs = np.zeros_like(dummy_preds)\n", "\n", "# dummy_metrics = evaluate_model(y_test, dummy_preds, dummy_probs)" ] }, { "cell_type": "markdown", "id": "7515b795", "metadata": {}, "source": [ "Train + Evaluate Random Forest on Each tier" ] }, { "cell_type": "code", "execution_count": 10, "id": "1cfd0e26", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "----- Training Random Forest on Tier_1_Basic (6 features) -----\n", "Fitting 5 folds for each of 20 candidates, totalling 100 fits\n", "Best params: {'clf__n_estimators': 150, 'clf__min_samples_split': 2, 'clf__max_depth': 5, 'clf__class_weight': 'balanced'}\n", "Completed in 184.7 seconds\n", "\n", "----- Training Random Forest on Tier_2_Clinical (10 features) -----\n", "Fitting 5 folds for each of 20 candidates, totalling 100 fits\n", "Best params: {'clf__n_estimators': 50, 'clf__min_samples_split': 2, 'clf__max_depth': 5, 'clf__class_weight': 'balanced'}\n", "Completed in 171.8 seconds\n", "\n", "----- Training Random Forest on Tier_3_Personalized (19 features) -----\n", "Fitting 5 folds for each of 20 candidates, totalling 100 fits\n", "Best params: {'clf__n_estimators': 100, 'clf__min_samples_split': 2, 'clf__max_depth': 5, 'clf__class_weight': 'balanced'}\n", "Completed in 179.4 seconds\n" ] } ], "source": [ "results = []\n", "best_models = {}\n", "\n", "for tier_name, features in tiers.items():\n", "\n", " print(f\"\\n----- Training Random Forest on {tier_name} ({len(features)} features) -----\")\n", " \n", " start_time = time.time()\n", "\n", " Xtr = X_train[features]\n", " Xte = X_test[features]\n", "\n", " model = build_model()\n", "\n", " bad_cols = Xtr.select_dtypes(include=[\"object\", \"string\", \"category\"]).columns.tolist()\n", " if bad_cols:\n", " raise ValueError(\n", " f\"[{tier_name}] Non-numeric columns found: {bad_cols}\"\n", " )\n", "\n", " # Fit WITH survey weights\n", " model.fit(Xtr, y_train, clf__sample_weight=w_train)\n", "\n", " # Save best model for fold-wise CV later\n", " best_models[tier_name] = model.best_estimator_\n", "\n", " preds = model.predict(Xte)\n", " probs = model.predict_proba(Xte)[:, 1]\n", "\n", " metrics = evaluate_model(y_test, preds, probs)\n", "\n", " # Bootstrapped CI\n", " ci_low, ci_high = bootstrap_auc_ci(model, Xte, y_test)\n", "\n", " metrics[\"tier\"] = tier_name\n", " metrics[\"threshold\"] = 0.50\n", " metrics[\"n_features\"] = len(features)\n", " metrics[\"AUC_CI_Low\"] = ci_low\n", " metrics[\"AUC_CI_High\"] = ci_high\n", " metrics[\"Best_Params\"] = str(model.best_params_)\n", "\n", " results.append(metrics)\n", " \n", " print(f\"Best params: {model.best_params_}\")\n", " print(f\"Completed in {time.time() - start_time:.1f} seconds\")" ] }, { "cell_type": "code", "execution_count": 11, "id": "7b8db549", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "===== Tier_1_Basic =====\n", "Fold Metric Tier Fold 1 Fold 2 Fold 3 Fold 4 \\\n", "0 Accuracy Tier_1_Basic 0.625036 0.661087 0.675536 0.655365 \n", "1 Sensitivity Tier_1_Basic 0.628285 0.636250 0.620000 0.622500 \n", "2 Specificity Tier_1_Basic 0.624616 0.664297 0.682714 0.659612 \n", "3 PPV Tier_1_Basic 0.177636 0.196753 0.201626 0.191171 \n", "4 NPV Tier_1_Basic 0.928674 0.933909 0.932892 0.931129 \n", "5 AUC Tier_1_Basic 0.692042 0.711625 0.718023 0.701979 \n", "6 F1 Tier_1_Basic 0.276966 0.300561 0.304294 0.292511 \n", "7 F2 Tier_1_Basic 0.416805 0.439779 0.438163 0.428941 \n", "\n", "Fold Fold 5 Mean Std Dev \n", "0 0.682787 0.659962 0.020029 \n", "1 0.613267 0.624060 0.007770 \n", "2 0.691761 0.664600 0.023193 \n", "3 0.204337 0.194305 0.009465 \n", "4 0.932694 0.931860 0.001824 \n", "5 0.715546 0.707843 0.009606 \n", "6 0.306537 0.296174 0.010724 \n", "7 0.437969 0.432331 0.008647 \n", "\n", "===== Tier_2_Clinical =====\n", "Fold Metric Tier Fold 1 Fold 2 Fold 3 Fold 4 \\\n", "0 Accuracy Tier_2_Clinical 0.833619 0.826753 0.834621 0.821888 \n", "1 Sensitivity Tier_2_Clinical 0.901126 0.893750 0.898750 0.902500 \n", "2 Specificity Tier_2_Clinical 0.824907 0.818094 0.826333 0.811470 \n", "3 PPV Tier_2_Clinical 0.399113 0.388376 0.400780 0.382213 \n", "4 NPV Tier_2_Clinical 0.984767 0.983492 0.984411 0.984709 \n", "5 AUC Tier_2_Clinical 0.931201 0.917699 0.928792 0.924019 \n", "6 F1 Tier_2_Clinical 0.553208 0.541462 0.554356 0.537003 \n", "7 F2 Tier_2_Clinical 0.720000 0.709185 0.719864 0.709373 \n", "\n", "Fold Fold 5 Mean Std Dev \n", "0 0.840607 0.831498 0.006513 \n", "1 0.908636 0.900952 0.004862 \n", "2 0.831826 0.822526 0.007050 \n", "3 0.410866 0.396270 0.010014 \n", "4 0.986021 0.984680 0.000811 \n", "5 0.931139 0.926570 0.005147 \n", "6 0.565861 0.550378 0.010217 \n", "7 0.731412 0.717967 0.008240 \n", "\n", "===== Tier_3_Personalized =====\n", "Fold Metric Tier Fold 1 Fold 2 Fold 3 \\\n", "0 Accuracy Tier_3_Personalized 0.841917 0.836624 0.833047 \n", "1 Sensitivity Tier_3_Personalized 0.904881 0.866250 0.913750 \n", "2 Specificity Tier_3_Personalized 0.833791 0.832795 0.822617 \n", "3 PPV Tier_3_Personalized 0.412671 0.401042 0.399672 \n", "4 NPV Tier_3_Personalized 0.985491 0.979666 0.986630 \n", "5 AUC Tier_3_Personalized 0.934691 0.922078 0.931718 \n", "6 F1 Tier_3_Personalized 0.566837 0.548259 0.556105 \n", "7 F2 Tier_3_Personalized 0.730598 0.703125 0.726785 \n", "\n", "Fold Fold 4 Fold 5 Mean Std Dev \n", "0 0.834907 0.845185 0.838336 0.004525 \n", "1 0.903750 0.906133 0.898953 0.016722 \n", "2 0.826010 0.837318 0.830506 0.005383 \n", "3 0.401667 0.418255 0.406661 0.007429 \n", "4 0.985164 0.985736 0.984537 0.002484 \n", "5 0.932145 0.939194 0.931965 0.005613 \n", "6 0.556154 0.572332 0.559937 0.008563 \n", "7 0.723000 0.734727 0.723647 0.010977 \n", "\n", "Saved: rf_5fold_threshold0.50_by_tier.csv\n" ] } ], "source": [ "# 5-Fold Cross-Validation on TRAIN set (no test leakage)\n", "skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)\n", "metrics_order = [\"Accuracy\", \"Sensitivity\", \"Specificity\", \"PPV\", \"NPV\", \"AUC\", \"F1\", \"F2\"]\n", "\n", "all_fold_dfs = []\n", "\n", "for tier_name, features in tiers.items():\n", " print(f\"\\n===== {tier_name} =====\")\n", " \n", " fold_rows = []\n", " best_model = best_models[tier_name]\n", " \n", " X_tier = X_train[features]\n", " w_tier = w_train\n", " \n", " for fold_i, (tr_idx, val_idx) in enumerate(skf.split(X_tier, y_train), start=1):\n", " Xtr_fold, Xval_fold = X_tier.iloc[tr_idx], X_tier.iloc[val_idx]\n", " ytr_fold, yval_fold = y_train.iloc[tr_idx], y_train.iloc[val_idx]\n", " wtr_fold = w_tier.iloc[tr_idx]\n", " \n", " model_clone = clone(best_model)\n", "\n", " assert not X_tier.select_dtypes(include=[\"object\", \"category\"]).any().any(), \\\n", " f\"{tier_name} still contains categorical dtypes\"\n", "\n", " if hasattr(model_clone, \"named_steps\"):\n", " model_clone.fit(Xtr_fold, ytr_fold, clf__sample_weight=wtr_fold)\n", " else:\n", " model_clone.fit(Xtr_fold, ytr_fold, sample_weight=wtr_fold)\n", " \n", " preds = model_clone.predict(Xval_fold)\n", " probs = model_clone.predict_proba(Xval_fold)[:, 1]\n", " \n", " fold_metrics = evaluate_model(yval_fold, preds, probs)\n", " fold_metrics[\"Fold\"] = f\"Fold {fold_i}\"\n", " fold_rows.append(fold_metrics)\n", " \n", " # Table: rows=Metric, cols=Fold1..Fold5 (TRANSPOSE like LightGBM)\n", " fold_df = pd.DataFrame(fold_rows).set_index(\"Fold\")[metrics_order].T\n", " fold_df[\"Mean\"] = fold_df.mean(axis=1)\n", " fold_df[\"Std Dev\"] = fold_df.std(axis=1)\n", " fold_df.insert(0, \"Tier\", tier_name)\n", " fold_df = fold_df.reset_index().rename(columns={\"index\": \"Metric\"})\n", " \n", " all_fold_dfs.append(fold_df)\n", " print(fold_df)\n", "\n", "# Combine and export\n", "cv_results_df = pd.concat(all_fold_dfs, ignore_index=True)\n", "cv_results_df.to_csv(\n", " RESULTS_DIR / \"rf_5fold_threshold0.50_by_tier.csv\",\n", " index=False\n", ")\n", "print(\"\\nSaved: rf_5fold_threshold0.50_by_tier.csv\")" ] }, { "cell_type": "markdown", "id": "83b9aad4", "metadata": {}, "source": [ "Convert to DataFrame + Export" ] }, { "cell_type": "code", "execution_count": 12, "id": "46eefea0", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " tier threshold n_features Accuracy Sensitivity \\\n", "0 Tier_1_Basic 0.5 6 0.682879 0.602603 \n", "1 Tier_2_Clinical 0.5 10 0.835203 0.874875 \n", "2 Tier_3_Personalized 0.5 19 0.834173 0.885886 \n", "\n", " Specificity PPV NPV AUC F1 F2 AUC_CI_Low \\\n", "0 0.693242 0.202285 0.931100 0.711132 0.302893 0.431727 0.694228 \n", "1 0.830081 0.399269 0.980913 0.923757 0.548306 0.706548 0.915998 \n", "2 0.827497 0.398649 0.982510 0.931227 0.549860 0.711873 0.924279 \n", "\n", " AUC_CI_High Best_Params \n", "0 0.726792 {'clf__n_estimators': 150, 'clf__min_samples_s... \n", "1 0.930662 {'clf__n_estimators': 50, 'clf__min_samples_sp... \n", "2 0.938016 {'clf__n_estimators': 100, 'clf__min_samples_s... \n" ] } ], "source": [ "results_df = pd.DataFrame(results)\n", "\n", "# Reorder columns to match LightGBM format\n", "front_cols = [\"tier\", \"threshold\", \"n_features\"]\n", "other_cols = [c for c in results_df.columns if c not in front_cols]\n", "results_df = results_df[front_cols + other_cols]\n", "\n", "results_df.to_csv(\n", " RESULTS_DIR / \"rf_two_rows_per_tier.csv\",\n", " index=False\n", ")\n", "\n", "print(results_df)" ] }, { "cell_type": "code", "execution_count": null, "id": "b43bf453", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "venv", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.13.9" } }, "nbformat": 4, "nbformat_minor": 5 }