4_Feature_Analyses.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import random\n",
    "import scipy.stats\n",
    "import pickle, os, time\n",
    "import itertools\n",
    "from datetime import datetime, timedelta\n",
    "from collections import Counter, defaultdict, namedtuple\n",
    "from PIL import Image\n",
    "import yaml\n",
    "from tqdm import tqdm\n",
    "\n",
    "from sklearn import preprocessing, model_selection, metrics, utils\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from tqdm import tqdm\n",
    "from joblib import Parallel, delayed\n",
    "from sklearn.base import clone\n",
    "\n",
    "import seaborn as sns\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "data_dir = '/data/GVHD/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%config InlineBackend.figure_format = 'svg'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Helper functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model training"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def train_model(Xtr, ytr):\n",
    "    np.random.seed(42)\n",
    "    random.seed(42)\n",
    "    \n",
    "    # Specify hyperparameters and cv parameters\n",
    "    base_estimator = LogisticRegression(\n",
    "        penalty='l2', \n",
    "        class_weight='balanced', \n",
    "        solver='liblinear'\n",
    "    )\n",
    "    param_grid = {\n",
    "        'C': [10. ** n for n in range(-6, 7)],\n",
    "        'penalty': ['l2'],\n",
    "    }\n",
    "    \n",
    "    cv_splits, cv_repeat = 5, 20\n",
    "    cv = model_selection.RepeatedStratifiedKFold(cv_splits, cv_repeat, random_state=0)\n",
    "    clf = model_selection.GridSearchCV(\n",
    "        clone(base_estimator), param_grid, \n",
    "        cv=cv, scoring='roc_auc', iid=False, n_jobs=5,\n",
    "    )\n",
    "    clf.fit(Xtr, ytr)\n",
    "    return clf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Evaluation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def boostrap_func_all(i, y_true, y_prob, threshold):\n",
    "    fpr, tpr, thresholds = metrics.roc_curve(y_true, y_prob)\n",
    "    y_true_b, y_prob_b = utils.resample(y_true, y_prob, replace=True, random_state=i)\n",
    "    y_pred_b = (y_prob_b > threshold)\n",
    "    tpr_cutoff = metrics.recall_score(y_true_b, y_pred_b)\n",
    "    idx = (np.abs(tpr - tpr_cutoff)).argmin()\n",
    "    \n",
    "    return (\n",
    "        metrics.roc_auc_score(y_true_b, y_prob_b), # AUC\n",
    "        tpr[idx], # sensitivity\n",
    "        1-fpr[idx], # specificity\n",
    "        metrics.precision_score(y_true_b, y_pred_b), # positive predictive value\n",
    "    )\n",
    "\n",
    "def boostrap_func_confusion(i, y_true, y_prob, threshold):\n",
    "    fpr, tpr, thresholds = metrics.roc_curve(y_true, y_prob)\n",
    "    y_true_b, y_prob_b = utils.resample(y_true, y_prob, replace=True, random_state=i)\n",
    "    y_pred_b = (y_prob_b > threshold)\n",
    "    tpr_cutoff = metrics.recall_score(y_true_b, y_pred_b)\n",
    "    idx = (np.abs(tpr - tpr_cutoff)).argmin()\n",
    "    \n",
    "    return metrics.confusion_matrix(y_true_b, y_pred_b).ravel()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def evaluate_model(clf, Xte, yte, Xtr=None, threshold_p=55, verbose=True):\n",
    "    y_true = yte\n",
    "    y_score = clf.decision_function(Xte)\n",
    "    y_prob = clf.predict_proba(Xte)[:,1]\n",
    "    \n",
    "    fpr, tpr, thresholds = metrics.roc_curve(y_true, y_score)\n",
    "    test_auc = metrics.roc_auc_score(y_true, y_score)\n",
    "    \n",
    "    # Picking a risk threshold based on training set if possible\n",
    "    if Xtr is not None:\n",
    "        if verbose: print('Risk threshold based on train set')\n",
    "        threshold = np.percentile(clf.predict_proba(Xtr)[:,1], threshold_p)\n",
    "    else:\n",
    "        if verbose: print('Risk threshold based on test set')\n",
    "        threshold = np.percentile(y_score, 55)\n",
    "    if verbose: print('p_Threshold', threshold)\n",
    "    \n",
    "    if verbose: print()\n",
    "    if verbose: print('Confusion matrix (95%CI lower, upper)')\n",
    "    y_pred = (y_prob > threshold)\n",
    "    conf_mat = metrics.confusion_matrix(y_true, y_pred)\n",
    "    if verbose: print(conf_mat) # Rows: actual 0, actual 1; Cols: predicted 0, predicted 1\n",
    "    \n",
    "    confmats = [boostrap_func_confusion(i, y_true, y_prob, threshold) for i in range(1000)]\n",
    "    confmats_ = np.asarray(confmats)\n",
    "    if verbose: print(np.percentile(confmats_, 2.5, axis=0).reshape(2,2))\n",
    "    if verbose: print(np.percentile(confmats_, 97.5, axis=0).reshape(2,2))\n",
    "\n",
    "    if verbose: print()\n",
    "    if verbose: print('scores')\n",
    "    tpr_ = metrics.recall_score(y_true, y_pred)\n",
    "    idx = (np.abs(tpr - tpr_)).argmin()\n",
    "    if verbose: print('AUROC={:.3f}'.format(test_auc))\n",
    "    if verbose: print('TPR={:.3f}, FPR={:.3f} PPV={:.3f}'.format(\n",
    "        tpr[idx], fpr[idx],\n",
    "        metrics.precision_score(y_true, y_pred)))\n",
    "\n",
    "    if verbose: print()\n",
    "    if verbose: print('scores (95%CI lower, upper)')\n",
    "    auc_scores, sensitivities, specificities, ppvs = zip(*[boostrap_func_all(i, y_true, y_prob, threshold) for i in range(1000)])\n",
    "    if verbose: print('Test AUC {:.3f} ({:.3f}, {:.3f})'.format(np.median(auc_scores), np.percentile(auc_scores, 2.5), np.percentile(auc_scores, 97.5)))\n",
    "    if verbose: print('Test AUC {:.3f} ± {:.3f}'.format(np.mean(auc_scores), np.std(auc_scores)))\n",
    "    if verbose: print('sens.\\t {:.1%} ({:.1%}, {:.1%})'.format(np.mean(sensitivities), np.percentile(sensitivities, 2.5), np.percentile(sensitivities, 97.5)))\n",
    "    if verbose: print('spec.\\t {:.1%} ({:.1%}, {:.1%})'.format(np.mean(specificities), np.percentile(specificities, 2.5), np.percentile(specificities, 97.5)))\n",
    "    if verbose: print('prec.\\t {:.1%} ({:.1%}, {:.1%})'.format(np.mean(ppvs), np.percentile(ppvs, 2.5), np.percentile(ppvs, 97.5)))\n",
    "    \n",
    "    return test_auc, auc_scores, fpr, tpr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data loading"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "pop = pd.read_csv(data_dir + 'population/d10_with_vitals.csv').set_index('BMT_ID')\n",
    "extracted_features = pd.read_csv('data/ts_features.csv', index_col='id')\n",
    "df_label = pop.join(pd.read_csv(data_dir + 'prep/label.csv', index_col='BMT_ID'), how='left')\n",
    "df_features = pd.read_csv('data/df_features.csv', index_col='id')\n",
    "feature_names = df_features.columns"
   ]
  },
  {
   "cell_type": "code",