scripts/07_calculate_metrics.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "aswNeob-p4vY"
   },
   "source": [
    "# Master Script 7: Calculate bootstrapping bias-corrected cross-validation (BBC-CV) area under the receiver operating characteristic curve (AUC) and classification metrics\n",
    "\n",
    "Shubhayu Bhattacharyay\n",
    "<br>\n",
    "University of Cambridge\n",
    "<br>\n",
    "Johns Hopkins University\n",
    "<br>\n",
    "email address: sb2406@cam.ac.uk"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Contents:\n",
    "### I. Initialization\n",
    "### II. Define functions for calculating metrics and ROC curve with BBC-CV\n",
    "### III. Calculate metrics for GCSm threshold-level detection\n",
    "### IV. Calculate metrics for GOSE (at discharge) threshold-level prediction\n",
    "### V. Calculate metrics for GOSE (at 12 months) threshold-level prediction\n",
    "### VI. Calculate precision-recall information for GOSE > 5 (at discharge) prediction with 6-hour observation window"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## I. Initialization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import necessary packages and methods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "executionInfo": {
     "elapsed": 470,
     "status": "ok",
     "timestamp": 1619452672623,
     "user": {
      "displayName": "Shubhayu Bhattacharyay",
      "photoUrl": "",
      "userId": "16704285935129215396"
     },
     "user_tz": -60
    },
    "id": "8D3-LDzSp4vf"
   },
   "outputs": [],
   "source": [
    "# Fundamental methods\n",
    "import os\n",
    "import re\n",
    "import sys\n",
    "import time\n",
    "import glob\n",
    "import random\n",
    "import warnings\n",
    "import itertools\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import pickle as cp\n",
    "import seaborn as sns\n",
    "import multiprocessing\n",
    "from scipy import interp\n",
    "from pathlib import Path\n",
    "import matplotlib.pyplot as plt\n",
    "from IPython.display import clear_output\n",
    "\n",
    "# tqdm method for progress monitoring\n",
    "from tqdm import tqdm\n",
    "\n",
    "# SciKit-Learn methods\n",
    "from sklearn.metrics import classification_report, roc_auc_score, roc_curve, auc, mean_squared_error, accuracy_score, confusion_matrix, precision_recall_curve, average_precision_score \n",
    "from sklearn.calibration import calibration_curve"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "JAJD9Ib5p4vg"
   },
   "source": [
    "### Establish necessary parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Establish number of cores for parallel processing (default: save 2 cores)\n",
    "NUM_CORES = multiprocessing.cpu_count() - 2\n",
    "\n",
    "# Establish number of resamples for bias-corrected bootstrapping of confidence intervals\n",
    "NUM_BOOTSTRAPS = 1000\n",
    "\n",
    "# Ignore potential warnings\n",
    "warnings.filterwarnings(action='ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## II. Define functions for calculating metrics and ROC curve with BBC-CV"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Function for calculating AUC and ROC axes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_bs_AUC(\n",
    "    curr_compiled_predictions,\n",
    "    curr_in_sample_UPIs,\n",
    "    metric='AUC',\n",
    "    roc_axes=True,\n",
    "    out_interp_fpr=np.linspace(0, 1, num=200),\n",
    "    ):\n",
    "\n",
    "    # Derive list of unique UPIs based on compiled prediction dataframe\n",
    "\n",
    "    uniq_UPIs = curr_compiled_predictions.UPI.unique()\n",
    "\n",
    "    # Get current out samples based on current bootstrap in-sample\n",
    "\n",
    "    curr_out_sample_UPIs = np.setdiff1d(uniq_UPIs, curr_in_sample_UPIs)\n",
    "\n",
    "    # Separate predictions into in-sample and out-sample based on current bootstrap resample\n",
    "\n",
    "    curr_in_sample_predictions = \\\n",
    "        curr_compiled_predictions[curr_compiled_predictions.UPI.isin(curr_in_sample_UPIs)]\n",
    "    curr_out_sample_predictions = \\\n",
    "        curr_compiled_predictions[curr_compiled_predictions.UPI.isin(curr_out_sample_UPIs)]\n",
    "\n",
    "    # Calculate in-sample AUCs to determine optimal model configuration\n",
    "\n",
    "    grouped_in_sample_predictions = \\\n",
    "        curr_in_sample_predictions.groupby('ConfigIdx', sort=False)\n",
    "    in_sample_AUCs = pd.DataFrame(np.empty((0, 2)), columns=['ConfigIdx'\n",
    "                                  , 'AUC'])\n",
    "\n",
    "    # Iterate through predictions group by configuration index and calculate group-specific AUC\n",
    "\n",
    "    for (name, group) in grouped_in_sample_predictions:\n",
    "        in_sample_AUCs = \\\n",
    "            in_sample_AUCs.append(pd.DataFrame({'ConfigIdx': name,\n",
    "                                  'AUC': roc_auc_score(group.TrueLabel,\n",
    "                                  group.Prob)}, index=[0]),\n",
    "                                  ignore_index=True)\n",
    "\n",
    "    # Fix datatype of `ConfigIdx`\n",
    "\n",
    "    in_sample_AUCs.ConfigIdx = in_sample_AUCs.ConfigIdx.astype('int')\n",
    "\n",
    "    # Isolate optimal configuration index for current sampling\n",
    "\n",
    "    curr_opt_ConfigIdx = in_sample_AUCs.nlargest(1, 'AUC', keep='first'\n",
    "            ).ConfigIdx.values[0]\n",
    "\n",
    "    # Filter out out-sample predictions of optimal configuration index\n",
    "\n",
    "    curr_out_sample_predictions = \\\n",
    "        curr_out_sample_predictions[curr_out_sample_predictions.ConfigIdx\n",
    "                                    == curr_opt_ConfigIdx]\n",
    "\n",
    "    # Calculate out-sample AUC\n",
    "\n",
    "    out_sample_AUC = \\\n",
    "        roc_auc_score(curr_out_sample_predictions.TrueLabel,\n",
    "                      curr_out_sample_predictions.Prob)\n",
    "\n",
    "    # Calculate out-sample ROC curve axes and interpolate to common FPR axis\n",
    "\n",
    "    if roc_axes:\n",
    "        (out_sample_fpr, out_sample_tpr, _) = \\\n",
    "            roc_curve(curr_out_sample_predictions.TrueLabel,\n",
    "                      curr_out_sample_predictions.Prob)\n",
    "        out_interp_tpr = np.interp(out_interp_fpr, out_sample_fpr,\n",
    "                                   out_sample_tpr)\n",
    "\n",
    "        # Return AUC, optimal configuration, and ROC curve axes\n",
    "\n",
    "        return (out_sample_AUC, curr_opt_ConfigIdx, out_interp_fpr,\n",
    "                out_interp_tpr)\n",
    "    else:\n",
    "\n",
    "        # Return just AUC and optimal configuration if ROC not requested\n",
    "\n",
    "        return (out_sample_AUC, curr_opt_ConfigIdx)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Function for calculating Precision Recall axes and associated AUC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_bs_PRC(\n",
    "    curr_compiled_predictions,\n",
    "    curr_in_sample_UPIs,\n",
    "    metric='AUPRC',\n",
    "    prc_axes=True,\n",
    "    out_interp_rec=np.linspace(0, 1, num=1000),\n",
    "    ):\n",
    "\n",
    "    # Derive list of unique UPIs based on compiled prediction dataframe\n",
    "\n",
    "    uniq_UPIs = curr_compiled_predictions.UPI.unique()\n",
    "\n",
    "    # Get current out samples based on current bootstrap in-sample\n",
    "\n",
    "    curr_out_sample_UPIs = np.setdiff1d(uniq_UPIs, curr_in_sample_UPIs)\n",
    "\n",
    "    # Separate predictions into in-sample and out-sample based on current bootstrap resample\n",
    "\n",
    "    curr_in_sample_predictions = \\\n",
    "        curr_compiled_predictions[curr_compiled_predictions.UPI.isin(curr_in_sample_UPIs)]\n",
    "    curr_out_sample_predictions = \\\n",
    "        curr_compiled_predictions[curr_compiled_predictions.UPI.isin(curr_out_sample_UPIs)]\n",
    "\n",
    "    # Calculate in-sample AUPRCs to determine optimal model configuration\n",
    "\n",
    "    grouped_in_sample_predictions = \\\n",
    "        curr_in_sample_predictions.groupby('ConfigIdx', sort=False)\n",
    "    in_sample_AUPRCs = pd.DataFrame(np.empty((0, 2)), columns=['ConfigIdx'\n",
    "                                  , 'AUPRC'])\n",
    "\n",
    "    # Iterate through predictions group by configuration index and calculate group-specific AUPRC\n",
    "\n",
    "    for (name, group) in grouped_in_sample_predictions:\n",
    "        \n",
    "        #(in_sample_prec, in_sample_rec, _) = precision_recall_curve(group.TrueLabel,group.Prob)\n",
    "        \n",
    "        in_sample_AUPRCs = \\\n",
    "            in_sample_AUPRCs.append(pd.DataFrame({'ConfigIdx': name,\n",
    "                                  'AUPRC': average_precision_score(group.TrueLabel,\n",
    "                                  group.Prob)}, index=[0]),\n",
    "                                  ignore_index=True)\n",
    "\n",
    "    # Fix datatype of `ConfigIdx`\n",
    "\n",
    "    in_sample_AUPRCs.ConfigIdx = in_sample_AUPRCs.ConfigIdx.astype('int')\n",
    "\n",
    "    # Isolate optimal configuration index for current sampling\n",
    "\n",
    "    curr_opt_ConfigIdx = in_sample_AUPRCs.nlargest(1, 'AUPRC', keep='first'\n",
    "            ).ConfigIdx.values[0]\n",
    "\n",
    "    # Filter out out-sample predictions of optimal configuration index\n",
    "\n",
    "    curr_out_sample_predictions = \\\n",
    "        curr_out_sample_predictions[curr_out_sample_predictions.ConfigIdx\n",
    "                                    == curr_opt_ConfigIdx]\n",
    "\n",
    "    # Calculate out-sample AUPRC\n",
    "\n",
    "    (out_sample_prec, out_sample_rec, _) = precision_recall_curve(curr_out_sample_predictions.TrueLabel,curr_out_sample_predictions.Prob)\n",
    "    \n",
    "    out_sample_AUPRC = \\\n",
    "        average_precision_score(curr_out_sample_predictions.TrueLabel,curr_out_sample_predictions.Prob)\n",
    "\n",
    "    # Calculate out-sample Precision-Recall curve axes and interpolate to common Recall axis\n",
    "\n",
    "    if prc_axes:\n",
    "        \n",
    "        out_interp_prec = np.interp(out_interp_rec,np.flip(out_sample_rec),np.flip(out_sample_prec))\n",
    "\n",
    "        # Return AUPRC, optimal configuration, and Precision-Recall curve axes\n",
    "\n",
    "        return (out_sample_AUPRC, curr_opt_ConfigIdx, out_interp_rec,\n",
    "                out_interp_prec)\n",
    "    else:\n",
    "\n",
    "        # Return just AUPRC and optimal configuration if Precision-Recall not requested\n",
    "\n",
    "        return (out_sample_AUPRC, curr_opt_ConfigIdx)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Function for calculating classification metrics\n",
    "Possible metrics: `precision`, `recall` (i.e., sensitivity), `f1-score`, and `specificity`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_bs_classification(curr_compiled_predictions,\n",
    "                                curr_in_sample_UPIs, metric='precision'\n",
    "                                ):\n",
    "\n",
    "    # Derive list of unique UPIs based on compiled prediction dataframe\n",
    "\n",
    "    uniq_UPIs = curr_compiled_predictions.UPI.unique()\n",
    "\n",
    "    # Get current out samples based on current bootstrap in-sample\n",
    "\n",
    "    curr_out_sample_UPIs = np.setdiff1d(uniq_UPIs, curr_in_sample_UPIs)\n",
    "\n",
    "    # Separate predictions into in-sample and out-sample based on current bootstrap resample\n",
    "\n",
    "    curr_in_sample_predictions = \\\n",
    "        curr_compiled_predictions[curr_compiled_predictions.UPI.isin(curr_in_sample_UPIs)]\n",
    "    curr_out_sample_predictions = \\\n",
    "        curr_compiled_predictions[curr_compiled_predictions.UPI.isin(curr_out_sample_UPIs)]\n",
    "\n",
    "    # Calculate in-sample metrics to determine optimal model configuration\n",
    "\n",
    "    grouped_in_sample_predictions = \\\n",
    "        curr_in_sample_predictions.groupby('ConfigIdx', sort=False)\n",
    "    in_sample_classifications = pd.DataFrame(np.empty((0, 2)),\n",
    "            columns=['ConfigIdx', 'Metric'])\n",
    "\n",
    "    # Iterate through predictions group by configuration index and calculate group-specific metric\n",
    "\n",
    "    for (name, group) in grouped_in_sample_predictions:\n",
    "        if metric == 'specificity':\n",
    "            curr_group_metric_value = \\\n",
    "                classification_report(group.TrueLabel, (group.Prob\n",
    "                    >= .5).astype('int').values, output_dict=True)['0'\n",
    "                    ]['recall']\n",
    "        else:\n",
    "            curr_group_metric_value = \\\n",
    "                classification_report(group.TrueLabel, (group.Prob\n",
    "                    >= .5).astype('int').values, output_dict=True)['1'\n",
    "                    ][metric]\n",
    "        in_sample_classifications = \\\n",
    "            in_sample_classifications.append(pd.DataFrame({'ConfigIdx': name,\n",
    "                'Metric': curr_group_metric_value}, index=[0]),\n",
    "                ignore_index=True)\n",
    "\n",
    "    # Fix datatype of `ConfigIdx`\n",
    "\n",
    "    in_sample_classifications.ConfigIdx = \\\n",
    "        in_sample_classifications.ConfigIdx.astype('int')\n",
    "\n",
    "    # Isolate optimal configuration index for current sampling\n",
    "\n",
    "    curr_opt_ConfigIdx = in_sample_classifications.nlargest(1, 'Metric'\n",
    "            , keep='first').ConfigIdx.values[0]\n",
    "\n",
    "    # Filter out out-sample predictions of optimal configuration index\n",
    "\n",
    "    curr_out_sample_predictions = \\\n",
    "        curr_out_sample_predictions[curr_out_sample_predictions.ConfigIdx\n",
    "                                    == curr_opt_ConfigIdx]\n",
    "\n",
    "    # Calculate out-sample metric\n",
    "\n",
    "    if metric == 'specificity':\n",
    "        out_sample_classifications = \\\n",
    "            classification_report(curr_out_sample_predictions.TrueLabel,\n",
    "                                  (curr_out_sample_predictions.Prob\n",
    "                                  >= .5).astype('int').values,\n",
    "                                  output_dict=True)['0']['recall']\n",
    "    else:\n",
    "        out_sample_classifications = \\\n",
    "            classification_report(curr_out_sample_predictions.TrueLabel,\n",
    "                                  (curr_out_sample_predictions.Prob\n",
    "                                  >= .5).astype('int').values,\n",
    "                                  output_dict=True)['1'][metric]\n",
    "\n",
    "    # Return out-sample metric value and optimal configuration\n",
    "\n",
    "    return (out_sample_classifications, curr_opt_ConfigIdx)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sub-function for each parallelization core"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def _calculate_bs_metrics(\n",
    "    metric_fun,\n",
    "    curr_compiled_predictions,\n",
    "    size,\n",
    "    metric,\n",
    "    progress_bar=True,\n",
    "    progress_bar_desc='',\n",
    "    ):\n",
    "\n",
    "    # Create random generator instance\n",
    "\n",
    "    rg = np.random.default_rng()\n",
    "\n",
    "    # Get list of unique UPIs\n",
    "\n",
    "    uniq_UPIs = curr_compiled_predictions.UPI.unique()\n",
    "\n",
    "    # Derive current bootstrapping sample\n",
    "\n",
    "    bootstrap_resamples = rg.choice(uniq_UPIs, (len(uniq_UPIs), size),\n",
    "                                    replace=True)\n",
    "\n",
    "    # Go through current resamples and ensure positive and negative cases exist for each resample\n",
    "\n",
    "    for bs_col in range(bootstrap_resamples.shape[1]):\n",
    "\n",
    "        # Get current in and out samples based on bootstrap resample\n",
    "\n",
    "        curr_in_sample_UPIs = np.unique(bootstrap_resamples[:, bs_col])\n",
    "        curr_out_sample_UPIs = np.setdiff1d(uniq_UPIs,\n",
    "                curr_in_sample_UPIs)\n",
    "\n",
    "        # Split current compiled prediction dataframe into in-sample and out-sample predictions\n",
    "\n",
    "        curr_in_sample_predictions = \\\n",
    "            curr_compiled_predictions[curr_compiled_predictions.UPI.isin(curr_in_sample_UPIs)]\n",
    "        curr_out_sample_predictions = \\\n",
    "            curr_compiled_predictions[curr_compiled_predictions.UPI.isin(curr_out_sample_UPIs)]\n",
    "\n",
    "        # If either the in or out sample contains only one outcome, resample current bootstrap until there is\n",
    "\n",
    "        if (len(curr_in_sample_predictions.TrueLabel.unique()) == 1) \\\n",
    "            | (len(curr_out_sample_predictions.TrueLabel.unique())\n",
    "               == 1):\n",
    "\n",
    "            fail_condition = True\n",
    "\n",
    "            while fail_condition:\n",
    "                temp_resample = rg.choice(uniq_UPIs, (len(uniq_UPIs),\n",
    "                        1), replace=True)\n",
    "                temp_in_sample_UPIs = np.unique(temp_resample[:, 0])\n",
    "                temp_out_sample_UPIs = np.setdiff1d(uniq_UPIs,\n",
    "                        temp_in_sample_UPIs)\n",
    "                temp_in_sample_predictions = \\\n",
    "                    curr_compiled_predictions[curr_compiled_predictions.UPI.isin(temp_in_sample_UPIs)]\n",
    "                temp_out_sample_predictions = \\\n",
    "                    curr_compiled_predictions[curr_compiled_predictions.UPI.isin(temp_out_sample_UPIs)]\n",
    "\n",
    "                if (len(temp_in_sample_predictions.TrueLabel.unique())\n",
    "                    == 2) \\\n",
    "                    & (len(temp_out_sample_predictions.TrueLabel.unique())\n",
    "                       == 2):\n",
    "                    bootstrap_resamples[:, bs_col] = temp_resample[:, 0]\n",
    "                    fail_condition = False\n",
    "\n",
    "    if progress_bar:\n",
    "        iterator = tqdm(range(size), desc=progress_bar_desc)\n",
    "    else:\n",
    "        iterator = range(size)\n",
    "\n",
    "    return [metric_fun(curr_compiled_predictions,\n",
    "            curr_in_sample_UPIs=np.unique(bootstrap_resamples[:,\n",
    "            bs_idx]), metric=metric) for bs_idx in iterator]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Surface function for BBC-CV metric calculation with parallel processing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def start_calculate_bs_metrics(\n",
    "    metric_fun,\n",
    "    curr_compiled_predictions,\n",
    "    num_bs,\n",
    "    n_cores,\n",
    "    metric='AUC',\n",
    "    progress_bar=True,\n",
    "    progress_bar_desc='',\n",
    "    ):\n",
    "\n",
    "    # Establish sizes of bootstrap replicates for each core\n",
    "\n",
    "    sizes = [num_bs // n_cores for _ in range(n_cores)]\n",
    "    sizes[-1] += num_bs - sum(sizes)\n",
    "\n",
    "    # Build arguments for metric sub-functions\n",
    "\n",
    "    arg_iterable = [(\n",
    "        metric_fun,\n",
    "        curr_compiled_predictions,\n",
    "        s,\n",
    "        metric,\n",
    "        progress_bar,\n",
    "        progress_bar_desc,\n",
    "        ) for s in sizes]\n",
    "\n",
    "    # Run metric sub-function in parallel\n",
    "\n",
    "    with multiprocessing.Pool(n_cores) as pool:\n",
    "        result = pool.starmap(_calculate_bs_metrics, arg_iterable)\n",
    "\n",
    "    return np.concatenate(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## III. Calculate metrics for GCSm threshold-level detection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 303
    },
    "executionInfo": {
     "elapsed": 25392911,
     "status": "error",
     "timestamp": 1619488950372,
     "user": {
      "displayName": "Shubhayu Bhattacharyay",
      "photoUrl": "",
      "userId": "16704285935129215396"
     },
     "user_tz": -60
    },
    "id": "W7JDp5NMp4vg",
    "outputId": "0b1f2573-ceae-441a-e6f5-f22111a37693"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Observation window: 0.05 hours completed.\n"
     ]
    }
   ],
   "source": [
    "# Extract directories of different observation windows\n",
    "\n",
    "obs_window_dirs = \\\n",
    "    glob.glob('../results/GCSm_threshold_prediction/*_h_obs_window/')\n",
    "\n",
    "# Iterate through observation window lengths\n",
    "\n",
    "for curr_obs_window_dir in obs_window_dirs[::-1]:\n",
    "\n",
    "    # Extract current observation window length from directory name\n",
    "\n",
    "    curr_obs_window_length = \\\n",
    "        float(re.search('GCSm_threshold_prediction/(.*)_h_obs_window/'\n",
    "              , curr_obs_window_dir).group(1))\n",
    "\n",
    "    # Status update on current observation window\n",
    "\n",
    "    print('Observation window: ' + str(curr_obs_window_length)\n",
    "        + ' hours started.')\n",
    "\n",
    "    # Extract compiled prediction file names in current observation window\n",
    "\n",
    "    thresh_pred_files = glob.glob(os.path.join(curr_obs_window_dir,\n",
    "                                  '*_compiled_predictions.csv'))\n",
    "\n",
    "    # Iterate through different threshold values\n",
    "\n",
    "    for curr_tresh_pred_file in thresh_pred_files:\n",
    "\n",
    "        # Extract current GCSm threshold from the file name\n",
    "\n",
    "        curr_GCSm_thresh = \\\n",
    "            re.search('_h_obs_window/(.*)_compiled_predictions.csv',\n",
    "                      curr_tresh_pred_file).group(1)\n",
    "\n",
    "        # Status update on current observation window\n",
    "\n",
    "        print('GCSm Threshold: ' + curr_GCSm_thresh + ' started.')\n",
    "\n",
    "        # Load current compiled prediction file for threshold\n",
    "\n",
    "        curr_compiled_predictions = pd.read_csv(curr_tresh_pred_file)\n",
    "\n",
    "        # Skip current GCSm threshold and observation window combination if no predictions are available\n",
    "\n",
    "        if curr_compiled_predictions.shape[0] == 0:\n",
    "            print('GCSm Threshold: ' + curr_GCSm_thresh\n",
    "                + ' skipped due to no predictions.')\n",
    "            continue\n",
    "\n",
    "        # Create new indicator for positive and negative cases\n",
    "\n",
    "        curr_compiled_predictions['PosCases'] = \\\n",
    "            curr_compiled_predictions.TrueLabel\n",
    "        curr_compiled_predictions['NegCases'] = \\\n",
    "            (curr_compiled_predictions.TrueLabel == 0).astype('int')\n",
    "\n",
    "        # Group predictions by UPI and calculate number of positive and negative cases by group\n",
    "\n",
    "        case_distribution = \\\n",
    "            pd.concat([curr_compiled_predictions.groupby('UPI'\n",
    "                      )['PosCases'].sum(),\n",
    "                      curr_compiled_predictions.groupby('UPI'\n",
    "                      )['NegCases'].sum()], axis=1)\n",
    "\n",
    "        # If there are not 2 or more patients for each case, skip current GCSm threshold and observation window combination\n",
    "\n",
    "        if ((case_distribution.PosCases >= 1).sum() < 2) \\\n",
    "            | ((case_distribution.NegCases >= 1).sum() < 2):\n",
    "            print('GCSm Threshold: ' + curr_GCSm_thresh\n",
    "                + ' skipped due to insufficient case distribution.')\n",
    "            continue\n",
    "\n",
    "        # Calculate AUC and ROC axes\n",
    "\n",
    "        bs_AUC_output = start_calculate_bs_metrics(\n",
    "            calculate_bs_AUC,\n",
    "            curr_compiled_predictions,\n",
    "            NUM_BOOTSTRAPS,\n",
    "            NUM_CORES,\n",
    "            metric='AUC',\n",
    "            progress_bar=True,\n",
    "            progress_bar_desc=curr_GCSm_thresh + ' AUC bootstrapping',\n",
    "            )\n",
    "        bs_AUC_df = pd.DataFrame(bs_AUC_output[:, 0:2],\n",
    "                                 columns=['Values', 'ConfigIdx'])\n",
    "        bs_AUC_df['Metrics'] = 'AUC'\n",
    "\n",
    "        # Calculate precision\n",
    "\n",
    "        bs_precision_output = start_calculate_bs_metrics(\n",
    "            calculate_bs_classification,\n",
    "            curr_compiled_predictions,\n",
    "            NUM_BOOTSTRAPS,\n",
    "            NUM_CORES,\n",
    "            metric='precision',\n",
    "            progress_bar=True,\n",
    "            progress_bar_desc=curr_GCSm_thresh\n",
    "                + ' precision bootstrapping',\n",
    "            )\n",
    "        bs_precision_df = pd.DataFrame(bs_precision_output,\n",
    "                columns=['Values', 'ConfigIdx'])\n",
    "        bs_precision_df['Metrics'] = 'precision'\n",
    "\n",
    "        # Calculate recall\n",
    "\n",
    "        bs_recall_output = start_calculate_bs_metrics(\n",
    "            calculate_bs_classification,\n",
    "            curr_compiled_predictions,\n",
    "            NUM_BOOTSTRAPS,\n",
    "            NUM_CORES,\n",
    "            metric='recall',\n",
    "            progress_bar=True,\n",
    "            progress_bar_desc=curr_GCSm_thresh + ' recall bootstrapping'\n",
    "                ,\n",
    "            )\n",
    "        bs_recall_df = pd.DataFrame(bs_recall_output, columns=['Values'\n",
    "                                    , 'ConfigIdx'])\n",
    "        bs_recall_df['Metrics'] = 'recall'\n",
    "\n",
    "        # Calculate f1-score\n",
    "\n",
    "        bs_f1_score_output = start_calculate_bs_metrics(\n",
    "            calculate_bs_classification,\n",
    "            curr_compiled_predictions,\n",
    "            NUM_BOOTSTRAPS,\n",
    "            NUM_CORES,\n",
    "            metric='f1-score',\n",
    "            progress_bar=True,\n",
    "            progress_bar_desc=curr_GCSm_thresh\n",
    "                + ' f1-score bootstrapping',\n",
    "            )\n",
    "        bs_f1_score_df = pd.DataFrame(bs_f1_score_output,\n",
    "                columns=['Values', 'ConfigIdx'])\n",
    "        bs_f1_score_df['Metrics'] = 'f1_score'\n",
    "\n",
    "        # Calculate specificity\n",
    "\n",
    "        bs_specificity_output = start_calculate_bs_metrics(\n",
    "            calculate_bs_classification,\n",
    "            curr_compiled_predictions,\n",
    "            NUM_BOOTSTRAPS,\n",
    "            NUM_CORES,\n",
    "            metric='specificity',\n",
    "            progress_bar=True,\n",
    "            progress_bar_desc=curr_GCSm_thresh\n",
    "                + ' specificity bootstrapping',\n",
    "            )\n",
    "        bs_specificity_df = pd.DataFrame(bs_specificity_output,\n",
    "                columns=['Values', 'ConfigIdx'])\n",
    "        bs_specificity_df['Metrics'] = 'specificity'\n",
    "\n",
    "        # Compile out-sample metrics into single dataframe\n",
    "\n",
    "        compiled_out_sample_metrics = pd.concat([bs_AUC_df,\n",
    "                bs_precision_df, bs_recall_df, bs_f1_score_df,\n",
    "                bs_specificity_df], ignore_index=True)\n",
    "\n",
    "        # Add metadata of current iteration to compiled dataframe\n",
    "\n",
    "        compiled_out_sample_metrics['Threshold'] = curr_GCSm_thresh\n",
    "        compiled_out_sample_metrics['ObsWindow'] = \\\n",
    "            curr_obs_window_length\n",
    "\n",
    "        # Reorder columns to intended order for subsequent R processing\n",
    "\n",
    "        compiled_out_sample_metrics = \\\n",
    "            compiled_out_sample_metrics[['ConfigIdx', 'Metrics',\n",
    "                'Values', 'Threshold', 'ObsWindow']]\n",
    "\n",
    "        # Save compiled metric dataframe in the correct observation window directory\n",
    "\n",
    "        compiled_out_sample_metrics.to_csv(os.path.join(curr_obs_window_dir,\n",
    "                curr_GCSm_thresh + '_compiled_metrics.csv'),\n",
    "                index=False)\n",
    "\n",
    "        # Extract ROC curve axis information and convert to dataframe\n",
    "\n",
    "        compiled_out_sample_ROC = pd.DataFrame({\n",
    "            'ConfigIdx': np.repeat(bs_AUC_output[:, 1], 200),\n",
    "            'FPR': np.concatenate(bs_AUC_output[:, 2]),\n",
    "            'TPR': np.concatenate(bs_AUC_output[:, 3]),\n",
    "            'Threshold': curr_GCSm_thresh,\n",
    "            'ObsWindow': curr_obs_window_length,\n",
    "            })\n",
    "\n",
    "        # Save the compiled ROC dataframe in the correct observation window directory\n",
    "\n",
    "        compiled_out_sample_ROC.to_csv(os.path.join(curr_obs_window_dir,\n",
    "                curr_GCSm_thresh + '_compiled_ROC.csv'), index=False)\n",
    "\n",
    "        # Clear output and signify end of threshold\n",
    "\n",
    "        clear_output()\n",
    "        print('Observation window: ' + str(curr_obs_window_length)\n",
    "            + ' hours started.')\n",
    "        print('GCSm Threshold: ' + curr_GCSm_thresh + ' completed.')\n",
    "\n",
    "    # Clear output and signify end of threshold\n",
    "\n",
    "    clear_output()\n",
    "    print('Observation window: ' + str(curr_obs_window_length)\n",
    "        + ' hours completed.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "executionInfo": {
     "elapsed": 811,
     "status": "ok",
     "timestamp": 1619453468689,
     "user": {
      "displayName": "Shubhayu Bhattacharyay",
      "photoUrl": "",
      "userId": "16704285935129215396"
     },
     "user_tz": -60
    },
    "id": "MT60z-FEp4vn"
   },
   "source": [
    "## IV. Calculate metrics for GOSE (at discharge) threshold-level prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 53
    },
    "executionInfo": {
     "elapsed": 5661170,
     "status": "ok",
     "timestamp": 1619053563006,
     "user": {
      "displayName": "Shubhayu Bhattacharyay",
      "photoUrl": "",
      "userId": "16704285935129215396"
     },
     "user_tz": -60
    },
    "id": "41SqxT9Qp4vn",
    "outputId": "83006717-a0de-4df5-8893-5ee9821999bd"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Observation window: 0.05 hours completed.\n"
     ]
    }
   ],
   "source": [
    "# Extract directories of different observation windows\n",
    "\n",
    "obs_window_dirs = \\\n",
    "    glob.glob('../results/GOSE_threshold_prediction/*_h_obs_window/')\n",
    "\n",
    "# Iterate through observation window lengths\n",
    "\n",
    "for curr_obs_window_dir in obs_window_dirs[::-1]:\n",
    "\n",
    "    # Extract current observation window length from directory name\n",
    "\n",
    "    curr_obs_window_length = \\\n",
    "        float(re.search('GOSE_threshold_prediction/(.*)_h_obs_window/'\n",
    "              , curr_obs_window_dir).group(1))\n",
    "\n",
    "    # Status update on current observation window\n",
    "\n",
    "    print('Observation window: ' + str(curr_obs_window_length)\n",
    "        + ' hours started.')\n",
    "\n",
    "    # Extract compiled prediction file names in current observation window\n",
    "\n",
    "    thresh_pred_files = glob.glob(os.path.join(curr_obs_window_dir,\n",
    "                                  '*_compiled_predictions.csv'))\n",
    "\n",
    "    # Iterate through different threshold values\n",
    "\n",
    "    for curr_tresh_pred_file in thresh_pred_files:\n",
    "\n",
    "        # Extract current GOSE threshold from the file name\n",
    "\n",
    "        curr_GOSE_thresh = \\\n",
    "            re.search('_h_obs_window/(.*)_compiled_predictions.csv',\n",
    "                      curr_tresh_pred_file).group(1)\n",
    "\n",
    "        # Status update on current observation window\n",
    "\n",
    "        print('GOSE Threshold: ' + curr_GOSE_thresh + ' started.')\n",
    "\n",
    "        # Load current compiled prediction file for threshold\n",
    "\n",
    "        curr_compiled_predictions = pd.read_csv(curr_tresh_pred_file)\n",
    "\n",
    "        # Skip current GOSE threshold and observation window combination if no predictions are available\n",
    "\n",
    "        if curr_compiled_predictions.shape[0] == 0:\n",
    "            print('GOSE Threshold: ' + curr_GOSE_thresh\n",
    "                + ' skipped due to no predictions.')\n",
    "            continue\n",
    "\n",
    "        # Create new indicator for positive and negative cases\n",
    "\n",
    "        curr_compiled_predictions['PosCases'] = \\\n",
    "            curr_compiled_predictions.TrueLabel\n",
    "        curr_compiled_predictions['NegCases'] = \\\n",
    "            (curr_compiled_predictions.TrueLabel == 0).astype('int')\n",
    "\n",
    "        # Group predictions by UPI and calculate number of positive and negative cases by group\n",
    "\n",
    "        case_distribution = \\\n",
    "            pd.concat([curr_compiled_predictions.groupby('UPI'\n",
    "                      )['PosCases'].sum(),\n",
    "                      curr_compiled_predictions.groupby('UPI'\n",
    "                      )['NegCases'].sum()], axis=1)\n",
    "\n",
    "        # If there are not 2 or more patients for each case, skip current GOSE threshold and observation window combination\n",
    "\n",
    "        if ((case_distribution.PosCases >= 1).sum() < 2) \\\n",
    "            | ((case_distribution.NegCases >= 1).sum() < 2):\n",
    "            print('GOSE Threshold: ' + curr_GOSE_thresh\n",
    "                + ' skipped due to insufficient case distribution.')\n",
    "            continue\n",
    "\n",
    "        # Calculate AUC and ROC axes\n",
    "\n",
    "        bs_AUC_output = start_calculate_bs_metrics(\n",
    "            calculate_bs_AUC,\n",
    "            curr_compiled_predictions,\n",
    "            NUM_BOOTSTRAPS,\n",
    "            NUM_CORES,\n",
    "            metric='AUC',\n",
    "            progress_bar=True,\n",
    "            progress_bar_desc=curr_GOSE_thresh + ' AUC bootstrapping',\n",
    "            )\n",
    "        bs_AUC_df = pd.DataFrame(bs_AUC_output[:, 0:2],\n",
    "                                 columns=['Values', 'ConfigIdx'])\n",
    "        bs_AUC_df['Metrics'] = 'AUC'\n",
    "\n",
    "        # Calculate precision\n",
    "\n",
    "        bs_precision_output = start_calculate_bs_metrics(\n",
    "            calculate_bs_classification,\n",
    "            curr_compiled_predictions,\n",
    "            NUM_BOOTSTRAPS,\n",
    "            NUM_CORES,\n",
    "            metric='precision',\n",
    "            progress_bar=True,\n",
    "            progress_bar_desc=curr_GOSE_thresh\n",
    "                + ' precision bootstrapping',\n",
    "            )\n",
    "        bs_precision_df = pd.DataFrame(bs_precision_output,\n",
    "                columns=['Values', 'ConfigIdx'])\n",
    "        bs_precision_df['Metrics'] = 'precision'\n",
    "\n",
    "        # Calculate recall\n",
    "\n",
    "        bs_recall_output = start_calculate_bs_metrics(\n",
    "            calculate_bs_classification,\n",
    "            curr_compiled_predictions,\n",
    "            NUM_BOOTSTRAPS,\n",
    "            NUM_CORES,\n",
    "            metric='recall',\n",
    "            progress_bar=True,\n",
    "            progress_bar_desc=curr_GOSE_thresh + ' recall bootstrapping'\n",
    "                ,\n",
    "            )\n",
    "        bs_recall_df = pd.DataFrame(bs_recall_output, columns=['Values'\n",
    "                                    , 'ConfigIdx'])\n",
    "        bs_recall_df['Metrics'] = 'recall'\n",
    "\n",
    "        # Calculate f1-score\n",
    "\n",
    "        bs_f1_score_output = start_calculate_bs_metrics(\n",
    "            calculate_bs_classification,\n",
    "            curr_compiled_predictions,\n",
    "            NUM_BOOTSTRAPS,\n",
    "            NUM_CORES,\n",
    "            metric='f1-score',\n",
    "            progress_bar=True,\n",
    "            progress_bar_desc=curr_GOSE_thresh\n",
    "                + ' f1-score bootstrapping',\n",
    "            )\n",
    "        bs_f1_score_df = pd.DataFrame(bs_f1_score_output,\n",
    "                columns=['Values', 'ConfigIdx'])\n",
    "        bs_f1_score_df['Metrics'] = 'f1_score'\n",
    "\n",
    "        # Calculate specificity\n",
    "\n",
    "        bs_specificity_output = start_calculate_bs_metrics(\n",
    "            calculate_bs_classification,\n",
    "            curr_compiled_predictions,\n",
    "            NUM_BOOTSTRAPS,\n",
    "            NUM_CORES,\n",
    "            metric='specificity',\n",
    "            progress_bar=True,\n",
    "            progress_bar_desc=curr_GOSE_thresh\n",
    "                + ' specificity bootstrapping',\n",
    "            )\n",
    "        bs_specificity_df = pd.DataFrame(bs_specificity_output,\n",
    "                columns=['Values', 'ConfigIdx'])\n",
    "        bs_specificity_df['Metrics'] = 'specificity'\n",
    "\n",
    "        # Compile out-sample metrics into single dataframe\n",
    "\n",
    "        compiled_out_sample_metrics = pd.concat([bs_AUC_df,\n",
    "                bs_precision_df, bs_recall_df, bs_f1_score_df,\n",
    "                bs_specificity_df], ignore_index=True)\n",
    "\n",
    "        # Add metadata of current iteration to compiled dataframe\n",
    "\n",
    "        compiled_out_sample_metrics['Threshold'] = curr_GOSE_thresh\n",
    "        compiled_out_sample_metrics['ObsWindow'] = \\\n",
    "            curr_obs_window_length\n",
    "\n",
    "        # Reorder columns to intended order for subsequent R processing\n",
    "\n",
    "        compiled_out_sample_metrics = \\\n",
    "            compiled_out_sample_metrics[['ConfigIdx', 'Metrics',\n",
    "                'Values', 'Threshold', 'ObsWindow']]\n",
    "\n",
    "        # Save compiled metric dataframe in the correct observation window directory\n",
    "\n",
    "        compiled_out_sample_metrics.to_csv(os.path.join(curr_obs_window_dir,\n",
    "                curr_GOSE_thresh + '_compiled_metrics.csv'),\n",
    "                index=False)\n",
    "\n",
    "        # Extract ROC curve axis information and convert to dataframe\n",
    "\n",
    "        compiled_out_sample_ROC = pd.DataFrame({\n",
    "            'ConfigIdx': np.repeat(bs_AUC_output[:, 1], 200),\n",
    "            'FPR': np.concatenate(bs_AUC_output[:, 2]),\n",
    "            'TPR': np.concatenate(bs_AUC_output[:, 3]),\n",
    "            'Threshold': curr_GOSE_thresh,\n",
    "            'ObsWindow': curr_obs_window_length,\n",
    "            })\n",
    "\n",
    "        # Save the compiled ROC dataframe in the correct observation window directory\n",
    "\n",
    "        compiled_out_sample_ROC.to_csv(os.path.join(curr_obs_window_dir,\n",
    "                curr_GOSE_thresh + '_compiled_ROC.csv'), index=False)\n",
    "\n",
    "        # Clear output and signify end of threshold\n",
    "\n",
    "        clear_output()\n",
    "        print('Observation window: ' + str(curr_obs_window_length)\n",
    "            + ' hours started.')\n",
    "        print('GOSE Threshold: ' + curr_GOSE_thresh + ' completed.')\n",
    "\n",
    "    # Clear output and signify end of threshold\n",
    "\n",
    "    clear_output()\n",
    "    print('Observation window: ' + str(curr_obs_window_length)\n",
    "        + ' hours completed.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "executionInfo": {
     "elapsed": 811,
     "status": "ok",
     "timestamp": 1619453468689,
     "user": {
      "displayName": "Shubhayu Bhattacharyay",
      "photoUrl": "",
      "userId": "16704285935129215396"
     },
     "user_tz": -60
    },
    "id": "MT60z-FEp4vn"
   },
   "source": [
    "## V. Calculate metrics for GOSE (at 12 months) threshold-level prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 53
    },
    "executionInfo": {
     "elapsed": 5661170,
     "status": "ok",
     "timestamp": 1619053563006,
     "user": {
      "displayName": "Shubhayu Bhattacharyay",
      "photoUrl": "",
      "userId": "16704285935129215396"
     },
     "user_tz": -60
    },
    "id": "41SqxT9Qp4vn",
    "outputId": "83006717-a0de-4df5-8893-5ee9821999bd"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Observation window: 1.0 hours completed.\n"
     ]
    }
   ],
   "source": [
    "# Extract directories of different observation windows\n",
    "\n",
    "obs_window_dirs = \\\n",
    "    glob.glob('../results/GOSE12m_threshold_prediction/*_h_obs_window/')\n",
    "\n",
    "# Iterate through observation window lengths\n",
    "\n",
    "for curr_obs_window_dir in obs_window_dirs:\n",
    "\n",
    "    # Extract current observation window length from directory name\n",
    "\n",
    "    curr_obs_window_length = \\\n",
    "        float(re.search('GOSE12m_threshold_prediction/(.*)_h_obs_window/'\n",
    "              , curr_obs_window_dir).group(1))\n",
    "\n",
    "    # Status update on current observation window\n",
    "\n",
    "    print('Observation window: ' + str(curr_obs_window_length)\n",
    "        + ' hours started.')\n",
    "\n",
    "    # Extract compiled prediction file names in current observation window\n",
    "\n",
    "    thresh_pred_files = glob.glob(os.path.join(curr_obs_window_dir,\n",
    "                                  '*_compiled_predictions.csv'))\n",
    "\n",
    "    # Iterate through different threshold values\n",
    "\n",
    "    for curr_tresh_pred_file in thresh_pred_files:\n",
    "\n",
    "        # Extract current GOSE12m threshold from the file name\n",
    "\n",
    "        curr_GOSE12m_thresh = \\\n",
    "            re.search('_h_obs_window/(.*)_compiled_predictions.csv',\n",
    "                      curr_tresh_pred_file).group(1)\n",
    "\n",
    "        # Status update on current observation window\n",
    "\n",
    "        print('GOSE12m Threshold: ' + curr_GOSE12m_thresh + ' started.')\n",
    "\n",
    "        # Load current compiled prediction file for threshold\n",
    "\n",
    "        curr_compiled_predictions = pd.read_csv(curr_tresh_pred_file)\n",
    "\n",
    "        # Skip current GOSE12m threshold and observation window combination if no predictions are available\n",
    "\n",
    "        if curr_compiled_predictions.shape[0] == 0:\n",
    "            print('GOSE12m Threshold: ' + curr_GOSE12m_thresh\n",
    "                + ' skipped due to no predictions.')\n",
    "            continue\n",
    "\n",
    "        # Create new indicator for positive and negative cases\n",
    "\n",
    "        curr_compiled_predictions['PosCases'] = \\\n",
    "            curr_compiled_predictions.TrueLabel\n",
    "        curr_compiled_predictions['NegCases'] = \\\n",
    "            (curr_compiled_predictions.TrueLabel == 0).astype('int')\n",
    "\n",
    "        # Group predictions by UPI and calculate number of positive and negative cases by group\n",
    "\n",
    "        case_distribution = \\\n",
    "            pd.concat([curr_compiled_predictions.groupby('UPI'\n",
    "                      )['PosCases'].sum(),\n",
    "                      curr_compiled_predictions.groupby('UPI'\n",
    "                      )['NegCases'].sum()], axis=1)\n",
    "\n",
    "        # If there are not 2 or more patients for each case, skip current GOSE12m threshold and observation window combination\n",
    "\n",
    "        if ((case_distribution.PosCases >= 1).sum() < 2) \\\n",
    "            | ((case_distribution.NegCases >= 1).sum() < 2):\n",
    "            print('GOSE12m Threshold: ' + curr_GOSE12m_thresh\n",
    "                + ' skipped due to insufficient case distribution.')\n",
    "            continue\n",
    "\n",
    "        # Calculate AUC and ROC axes\n",
    "\n",
    "        bs_AUC_output = start_calculate_bs_metrics(\n",
    "            calculate_bs_AUC,\n",
    "            curr_compiled_predictions,\n",
    "            NUM_BOOTSTRAPS,\n",
    "            NUM_CORES,\n",
    "            metric='AUC',\n",
    "            progress_bar=True,\n",
    "            progress_bar_desc=curr_GOSE12m_thresh + ' AUC bootstrapping',\n",
    "            )\n",
    "        bs_AUC_df = pd.DataFrame(bs_AUC_output[:, 0:2],\n",
    "                                 columns=['Values', 'ConfigIdx'])\n",
    "        bs_AUC_df['Metrics'] = 'AUC'\n",
    "\n",
    "        # Calculate precision\n",
    "\n",
    "        bs_precision_output = start_calculate_bs_metrics(\n",
    "            calculate_bs_classification,\n",
    "            curr_compiled_predictions,\n",
    "            NUM_BOOTSTRAPS,\n",
    "            NUM_CORES,\n",
    "            metric='precision',\n",
    "            progress_bar=True,\n",
    "            progress_bar_desc=curr_GOSE12m_thresh\n",
    "                + ' precision bootstrapping',\n",
    "            )\n",
    "        bs_precision_df = pd.DataFrame(bs_precision_output,\n",
    "                columns=['Values', 'ConfigIdx'])\n",
    "        bs_precision_df['Metrics'] = 'precision'\n",
    "\n",
    "        # Calculate recall\n",
    "\n",
    "        bs_recall_output = start_calculate_bs_metrics(\n",
    "            calculate_bs_classification,\n",
    "            curr_compiled_predictions,\n",
    "            NUM_BOOTSTRAPS,\n",
    "            NUM_CORES,\n",
    "            metric='recall',\n",
    "            progress_bar=True,\n",
    "            progress_bar_desc=curr_GOSE12m_thresh + ' recall bootstrapping'\n",
    "                ,\n",
    "            )\n",
    "        bs_recall_df = pd.DataFrame(bs_recall_output, columns=['Values'\n",
    "                                    , 'ConfigIdx'])\n",
    "        bs_recall_df['Metrics'] = 'recall'\n",
    "\n",
    "        # Calculate f1-score\n",
    "\n",
    "        bs_f1_score_output = start_calculate_bs_metrics(\n",
    "            calculate_bs_classification,\n",
    "            curr_compiled_predictions,\n",
    "            NUM_BOOTSTRAPS,\n",
    "            NUM_CORES,\n",
    "            metric='f1-score',\n",
    "            progress_bar=True,\n",
    "            progress_bar_desc=curr_GOSE12m_thresh\n",
    "                + ' f1-score bootstrapping',\n",
    "            )\n",
    "        bs_f1_score_df = pd.DataFrame(bs_f1_score_output,\n",
    "                columns=['Values', 'ConfigIdx'])\n",
    "        bs_f1_score_df['Metrics'] = 'f1_score'\n",
    "\n",
    "        # Calculate specificity\n",
    "\n",
    "        bs_specificity_output = start_calculate_bs_metrics(\n",
    "            calculate_bs_classification,\n",
    "            curr_compiled_predictions,\n",
    "            NUM_BOOTSTRAPS,\n",
    "            NUM_CORES,\n",
    "            metric='specificity',\n",
    "            progress_bar=True,\n",
    "            progress_bar_desc=curr_GOSE12m_thresh\n",
    "                + ' specificity bootstrapping',\n",
    "            )\n",
    "        bs_specificity_df = pd.DataFrame(bs_specificity_output,\n",
    "                columns=['Values', 'ConfigIdx'])\n",
    "        bs_specificity_df['Metrics'] = 'specificity'\n",
    "\n",
    "        # Compile out-sample metrics into single dataframe\n",
    "\n",
    "        compiled_out_sample_metrics = pd.concat([bs_AUC_df,\n",
    "                bs_precision_df, bs_recall_df, bs_f1_score_df,\n",
    "                bs_specificity_df], ignore_index=True)\n",
    "\n",
    "        # Add metadata of current iteration to compiled dataframe\n",
    "\n",
    "        compiled_out_sample_metrics['Threshold'] = curr_GOSE12m_thresh\n",
    "        compiled_out_sample_metrics['ObsWindow'] = \\\n",
    "            curr_obs_window_length\n",
    "\n",
    "        # Reorder columns to intended order for subsequent R processing\n",
    "\n",
    "        compiled_out_sample_metrics = \\\n",
    "            compiled_out_sample_metrics[['ConfigIdx', 'Metrics',\n",
    "                'Values', 'Threshold', 'ObsWindow']]\n",
    "\n",
    "        # Save compiled metric dataframe in the correct observation window directory\n",
    "\n",
    "        compiled_out_sample_metrics.to_csv(os.path.join(curr_obs_window_dir,\n",
    "                curr_GOSE12m_thresh + '_compiled_metrics.csv'),\n",
    "                index=False)\n",
    "\n",
    "        # Extract ROC curve axis information and convert to dataframe\n",
    "\n",
    "        compiled_out_sample_ROC = pd.DataFrame({\n",
    "            'ConfigIdx': np.repeat(bs_AUC_output[:, 1], 200),\n",
    "            'FPR': np.concatenate(bs_AUC_output[:, 2]),\n",
    "            'TPR': np.concatenate(bs_AUC_output[:, 3]),\n",
    "            'Threshold': curr_GOSE12m_thresh,\n",
    "            'ObsWindow': curr_obs_window_length,\n",
    "            })\n",
    "\n",
    "        # Save the compiled ROC dataframe in the correct observation window directory\n",
    "\n",
    "        compiled_out_sample_ROC.to_csv(os.path.join(curr_obs_window_dir,\n",
    "                curr_GOSE12m_thresh + '_compiled_ROC.csv'), index=False)\n",
    "\n",
    "        # Clear output and signify end of threshold\n",
    "\n",
    "        clear_output()\n",
    "        print('Observation window: ' + str(curr_obs_window_length)\n",
    "            + ' hours started.')\n",
    "        print('GOSE12m Threshold: ' + curr_GOSE12m_thresh + ' completed.')\n",
    "\n",
    "    # Clear output and signify end of threshold\n",
    "\n",
    "    clear_output()\n",
    "    print('Observation window: ' + str(curr_obs_window_length)\n",
    "        + ' hours completed.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## VI. Calculate precision-recall information for GOSE > 5 (at discharge) prediction with 6-hour observation window"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "GOSE.gt.5 PRC bootstrapping: 100%|██████████| 100/100 [00:10<00:00,  9.58it/s]\n",
      "GOSE.gt.5 PRC bootstrapping: 100%|██████████| 100/100 [00:10<00:00,  9.56it/s]\n",
      "GOSE.gt.5 PRC bootstrapping: 100%|██████████| 100/100 [00:10<00:00,  9.63it/s]\n",
      "GOSE.gt.5 PRC bootstrapping: 100%|██████████| 100/100 [00:10<00:00,  9.66it/s]\n",
      "GOSE.gt.5 PRC bootstrapping: 100%|██████████| 100/100 [00:10<00:00,  9.70it/s]\n",
      "GOSE.gt.5 PRC bootstrapping: 100%|██████████| 100/100 [00:10<00:00, 12.07it/s]\n",
      "GOSE.gt.5 PRC bootstrapping: 100%|██████████| 100/100 [00:10<00:00,  9.62it/s]\n",
      "GOSE.gt.5 PRC bootstrapping: 100%|██████████| 100/100 [00:10<00:00,  9.80it/s]\n",
      "GOSE.gt.5 PRC bootstrapping: 100%|██████████| 100/100 [00:10<00:00,  9.67it/s]\n",
      "GOSE.gt.5 PRC bootstrapping: 100%|██████████| 100/100 [00:10<00:00,  9.69it/s]\n"
     ]
    }
   ],
   "source": [
    "# Extract directories of different observation windows\n",
    "\n",
    "curr_compiled_predictions = pd.read_csv('../results/GOSE_threshold_prediction/06.00_h_obs_window/GOSE.gt.5_compiled_predictions.csv')\n",
    "\n",
    "# Create new indicator for positive and negative cases\n",
    "\n",
    "curr_compiled_predictions['PosCases'] = \\\n",
    "    curr_compiled_predictions.TrueLabel\n",
    "curr_compiled_predictions['NegCases'] = \\\n",
    "    (curr_compiled_predictions.TrueLabel == 0).astype('int')\n",
    "\n",
    "# Group predictions by UPI and calculate number of positive and negative cases by group\n",
    "\n",
    "case_distribution = \\\n",
    "    pd.concat([curr_compiled_predictions.groupby('UPI'\n",
    "              )['PosCases'].sum(),\n",
    "              curr_compiled_predictions.groupby('UPI'\n",
    "              )['NegCases'].sum()], axis=1)\n",
    "\n",
    "# Calculate precision recall information in parallel\n",
    "\n",
    "bs_PRC_output = start_calculate_bs_metrics(\n",
    "    calculate_bs_PRC,\n",
    "    curr_compiled_predictions,\n",
    "    NUM_BOOTSTRAPS,\n",
    "    NUM_CORES,\n",
    "    metric='AUPRC',\n",
    "    progress_bar=True,\n",
    "    progress_bar_desc='GOSE.gt.5' + ' PRC bootstrapping',\n",
    "    )\n",
    "bs_AUPRC_df = pd.DataFrame(bs_PRC_output[:, 0:2],\n",
    "                           columns=['Values', 'ConfigIdx'])\n",
    "bs_AUPRC_df['Metrics'] = 'AUPRC'\n",
    "\n",
    "# Add metadata of current iteration to compiled dataframe\n",
    "\n",
    "bs_AUPRC_df['Threshold'] = 'GOSE.gt.5'\n",
    "bs_AUPRC_df['ObsWindow'] = 6\n",
    "\n",
    "# Reorder columns to intended order for subsequent R processing\n",
    "\n",
    "bs_AUPRC_df = \\\n",
    "    bs_AUPRC_df[['ConfigIdx', 'Metrics',\n",
    "        'Values', 'Threshold', 'ObsWindow']]\n",
    "\n",
    "# Save AUPRC dataframe in the correct observation window directory\n",
    "\n",
    "bs_AUPRC_df.to_csv('../results/GOSE_threshold_prediction/GOSE.gt.5_AUPRC.csv',\n",
    "                   index=False)\n",
    "\n",
    "# Extract Precision-Recall curve axis information and convert to dataframe\n",
    "\n",
    "compiled_out_sample_ROC = pd.DataFrame({\n",
    "    'ConfigIdx': np.repeat(bs_PRC_output[:, 1], 1000),\n",
    "    'Recall': np.concatenate(bs_PRC_output[:, 2]),\n",
    "    'Precision': np.concatenate(bs_PRC_output[:, 3]),\n",
    "    'Threshold': 'GOSE.gt.5',\n",
    "    'ObsWindow': 6,\n",
    "    })\n",
    "\n",
    "# Save the Precision Recall dataframe in the correct directory\n",
    "\n",
    "compiled_out_sample_ROC.to_csv('../results/GOSE_threshold_prediction/GOSE.gt.5_precision_recall_curve.csv',\n",
    "                               index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.8366745283018868"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Load baseline characteristics\n",
    "baseline_characteristics = pd.read_csv('../clinical_data/patient_baseline_characteristics.csv')\n",
    "\n",
    "# Load outcomes\n",
    "outcomes = pd.read_csv('../clinical_data/patient_outcomes.csv')\n",
    "\n",
    "# Merge baseline characteristics with outcomes\n",
    "auc_df = pd.merge(baseline_characteristics[['UPI','APACHEIIMortalityRisk']],outcomes[['UPI','DiedDuringHospitalStay']],on='UPI',how='left')\n",
    "\n",
    "roc_auc_score(auc_df.DiedDuringHospitalStay,auc_df.APACHEIIMortalityRisk)\n",
    "#accuracy_score(auc_df.DiedDuringHospitalStay,(auc_df.APACHEIIMortalityRisk > .5).astype('int'))"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "07_calculate_metrics_v2.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "sb_venv",
   "language": "python",
   "name": "sb_venv"
  },
  "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}