ALD-ML/ALD_ML.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "# Proteomics-based biomarker identification for alcohol-related liver disease (ALD)\n",
    "\n",
    "**Outline**\n",
    "\n",
    "1. Prediction targets are three pathological features of liver disease:\n",
    "    1. Fibrosis\n",
    "    1. Inflammation\n",
    "    2. Steatosis\n",
    "    \n",
    "2. Two data sets with \n",
    "    1. clinical markers\n",
    "    2. Plasma proteomics \n",
    "    \n",
    "**Highlighted Contents**\n",
    "> In order to jump to highlighted sections, use a table of contents plugin (toc) for a structured view ([lab](https://github.com/jupyterlab/jupyterlab-toc) | [notebook](https://jupyter-contrib-nbextensions.readthedocs.io/en/latest/nbextensions/toc2/README.html)). \n",
    "> Without data, but for easier navigation try to go to [colab](https://colab.research.google.com/).\n",
    "\n",
    "2. Explore datasets\n",
    "    1. Proteomics data and clinical data\n",
    "        - available from the authors upon reasonable request. Intersted researchers may contact Odense Patient Data Exploratory Network open@rsyd.dk to obtain relevant data\n",
    "        \n",
    "3. Models\n",
    "    1. (3.2) Individual Models for three endpoints fibrosis, steatosis and inflammation\n",
    "        - Cross-validation results\n",
    "    2. (3.4) Final Model\n",
    "        - final model used for DeLong-Test comparison, clinical follow-up evaluation, and classification in the validation cohort\n",
    "\n",
    "\n",
    "> Some data is hidden from the public output until it is cleared. `#hide`Odense Patient Data Exploratory Network open@rsyd.dk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "import os\n",
    "from pathlib import Path\n",
    "CPUS = os.cpu_count()\n",
    "RANDOMSTATE = 29\n",
    "\n",
    "FOLDER_DATA_RAW = 'data/raw'\n",
    "DATAFOLDER = 'data/processed'\n",
    "os.makedirs('data/processed', exist_ok=True)\n",
    "TABLEFOLDER = 'tables'\n",
    "RESULT_FOLDER = 'results'\n",
    "FIGURE_FOLDER = Path('figures')\n",
    "FIGURE_FOLDER.mkdir(exist_ok=True)\n",
    "\n",
    "import logging\n",
    "logger = logging.getLogger()\n",
    "logger.setLevel(logging.INFO)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "import sklearn\n",
    "import sklearn.linear_model as skllm\n",
    "from src.pandas import combine_value_counts\n",
    "import sklearn.metrics as sklm\n",
    "import sklearn.model_selection as sklms\n",
    "from sklearn.model_selection import cross_val_predict\n",
    "from sklearn.metrics import roc_auc_score, roc_curve\n",
    "\n",
    "import ipywidgets as widgets\n",
    "import tqdm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "import src"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "# Explore datasets\n",
    "\n",
    "Diagnostic comparators (existing best-in-class) biomarkers\n",
    "- Fibrosis markers: transient elastography, 2-dimensional shear wave elastography, ELF test, FibroTest, FIB4 score, APRI score, Forns score, ProC3\n",
    "- Inflammation markers: M30=caspase-cleaved cytokeratin-18 fragments, M65=total CK18, AST:ALT ratio, ProC3\n",
    "- Steatosis: Controlled attenuation parameter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.set_option('max_columns', 9)\n",
    "\n",
    "_folder = FOLDER_DATA_RAW\n",
    "_index_col = 'Sample ID'\n",
    "\n",
    "files = [file for file in os.listdir(_folder) if '.csv' in file]\n",
    "\n",
    "if not files:\n",
    "    print(\"No processed files found.\")\n",
    "else:\n",
    "    w_data = widgets.Dropdown(options=files)\n",
    "\n",
    "    show_data = src.widgets.create_show_data(index_col=_index_col, datafolder=_folder)\n",
    "    out = widgets.interactive_output(show_data, controls={'file':w_data})\n",
    "\n",
    "    data = show_data.__closure__[0].cell_contents\n",
    "    w_cols = widgets.SelectMultiple(options=list(data.columns))\n",
    "\n",
    "    show_selected_proteins = src.widgets.create_show_selected_proteins(data=data)\n",
    "\n",
    "    out_sel = widgets.interactive_output(show_selected_proteins, {'columns': w_cols})\n",
    "    out_sel = widgets.interactive_output(show_selected_proteins, {'columns': w_cols})\n",
    "\n",
    "    # Updater\n",
    "    def widget_updater(other_widget):\n",
    "        \"\"\"Picks first element from closure. other_widget is not used directly\"\"\"\n",
    "        data = show_data.__closure__[0].cell_contents\n",
    "        w_cols.options = list(data.columns)\n",
    "        show_selected_proteins.__closure__[0].cell_contents = data\n",
    "\n",
    "    _ = widgets.interactive_output(widget_updater, {'other_widget': w_data})\n",
    "\n",
    "    display(widgets.VBox([w_data, out, w_cols, out_sel]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "## Proteomics data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "### Load Complete proteomics data\n",
    "\n",
    "Full preprocessing pipeline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "annotation_file = pd.read_csv(os.path.join(FOLDER_DATA_RAW, 'Experiment annotation file_upgrade.csv'), index_col = [0])\n",
    "annotation_file_plasma = annotation_file[annotation_file['Sample type'] == 'Plasma']\n",
    "annotation_file_plasma.index = pd.Index(annotation_file_plasma.index, dtype=int)\n",
    "display(annotation_file_plasma.head())\n",
    "annotation_file_plasma.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The annotation file holds the filename for the the processed raw data by Skyline and some annotation, e.g. the `Sample ID`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The mapping for Protein ID to the gene ID is given by `report_plasma`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "report_plasma = pd.read_csv(os.path.join(FOLDER_DATA_RAW, 'Protein_ALDupgrade_Report.csv'), na_values='Filtered')\n",
    "report_plasma.rename({'PG.Genes': 'Gene names', \n",
    "                      'PG.ProteinAccessions': 'Protein ID', \n",
    "                      'PG.ProteinDescriptions':'Protein names'}, inplace= True, axis=1)\n",
    "IDmapping_ProteinID_to_Proteinname=dict(zip(report_plasma['Protein ID'], report_plasma['Protein names']))\n",
    "report_plasma=report_plasma.drop('Protein names', axis=1)\n",
    "report_plasma.head() "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create mapping object (see if proteins are unique -> get function from other project)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "columns_ = ['Protein ID', 'Gene names']\n",
    "ids_ = report_plasma[columns_].apply(lambda series_: series_.str.split(';'))\n",
    "ids_.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "def length_(x):\n",
    "    try:\n",
    "        return len(x)\n",
    "    except:\n",
    "        return 0\n",
    "    \n",
    "count_groups_proteins = ids_.apply(lambda series_: series_.apply(length_))\n",
    "count_groups_proteins"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "from src.pandas import combine_value_counts\n",
    "combine_value_counts(count_groups_proteins)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Protein IDs are always set\n",
    "- two proteins have no annotations (Gene name count of 0 appears twice)\n",
    "- the are some protein names which are mapped to the same gene. \n",
    "Let have a look at cases where a set of proteins was not mapped uniquely to one gene:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ids_.loc[count_groups_proteins['Gene names'] != 1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "from src.imputation import imputation_normal_distribution, log2, NP_LOG_FCT, IMPUTATION_MEAN_SHIFT, IMPUTATION_STD_SHRINKAGE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "experimental_columns = annotation_file_plasma['Sample ID']\n",
    "report_plasma[columns_] = report_plasma[columns_].apply(lambda series_: series_.str.split(';').str[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "report_plasma.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "map_filenames_ids = dict(zip(annotation_file['File name'], annotation_file['Sample ID']))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "- Remove some measurements which are not intensities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "data_plasma_raw = report_plasma.copy()\n",
    "data_plasma_raw.drop(data_plasma_raw.filter(regex='StrippedSequences').columns, axis=1, inplace = True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "- rename column names to sample ID from annotation file\n",
    "- set index to proteins"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "data_plasma_raw = data_plasma_raw.rename(mapper = map_filenames_ids, axis=1)\n",
    "IDmapping_UniprotID_to_Genename = dict(zip(data_plasma_raw['Protein ID'], data_plasma_raw['Gene names']))\n",
    "data_plasma_raw = data_plasma_raw.set_index('Protein ID').drop('Gene names', axis = 1)\n",
    "data_plasma_raw.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "mask = data_plasma_raw.notna().sum(axis=1) > 603 * 0.6\n",
    "mask.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "- Filter at protein level for 60% data completeness across all runs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "DATA_COMPLETENESS = 0.6\n",
    "data_plasma_filtered = data_plasma_raw.dropna(axis=0, thresh = data_plasma_raw.shape[1] * DATA_COMPLETENESS)\n",
    "# data_plasma_filtered #hide"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Filter at sample level for at least 200 proteins quantified"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask_filtered_out = data_plasma_filtered.notna().sum() < 200\n",
    "if mask_filtered_out.sum()>0:\n",
    "    data_plasma_filtered.loc[:, list(mask_filtered_out)].describe().loc['count'].astype(int).sort_values()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if mask_filtered_out.sum()>0:\n",
    "    data_plasma_raw.loc[:, mask_filtered_out].describe().loc['count'].astype(int).sort_values()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Cutoff of 118 is next one where another sample would be discarded."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false",
    "tags": []
   },
   "outputs": [],
   "source": [
    "MIN_N_PROTEIN_GROUPS = 200\n",
    "print(f\"Min No. of Protein-Groups in single sample: {MIN_N_PROTEIN_GROUPS}, i.e. a fraction of {MIN_N_PROTEIN_GROUPS/len(data_plasma_filtered)}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "data_plasma_filtered = data_plasma_filtered.dropna(axis=1, thresh = MIN_N_PROTEIN_GROUPS)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "assert (data_plasma_filtered.dtypes != float).sum() == 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "data_plasma_filtered_log = data_plasma_filtered.apply(log2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "##### Imputation\n",
    "\n",
    "- imputation is done before coefficient of variation (CV)\n",
    "- is this sensible?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "SCALE_DATA = False\n",
    "if SCALE_DATA:\n",
    "    from sklearn.preprocessing import StandardScaler\n",
    "\n",
    "    scaler = StandardScaler()\n",
    "\n",
    "    data_plasma_filtered_log_imputed_np = scaler.fit_transform(data_plasma_filtered_log.values)\n",
    "    data_plasma_filtered_log_imputed = data_plasma_filtered_log.copy()\n",
    "    data_plasma_filtered_log_imputed.loc[:,:] = np.nan_to_num(data_plasma_filtered_log_imputed_np)\n",
    "else:\n",
    "    data_plasma_filtered_log_imputed = data_plasma_filtered_log.apply(imputation_normal_distribution)\n",
    "    assert data_plasma_filtered_log_imputed.loc['A0A024R6I7', 'Plate5_D5'] - 10.355115 < 0.0001, 'Imputed value changed in comparison to previous run'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_plasma_filtered_log_imputed.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "file = Path('data/processed/plasma_processed.csv')\n",
    "file.parent.mkdir(parents=True, exist_ok=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "import logging\n",
    "try:\n",
    "    data_plasma_filtered_log_imputed.to_csv(file.absolute())\n",
    "except PermissionError as e:\n",
    "    logging.warning(f\"No write permission to directory: {e}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "qc_plasma = annotation_file_plasma[annotation_file_plasma['Group2'] == 'QC']['Sample ID']\n",
    "df_qc = data_plasma_filtered.copy()[qc_plasma]\n",
    "coef_of_variation = lambda x: np.std(x) / np.mean(x)\n",
    "proteins_cv = df_qc.apply(coef_of_variation, axis = 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false",
    "tags": []
   },
   "outputs": [],
   "source": [
    "CV_COEFFICIENT = 0.4\n",
    "cv_selected = proteins_cv < CV_COEFFICIENT\n",
    "print(f\"Selected proteins # {cv_selected.sum()} of a total of # {len(cv_selected)}!\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "df_qc = df_qc.assign(cv = proteins_cv)\n",
    "qc_cv = df_qc[cv_selected].index\n",
    "\n",
    "df = data_plasma_filtered_log_imputed.copy()\n",
    "df = df.rename_axis('Sample ID', axis=1).T\n",
    "df_cv = df[qc_cv]\n",
    "data_proteomics = df_cv.drop('A0A0B4J2D9', axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "PROTEOM  = 'data_ml_proteomics_cleaned.csv'\n",
    "\n",
    "data_proteomics.to_csv(os.path.join(DATAFOLDER, PROTEOM))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"A maximum of {1} proteins in {0} samples can be used for proteomic models\".format(*data_proteomics.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Low intensities** below 8 (in log-scale)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "intensities_below_8 = data_proteomics[data_proteomics < 8].dropna(how='all').dropna(how='all', axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "- Data proteomics is the summary of the following processing steps:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "summary_protein_preprocessing = [(\"Proportion protein has to be shared between samples\" , DATA_COMPLETENESS),\n",
    "                                 (\"Minimum number of protein in single sample\", MIN_N_PROTEIN_GROUPS),\n",
    "                                 (\"Maximum coefficient of variation (CV) for protein intensities\", CV_COEFFICIENT),\n",
    "                                 (\"Logarithm employed for transformation\", NP_LOG_FCT),\n",
    "                                 (\"Imputation: Mean-Shift\", IMPUTATION_MEAN_SHIFT), \n",
    "                                 (\"Imputation: Std-Dev. shrinkage\", IMPUTATION_STD_SHRINKAGE)\n",
    "                                ]\n",
    "\n",
    "for descr, value in summary_protein_preprocessing:\n",
    "    print('{}: {}'.format(descr, value))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "### Load Protein GeneID Mapping\n",
    "\n",
    "- UniProtID to Gene name mapping"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "key_ProteinID = pd.read_csv(os.path.join(FOLDER_DATA_RAW, 'ID_matching_key_ALDupgrade.csv'), \n",
    "                            index_col=\"Protein ID\")\n",
    "key_ProteinID.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Note that it is possible multiple protein IDs are mapped to the same gene."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ids_.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert len(key_ProteinID) == len(ids_), \"Both references should match at least in the number of proteins. \""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "## Clinical data\n",
    "### Load Complete clinical data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "CLINICAL = 'df_cli_164.csv'\n",
    "COL_ID = 'Sample ID'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "CLINICAL = 'df_cli_164.csv'\n",
    "COL_ID = 'Sample ID'\n",
    "\n",
    "f_data_clinic = os.path.join(FOLDER_DATA_RAW, CLINICAL)\n",
    "data_cli = pd.read_csv(f_data_clinic, index_col=COL_ID)\n",
    "data_cli = data_cli[data_cli['kleiner']!=0.5]\n",
    "# data_cli #hide"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "annotation_file[annotation_file['Sample ID']=='Plate6_F2']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "w_cols_cli = widgets.SelectMultiple(options=list(data_cli.columns))\n",
    "\n",
    "def show_selected_markers(columns):\n",
    "    if len(columns)> 0:\n",
    "        display(data_cli[list(w_cols_cli.value)])\n",
    "        display(data_cli[list(w_cols_cli.value)].describe())\n",
    "    else:\n",
    "        print('Select clinical markers')\n",
    "\n",
    "out_cli = widgets.interactive_output(show_selected_markers, {'columns': w_cols_cli})\n",
    "widgets.VBox([w_cols_cli, out_cli])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "### Selected Clinical markers\n",
    "\n",
    "Diagnostic comparators (existing best-in-class) biomarkers\n",
    "- state-of-the-art (**SOTA**) Fibrosis markers: \n",
    "    - `te`: transient elastography (sona liver scan)\n",
    "    - `swe`: 2-dimensional shear wave elastography\n",
    "    - `elf`: ELF test\n",
    "    - `ft`: FibroTest\n",
    "    - `fib4`: FIB4 score\n",
    "    - `apri`: APRI score\n",
    "    - `forns`: Forns score\n",
    "    - `p3np`: ProC3\n",
    "- Inflammation markers:\n",
    "    - M30=caspase-cleaved cytokeratin-18 fragments\n",
    "    - M65=total CK18\n",
    "    - AST:ALT ratio\n",
    "    - ProC3\n",
    "- Steatosis: Controlled attenuation parameter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "#SOTA_fibrosis = ['te', 'swe', 'elf', 'ft', 'fib4', 'apri', 'forns', 'p3np']\n",
    "SOTA_fibrosis = ['elf', 'ft', 'fib4', 'apri', 'forns', 'p3np']\n",
    "data_cli.groupby('kleiner')[SOTA_fibrosis].count()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "pd.set_option('max_columns', 20)\n",
    "FEATURES_ML = ['nas_steatosis_ordinal', 'nas_inflam', 'kleiner', \n",
    "          'fib4', 'elf', 'ft', 'te', 'swe', 'aar','ast',\n",
    "          'apri','forns','m30', 'm65', 'meld', 'p3np', 'timp1', 'cap' ]\n",
    "# data_cli[FEATURES_ML].head() #hide"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "data_cli.groupby('group2')[FEATURES_ML].count()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "SOTA_fibrosis = ['te', 'swe', 'elf', 'ft', 'fib4', 'apri', 'forns', 'p3np']\n",
    "data_cli.groupby('kleiner')[SOTA_fibrosis].median()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Selected Demographics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "demographics = data_cli[['age', 'bmi', 'gender_num']] # 1 is male\n",
    "demographics.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SELECTED_DEMOGRAPHICS = ['age', 'gender_num']\n",
    "data_cli[SELECTED_DEMOGRAPHICS].head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "### Targets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "fibrosis_score = data_cli.kleiner\n",
    "inflamation_score = data_cli.nas_inflam\n",
    "steatosis_score = data_cli.nas_steatosis_ordinal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "TARGETS = ['kleiner', 'nas_steatosis_ordinal', 'nas_inflam']\n",
    "Y = data_cli[TARGETS]\n",
    "Y.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from src.pandas import combine_value_counts\n",
    "\n",
    "freq_targets = combine_value_counts(Y)\n",
    "freq_targets.loc['Total',:] = freq_targets.sum()\n",
    "freq_targets.to_excel(os.path.join(TABLEFOLDER, 'freq_endpoints_unique_values.xlsx'))\n",
    "freq_targets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "Several binary features can be created.\n",
    "\n",
    "target      | Scale   | unique values              | Binarization                 |  N samples\n",
    "-----       | --------| ---------------            | -------------------------    |  ---------\n",
    "fibrosis    | five    | F0, F1, F2, F3, F4         | (F0,F1) vs (F2, F3, F4)      |  360\n",
    "fibrosis    | five    | F0, F1, F2, F3, F4         | (F0,F1,F2) vs (F3, F4)       |  360\n",
    "inflamation | seven   | I0, I1, I2, I3, I4, I5     | (I0, I1) vs (I2, I3, I4, I5) |  352\n",
    "steatosis   | five    | S0, S1, S2, S3             | (S0) vs (S1, S2, S3)         |  352\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "Variable naming: `<target>_greater-equal_<value>`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "from src.pandas import create_dichotome\n",
    "kleiner_ge_2     = create_dichotome(Y['kleiner'], 2)\n",
    "kleiner_ge_3     = create_dichotome(Y['kleiner'], 3)\n",
    "steatosis_ge_1   = create_dichotome(Y['nas_steatosis_ordinal'], 1)\n",
    "inflamation_ge_2 = create_dichotome(Y['nas_inflam'], 2)\n",
    "\n",
    "end_points = ['F2', 'F3', 'S1', 'I2']\n",
    "dichotomies = [kleiner_ge_2, kleiner_ge_3, steatosis_ge_1, inflamation_ge_2]\n",
    "targets_dict  = {k: v for k, v in zip(end_points, dichotomies)}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Frequencies of binary variables:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "freq_targets = pd.DataFrame(\n",
    "    {'kleiner>=2': kleiner_ge_2.value_counts(dropna=False, sort=False),\n",
    "     'kleiner>=3': kleiner_ge_3.value_counts(dropna=False, sort=False),\n",
    "     'steatosis>=1' : steatosis_ge_1.value_counts(dropna=False, sort=False),\n",
    "     'inflamation>=2':inflamation_ge_2.value_counts(dropna=False, sort=False)\n",
    "    })\n",
    "freq_targets.loc['total'] = freq_targets.sum()\n",
    "freq_targets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "### Clinical Cutoffs for targets\n",
    "\n",
    "Cutoff for binary grouping of targets\n",
    "\n",
    "target      | Scale   | unique values            | N samples\n",
    "----------- | ------- | ----------------------   | -------\n",
    "fibrosis    | five    | F0, F1, F2, F3, F4       | 360\n",
    "steatosis   | five    | S0, S1, S2, S3           | 352\n",
    "inflamation | seven   | I0, I1, I2, I3, I4, I5   | 352\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "file_cutoff_clinic = os.path.join(FOLDER_DATA_RAW, \"clinical_marker_test_cut-offs.xlsx\")\n",
    "cutoffs_clinic = pd.read_excel(file_cutoff_clinic, sheet_name=\"cutoffs\", index_col='marker')\n",
    "cutoffs_clinic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "markers_to_drop = []\n",
    "for marker in cutoffs_clinic.index:\n",
    "    if marker not in data_cli.columns:\n",
    "        print(f\"{marker}: Missing in clinics data.\")\n",
    "        markers_to_drop.append(marker)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "`proc3` is not in data_clinic. drop this from the list of cutoffs! (Cutoff can be learned later)\n",
    "Rename columns to desired endpoint name."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "if markers_to_drop:\n",
    "    cutoffs_clinic.drop(labels=markers_to_drop, inplace=True)\n",
    "cutoffs_clinic.columns = ['F2', 'F3', 'I2', 'S1']\n",
    "cutoffs_clinic"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "Extract certain cutoff for binary targets defined by column name:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "cutoffs_clinic['F2'].dropna().to_dict()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "See statistics (e.g. median) of SOTA-markers for clinical fibrosis assessment (represented by categories 0 to 4)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "SOTA_fibrosis = ['te', 'swe', 'elf', 'ft', 'fib4', 'apri', 'forns', 'p3np']\n",
    "data_cli.groupby('kleiner')[SOTA_fibrosis].median()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "### Handle missing features of clinical data (Global Missing Pattern)\n",
    "\n",
    "> No imputation of clinical features for now as only single clinical features are used in \"univariate\" models. Imputation is only sensible if several types of information are combined. Then one could use [`sklearn.impute.simpleImputer`](https://scikit-learn.org/stable/modules/impute.html)'s default `'mean'` strategy or alternatively one could replace missing values with zeros on the standardised data to zero mean and standard deviation of one.\n",
    "\n",
    "Features are present to widely different degree. In order to be able to define global splits with the same pattern of missings over the features and targets by samples (here: patients), we define a missing pattern for stratification."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "FEATURES_CLINIC = cutoffs_clinic.index\n",
    "data_cli[FEATURES_CLINIC].describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We keep only samples for which any target is present. The other could be later used for verification of model prediction in the clinic."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "patient_ids_w_target = data_cli[TARGETS].dropna(how='all').index\n",
    "print(f\"No. of samples without target variable: {len(data_cli) -len(patient_ids_w_target)} \")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now define the set of variables of which we want to define missingness patterns:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "FEATURES_CLINIC_ALL = list(FEATURES_CLINIC) + SELECTED_DEMOGRAPHICS + TARGETS\n",
    "data_cli.loc[patient_ids_w_target, FEATURES_CLINIC_ALL].describe().sort_values(by=\"count\", ascending=False, axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false",
    "tags": []
   },
   "outputs": [],
   "source": [
    "def ordered_missing_table(data:pd.DataFrame):\n",
    "    \"\"\"Order dataframe by data completeness (first column has most features) \n",
    "    and then return an encoding of completeness (1 = available, 0 0 not available)\"\"\"\n",
    "       \n",
    "    data_missing_table = data.notna().astype(int)\n",
    "    var_ordered_by_completness = list(data.describe().loc['count'].sort_values(ascending=False).index)\n",
    "    data_missing_table = data_missing_table.sort_values(by=var_ordered_by_completness)[var_ordered_by_completness]\n",
    "    return data_missing_table.replace(0, pd.NA).convert_dtypes()\n",
    "\n",
    "print(\"Used features: {}\".format(\", \".join(FEATURES_CLINIC_ALL)))\n",
    "data_cli_missing_table = ordered_missing_table(data=data_cli.loc[patient_ids_w_target, FEATURES_CLINIC_ALL])\n",
    "data_cli_missing_table = data_cli_missing_table.dropna(how='all', axis=0).dropna(how='all', axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_cli_missing_table.describe().loc['count'].astype(int)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare both available data for proteomics and clinical features. We will add the availability of proteomics data as another feature to our missingness patterns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_proteomics.isna().any(axis=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "in_both = data_proteomics.index.intersection(data_cli_missing_table.index)\n",
    "samples_wo_proteomics_data = data_cli_missing_table.index.difference(in_both)\n",
    "print(\"{} diagnosed patients have no valid proteome measure: {}\".format(\n",
    "    len(samples_wo_proteomics_data), \n",
    "    \", \".join(samples_wo_proteomics_data)\n",
    "))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "HAS_QUANT_PROT = 'has_prot'\n",
    "data_cli_missing_table[HAS_QUANT_PROT] = pd.Series(1, index=data_proteomics.index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_cli_missing_table.dropna(how='all').describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_cli_missing_table = ordered_missing_table(data_cli_missing_table)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_cli_missing_strings = data_cli_missing_table.fillna(value=0)\n",
    "data_cli_missing_strings = data_cli_missing_strings.astype(str)\n",
    "stratifier = data_cli_missing_strings.apply(lambda x : x.str.cat(), axis=1)\n",
    "# display(stratifier.head()) #hide\n",
    "stratifier_tab = stratifier.value_counts()\n",
    "stratifier_tab"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will have to get ride of the singletons (unique value only once observed). Possibly the grouping could be extended to the values up to 5."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "unique_missing_patterns = list(stratifier_tab.index)\n",
    "\n",
    "def match_observed(seq1, seq2):\n",
    "    return sum(pos1 == pos2 for pos1, pos2 in zip(seq1, seq2))\n",
    "\n",
    "assert match_observed(\"111111110011100001\", \"111111110011110010\") == 15, \"Failed\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stratifier.value_counts().min()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def update_stratifier(stratifier_var:pd.Series, threshold:int=None, verbose:bool=False):\n",
    "    \"\"\"Takes a stratifier variable, and assigns the pattern the less \n",
    "    often observed (defined by threshold or the minimum) to the closest other missing pattern.\n",
    "    Clossness is defined by the number of features which are present/absent for the samples.\"\"\"\n",
    "    stratifier_var =stratifier_var.copy()\n",
    "    stratifier_tab = stratifier_var.value_counts()\n",
    "    current_minimum = stratifier_tab.min()\n",
    "    if threshold is not None  and current_minimum >= threshold:\n",
    "        logger.info(\"Threshold already reached.\")\n",
    "        return stratifier_var\n",
    "    unique_missing_patterns = list(stratifier_tab.index)\n",
    "    list_single_missing_patterns = stratifier_tab[stratifier_tab <= current_minimum].index\n",
    "    for single_missing_pattern in list_single_missing_patterns:\n",
    "        if verbose:\n",
    "            logger.info(f\"Find match for: {single_missing_pattern}\")\n",
    "        closest = 0\n",
    "        for i, other_seq in enumerate(unique_missing_patterns):\n",
    "            if not other_seq in list_single_missing_patterns:\n",
    "                relatedness = match_observed(single_missing_pattern, other_seq)\n",
    "                if relatedness > closest:\n",
    "                    closest = relatedness\n",
    "                    best = other_seq\n",
    "        stratifier_var[stratifier_var == single_missing_pattern] = best\n",
    "        if verbose:\n",
    "            logger.info(f\"Best match is : {best}\")\n",
    "    if threshold is not None:\n",
    "        stratifier_tab = stratifier_var.value_counts()\n",
    "        new_minimum = stratifier_tab.min()\n",
    "        if new_minimum < threshold:\n",
    "            stratifier_var = update_stratifier(stratifier_var, threshold=threshold)        \n",
    "    return stratifier_var\n",
    "# stratifier = update_stratifier(stratifier).value_counts()\n",
    "stratifier = update_stratifier(stratifier, threshold=5, verbose=True)\n",
    "stratifier.value_counts()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Global stratification based on string. It won't be possible to distribute unique cases, which is why they were assigned to their closest papern. Then split models between endpoints are comparable as they are subsets of the global splits. This  garuantess:\n",
    "1. By endpoint: Metrics for marker on one test set does not contain patients which are in the training set of a different marker.\n",
    "2. By marker: Metrics for different endpoints on the test set does not contain training samples of a different endpoint."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import RepeatedStratifiedKFold\n",
    "\n",
    "CV_FOLDS = 5\n",
    "CV_REPEATS = 10\n",
    "\n",
    "RANDOM_SEED = 123\n",
    "\n",
    "rskf = RepeatedStratifiedKFold(n_splits=CV_FOLDS, n_repeats=CV_REPEATS, random_state=RANDOM_SEED) \n",
    "splits = list(rskf.split(data_cli_missing_table, stratifier))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cv_train_test_indices = list()\n",
    "for train_indices, test_indices in splits:\n",
    "    cv_train_test_indices.append(\n",
    "     (stratifier.index[train_indices], stratifier.index[test_indices])\n",
    "    )\n",
    "cv_train_test_indices[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "## Visualization of data\n",
    "\n",
    "Look at UMAPs with labels from disease categories.\n",
    "  - Does the assigned disease correspond to certain groups\n",
    " \n",
    "For clinical data, on could look at a selection of scatter plots in order to see if it is feasible to separate some groups based on two features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#ToDo\n",
    "#import umap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "# Models\n",
    "\n",
    "Different _experimental_ setups for prediction models will be compared. First, for the target **fibrosis**. Fibrosis is reported on a five-point scale from stage F0 to F4.\n",
    "\n",
    "ML setup binary    | HP  | F0  | F1  | F2  | F3  | F4\n",
    "--- | --- | ---    | --- | --- | --- | ---\n",
    "HP-F0-F2 vs F3-F4  | c   | c   | c   | c   | t   | t    \n",
    "F0-F2 vs F3-F4 (advanced)    |     | c   | c   | c   | t   | t\n",
    "F0-F1 vs F2-F4 (significant)    |     | c   | c   | t   | t   | t\n",
    "\n",
    "In the table, c stands for control  and t for target. The clinical relevance is to distinguish different \n",
    "stages of disease. The question is wheater one should include a healthy, untested patient cohort can help building a \n",
    "classification model, as e.g. for fibrosis the general prevalence in the population is between 6 to 7 percent. Alternatively a _multi-task model_ with having 5 classes/end-points can be fit.\n",
    "\n",
    "\n",
    "In addition to fibrosis, the endpoints **steatosis** and **inflamation** can be predicted.\n",
    "\n",
    "target      | Scale   | unique values              | N samples\n",
    "-----       | --------| ---------------            | -------\n",
    "fibrosis    | five    | F0, F1, F2, F3, F4         | \n",
    "steatosis   | five    | S0, S1, S2, S3, S4         | \n",
    "inflamation | seven   | I0, I1, I2, I3, I4, I5, I6 | \n",
    "\n",
    "\n",
    "What is population of interest?\n",
    "- population at risk\n",
    "- general population (which we do not have as a \"random\" sample)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "## Preparation: Classifiers and Evaluation Functionality"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "### Predefined (Tree-based classifiers, SVMs, GLMs)\n",
    "- Select Classifier by cross-validation using [sklearn functionality](https://scikit-learn.org/stable/model_selection.html#model-selection)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "clf_lr    = skllm.LogisticRegression(random_state=0, solver='liblinear')\n",
    "clf_lr_key = 'Logistic'\n",
    "\n",
    "# specify more sklearn classifiers if you need\n",
    "clf_sklearn = {clf_lr_key: clf_lr}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "[Refitting](https://scikit-learn.org/stable/tutorial/basic/tutorial.html#refitting-and-updating-parameters) the same estimator by invocing it `fit`-method overwrites the previously learned weights."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "### Custom Threshold-based classification\n",
    "Create a classifier based on the threshold which is compatible with the basic scikit-learn functionality, see [instructions](https://scikit-learn.org/stable/developers/develop.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Example for using the cutoff of Fibrosis >=2 from the cutoff-table:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cutoffs_clinic.loc['te','F2']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false",
    "tags": []
   },
   "outputs": [],
   "source": [
    "from src.threshold_classifier import ThresholdClassifier\n",
    "clf_te = ThresholdClassifier(threshold={'te':7.0})\n",
    "print(clf_te.threshold)\n",
    "clf_te.fit(data_cli.fillna(value=0))\n",
    "y_pred = clf_te.predict(data_cli)\n",
    "clf_te.predict_proba(data_cli)[:4] # no scores, either 0 or 1 as cutoff is just compared"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "Note: Having only one feature for threshold classification does make the definiton of a cutoff unnecessary. AUC-ROC statistics are not meaningful withouth scores."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "# from sklearn.utils.estimator_checks import check_estimator\n",
    "# check_estimator(ThresholdClassifier)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "true"
   },
   "source": [
    "### Confusion Matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "from src.scoring import ConfusionMatrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "How to use it, e.g. for using a clinical marker cutoff for fibrosis:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "y_true = data_cli.kleiner > 2.0\n",
    "\n",
    "# y_pred defined as Threshold-example\n",
    "cm_f2_te = ConfusionMatrix(y_true, y_pred)\n",
    "print(\"As DataFrame:\")\n",
    "display(cm_f2_te.as_dataframe)\n",
    "print(\"Plain:\\n\",cm_f2_te)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "### Cross-Validation Procedure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "*Note on Cross-Validation Procedure*\n",
    "- Comparing the performance on random splits of the entire data will lead to overconfident predicitons.\n",
    "- Performing the Cross-Validation only on a `Train`-split would allow to have a better evaluation on the test dataset. \n",
    "- Cutoff calibration would need a validation split\n",
    "\n",
    "\n",
    "##### Cutoff Specification\n",
    "> in clinical setting, false-alarms are preferrable than missed detections. Yes we should find a way to customize the cutoffs\n",
    "> to have a high sensitivity but also decent specificity, but I guess it risks over-tuning on this specific dataset?  \n",
    "> Can one ramp over and find the optimal based on F1 score? Would MCC be a better alternative?  \n",
    "> (Author?)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "annotation_file[annotation_file['Sample type']=='Plasma']['Group2'].value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "annotation_file['Groups'].value_counts()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Selected Metrics for Binary Classification evaluation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "from sklearn.model_selection import cross_validate\n",
    "from sklearn.metrics import roc_curve\n",
    "scoring = ['precision', 'recall', 'f1', 'balanced_accuracy', 'roc_auc']\n",
    "scoring"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We build a dictionary of the scoring functions for later use:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Define population in order to obtain comparable splits of the data\n",
    "\n",
    "The common functionality provided by `sklearn` does not allow for [nested stratification](https://stackoverflow.com/a/45526792/9684872). `RepeatedKFold` is splitting based on the data in the target variable.\n",
    "One solution is to encode a sample with missing target or feature values explicitly into the target variable, but this is not feasible for many different sets of feature and target variables (here each marker-variable with each target is a set).\n",
    "\n",
    "Data has to be aligned for computation. Due to missing values on some features, the runs are not directly comparable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_y = kleiner_ge_2\n",
    "_X = data_cli.te.to_frame().fillna(0)\n",
    "in_both = _y.index.intersection(_X.index)\n",
    "_X = _X.loc[in_both]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "e.g. for clinical marker"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false",
    "tags": []
   },
   "outputs": [],
   "source": [
    "from src.cross_validation import run_cv_binary_simple, _get_cv_means\n",
    "clf = {**{'f2_te': clf_te}, **clf_sklearn}\n",
    "print(\"Klassifiers:\", \", \".join(clf.keys()))\n",
    "\n",
    "result_dict = run_cv_binary_simple(clf, X=_X, y=_y, cv=5, scoring=scoring, return_estimator=True)\n",
    "result_dict.keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_get_cv_means(result_dict).sort_values(('test_f1', 'mean'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "> NOTE: The ROC_AUC value is misleading in case of the ThresholdClassification `f2_te` as the predictor does not yield probabilites ($y_{predicted} \\in \\{0,1\\}$)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Try to visualize Decision"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Univariate Logistic Regression\n",
    "For the univariate logistic regression\n",
    "$$ ln \\frac{p}{1-p} = \\beta_0 + \\beta_1 \\cdot x $$\n",
    "the cutoff `c=0.5` corresponds a feature value of: \n",
    "$$ x = - \\frac{\\beta_0}{\\beta_1} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for lr_est in result_dict['Logistic']['estimator']:\n",
    "    # lr_0 = result_dict['Logistic']['estimator'][0]    \n",
    "    print(f\"Custom cutoff defined by Logistic regressor: {- float(lr_est.intercept_) / float(lr_est.coef_):.2f} \")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "#### Rebuilded `run_cv_binary` to get roc_curve value\n",
    "\n",
    "The re-implemented interface for `run_cv_binary` has a a similar interface as sklearns [`cross_validate`](https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.cross_validate.html). The `group` parameter is missing as it's not used currently in this setup."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from src.cross_validation import run_cv_binary\n",
    "#run_cv_binary?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_dict, roc_curve_results, precision_recall_results = run_cv_binary(clf, X=_X, y=_y, cv=cv_train_test_indices, prefix='F2_', verbose=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# mean over flattend array\n",
    "assert np.mean(pd.DataFrame(results_dict).loc['y_test', 'F2_Logistic']) == np.mean(np.array(pd.DataFrame(results_dict).loc['y_test', 'F2_Logistic']).flatten())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "Display CV results (metrics):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "_get_cv_means(results_dict)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Averaging the predictions\n",
    "Get predictions for samples as average of the predictions on test set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_df = pd.DataFrame(index=_y.index)\n",
    "for _i, _y_pred in enumerate(results_dict['F2_Logistic']['y_test']):\n",
    "    _df[f'run_{_i:02}'] = _y_pred\n",
    "display(_df.head())\n",
    "results = _df.mean(axis=1).to_frame(name='mean')\n",
    "results['std'] = _df.std(axis=1, skipna=True)\n",
    "results['n_pred'] = _df.notna().sum(axis=1).astype(int)\n",
    "results.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "del results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Confidence-Intervals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result_metrics = _get_cv_means(results_dict)\n",
    "result_metrics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_95CI_df(df_metrics, selected_metrics=None):\n",
    "    \"\"\"Expects output from _get_cv_means.\"\"\"\n",
    "    if selected_metrics is None:\n",
    "        selected_metrics = result_metrics.columns.levels[0]\n",
    "    else:\n",
    "        assert set(selected_metrics) in set(result_metrics.columns.levels[0])\n",
    "    \n",
    "    key_lower_CI = 'lower'\n",
    "    key_upper_CI = 'upper'\n",
    "    \n",
    "    def _create_95CI_df(df_metric, mean_col='mean', std_col='std'):\n",
    "        \"\"\"Create from a DataFrame of results the 95% CI.\n",
    "        Lower and upper bound.\"\"\"\n",
    "        df_CI = pd.DataFrame(index=df_metric.index)\n",
    "        df_CI[key_lower_CI] = df_metric[mean_col] - 2*df_metric[std_col]\n",
    "        df_CI[key_upper_CI] = df_metric[mean_col] + 2*df_metric[std_col]\n",
    "        return df_CI\n",
    "\n",
    "    df_95CI = pd.DataFrame(index=df_metrics.index, columns=pd.MultiIndex.from_product([selected_metrics, ['lower', 'upper']]))\n",
    "    for _metric in selected_metrics:\n",
    "        df_95CI[_metric] = _create_95CI_df(df_metrics[_metric])\n",
    "    return df_95CI\n",
    "\n",
    "create_95CI_df(result_metrics)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "#### Extension: Learn imputation on fold\n",
    "Include Preprocessing (here: imputation into the pipeline). The imputation of the proteomics data would then be based only on moments learned the training data (splits) for the Gaussian distribution of each peptide.\n",
    "\n",
    "> write custom [`FunctionTransformer`](https://scikit-learn.org/stable/modules/preprocessing.html#custom-transformers) to included preprocessing.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "# #ToDo: Possible extension\n",
    "from sklearn.pipeline import make_pipeline\n",
    "from sklearn.preprocessing import FunctionTransformer\n",
    "\n",
    "def custom_preprocessing_function(X):\n",
    "    \"\"\"Operate on a a set of rows from the dataset.\n",
    "    here: apply imputation to log-transformed values?\n",
    "    \"\"\"\n",
    "    raise NotImplementedError\n",
    "\n",
    "# clf = make_pipeline(FunctionTransformer(custom_preprocessing_function), svm.SVC(C=1))\n",
    "# cross_val_score(clf, X, y, cv=cv)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "## Individual Models for three endpoints fibrosis, steatosis and inflammation\n",
    "Strategy for feature selection and model performance validation: 3 models to be built, fibrosis (F0-1 vs. F2-4, and F0-2 vs. F3-4), inflammation (0-1 vs. 2-5) and steatosis (0 vs. >0). Then compare each model with their respective existing best-in-class markers according to their standard cut-offs in clinic. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from src.sklearn import FeatureSelector\n",
    "    \n",
    "feature_selected = FeatureSelector(k=10, protein_gene_data=key_ProteinID)\n",
    "_proteins_selected_f2 = feature_selected.fit(data_proteomics, kleiner_ge_2, 'F2')\n",
    "_proteins_selected_f2.columns = ['F2 (k=10)']\n",
    "_proteins_selected_f2['F2 (k=5)'] = feature_selected.get_k_best('F2', 5)\n",
    "_proteins_selected_f2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Screen for optimized number of features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "[Feature selection](https://scikit-learn.org/stable/modules/feature_selection.html#feature-selection) based on mutual information"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.exceptions import UndefinedMetricWarning\n",
    "import warnings; warnings.simplefilter('ignore', UndefinedMetricWarning)\n",
    "\n",
    "RECALCULATE_FEATURES = False\n",
    "RESULT_FEATURE_COMPARISON = os.path.join(DATAFOLDER, 'summary_n_features_msproteomics_upgrade.pkl')\n",
    "def main_n_features_comparison(n_features_max=50):\n",
    "    \"compare performance using an grid of features\"      \n",
    "    from tqdm.notebook import tqdm as tqdm\n",
    "    from time import perf_counter as pc\n",
    "    t0 = pc()\n",
    "    summary = []\n",
    "    feature_selected = FeatureSelector(k=n_features_max, protein_gene_data=key_ProteinID)\n",
    "    _ = feature_selected.fit(data_proteomics, kleiner_ge_2, 'F2')\n",
    "    _ = feature_selected.fit(data_proteomics, kleiner_ge_3, 'F3')\n",
    "    _ = feature_selected.fit(data_proteomics, steatosis_ge_1, 'S1')\n",
    "    _ = feature_selected.fit(data_proteomics, inflamation_ge_2, 'I2')\n",
    "    for n_features in tqdm(range(1,n_features_max)): \n",
    "#         feature_selected = FeatureSelector(k=n_features, protein_gene_data=key_ProteinID)\n",
    "        proteins_selected_f2 = feature_selected.get_k_best('F2', n_features)\n",
    "        proteins_selected_f3 = feature_selected.get_k_best('F3', n_features)\n",
    "        proteins_selected_s1 = feature_selected.get_k_best('S1', n_features)\n",
    "        proteins_selected_I2 = feature_selected.get_k_best('I2', n_features)\n",
    "        test_cases = {}\n",
    "        test_cases['F2'] = {'proteins': proteins_selected_f2, 'y':kleiner_ge_2}\n",
    "        test_cases['F3'] = {'proteins': proteins_selected_f3, 'y':kleiner_ge_3}\n",
    "        test_cases['S1'] = {'proteins': proteins_selected_s1, 'y':steatosis_ge_1}\n",
    "        test_cases['I2'] = {'proteins': proteins_selected_I2, 'y':inflamation_ge_2}\n",
    "        for test_case in test_cases.keys():\n",
    "            _clf_key = 'LR'\n",
    "            _clf = skllm.LogisticRegression(random_state=0, solver='liblinear')\n",
    "            proteins_selected = test_cases[test_case]['proteins']\n",
    "            y = test_cases[test_case]['y']\n",
    "            _X = data_proteomics[proteins_selected.index]\n",
    "            in_both = y.index.intersection(_X.index)\n",
    "            _X = _X.loc[in_both]\n",
    "            _y = y.loc[in_both]\n",
    "#             result = cross_validate(_clf, X=_X, y=_y, groups=_y, cv=RepeatedStratifiedKFold(n_splits=CV_FOLDS, n_repeats=CV_REPEATS, random_state=RANDOM_SEED) , scoring=scoring)\n",
    "            result, _, _ = run_cv_binary({_clf_key:_clf}, X=_X, y=_y, cv=cv_train_test_indices, prefix=f'{test_cases[\"F2\"]}_', verbose=False)\n",
    "            _key = list(result.keys()).pop()\n",
    "            result = pd.DataFrame(result[_key])\n",
    "            result['name'] = _clf.__class__.__name__\n",
    "            result['n_features'] = n_features\n",
    "            result['test_case'] = test_case\n",
    "            summary.append(result)\n",
    "    summary = [pd.DataFrame(_) for _ in summary]\n",
    "    summary_n_features = pd.concat(summary)\n",
    "    summary_n_features.to_pickle(RESULT_FEATURE_COMPARISON) # long format\n",
    "    print(f\"Finished. Elapsed seconds {pc()-t0:.2f}\")\n",
    "    return summary_n_features\n",
    "\n",
    "if not RECALCULATE_FEATURES:\n",
    "    try:\n",
    "        summary_n_features = pd.read_pickle(RESULT_FEATURE_COMPARISON)\n",
    "    except FileNotFoundError:\n",
    "        summary_n_features = main_n_features_comparison()\n",
    "else:\n",
    "    summary_n_features = main_n_features_comparison() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(5,5))\n",
    "sns.lineplot(x='n_features',y='roc_auc',hue='test_case', data=summary_n_features)\n",
    "plt.ylim([0.5,1])\n",
    "plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.title('Number of Features vs roc auc')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig=plt.figure(figsize=(14.5, 4))\n",
    "metrics = ['f1', 'balanced_accuracy', 'roc_auc']\n",
    "for i in range(3):\n",
    "    plt.subplot(1, 3, i+1)\n",
    "    sns.lineplot(x='n_features',y=metrics[i],hue='test_case', data=summary_n_features[summary_n_features['test_case']!='F3'], palette=['darkblue', 'gray', 'darkred'])\n",
    "    plt.ylim(0, 1)\n",
    "    #plt.title('Number of Features vs {}'.format(metrics[i]), fontsize=14)\n",
    "    plt.ylabel(metrics[i], fontsize=14)\n",
    "    plt.xlabel('Number of features', fontsize=14)\n",
    "    plt.xticks(fontsize=14);\n",
    "    plt.yticks(fontsize=14);\n",
    "#plt.savefig('figures/1_vs_panel.png', dpi=120, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "combined = summary_n_features.groupby(['test_case','n_features']).mean()\n",
    "\n",
    "best = combined.sort_values(by='f1', ascending=False).groupby('test_case').head(1)\n",
    "best"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "single_protein = combined[combined['num_feat']==1]\n",
    "single_protein"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "one_vs_panel = pd.concat([best, single_protein]).sort_values(by='test_case')\n",
    "one_vs_panel.to_csv('tables/1_vs_panel.csv')\n",
    "one_vs_panel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "best_dict = {}\n",
    "for i, j in best.index.to_list():\n",
    "    best_dict[i] = j\n",
    "best_dict"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "### Top k selected proteins for prediction\n",
    "[Feature selection](https://scikit-learn.org/stable/modules/feature_selection.html#feature-selection) based on mutual information. \n",
    "\n",
    "Each endpoint will yield different `top-k` proteins. An aggregation strategy in the simplest form is to combine the top-k. Maybe there is also some kind of rank-algorithm combining the top-k minimizing the overall rank?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "protein_panels={}\n",
    "\n",
    "for end_point, dichotomy in targets_dict.items():\n",
    "    k=best_dict[end_point]\n",
    "    feature_selected = FeatureSelector(k=k, protein_gene_data=key_ProteinID)\n",
    "    protein_panels[end_point] = feature_selected.fit(data_proteomics, dichotomy, end_point)\n",
    "\n",
    "    \n",
    "proteins_selected_f2=protein_panels['F2']\n",
    "proteins_selected_f3=protein_panels['F3']\n",
    "proteins_selected_s1=protein_panels['S1']\n",
    "proteins_selected_I2=protein_panels['I2']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "proteins_selected_f2.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "proteins_selected_f3.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "proteins_selected_s1.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "proteins_selected_I2.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# don't label by gene to retain scoring information (to sort at least by endpoint for importance)\n",
    "protein_panels_scores={}\n",
    "for end_point, dichotomy in targets_dict.items():\n",
    "    k=best_dict[end_point]\n",
    "    feature_selected = FeatureSelector(k=k)\n",
    "    protein_panels_scores[end_point] = feature_selected.fit(data_proteomics, dichotomy, end_point)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_feature_comp(protein_panel:dict, exclude=[], order=None):\n",
    "    \"\"\"Custom processor for dictonary holding multual information DataFrame per endpoint\n",
    "        from above.\"\"\"\n",
    "    for i, (endpoint, _df) in enumerate(protein_panel.items()):\n",
    "        if endpoint not in exclude:\n",
    "            if i == 0:\n",
    "                df_protein_panel = _df\n",
    "            else:\n",
    "                df_protein_panel = df_protein_panel.join(_df, how='outer')\n",
    "    if order:\n",
    "        df_protein_panel = df_protein_panel[order]\n",
    "    mask = df_protein_panel.isna()   \n",
    "    # df_protein_panel.where(mask, other=1).fillna(0).sort_values(by=list, ascending=False)\n",
    "    return df_protein_panel.sort_values(by=list(df_protein_panel.columns), ascending=False).fillna('-')\n",
    "\n",
    "df_protein_panel = get_feature_comp(protein_panel=protein_panels_scores, exclude=['F3'], order=['F2', 'I2', 'S1'])\n",
    "df_protein_panel['Gene Name'] = key_ProteinID.loc[df_protein_panel.index]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.set_option('precision', 3)\n",
    "display(df_protein_panel.head(10))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(df_protein_panel)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### mrmr feature selction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mrmr import mrmr_classif"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from sklearn.model_selection import RepeatedStratifiedKFold, cross_validate\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from mrmr import mrmr_classif\n",
    "\n",
    "RESULT_FEATURE_COMPARISON_MRMR = os.path.join(DATAFOLDER, 'summary_n_features_proteomics_mrmr_upgrade.pkl')\n",
    "RECALCULATE_FEATURES = False\n",
    "\n",
    "def n_features_comparison_mrmr(n_features_max=50):\n",
    "    summary = []\n",
    "    test_cases = {}\n",
    "    test_cases['F2'] = {'y':kleiner_ge_2}\n",
    "    test_cases['F3'] = { 'y':kleiner_ge_3}\n",
    "    test_cases['S1'] = { 'y':steatosis_ge_1}\n",
    "    test_cases['I2'] = { 'y':inflamation_ge_2}    \n",
    "    for test_case in test_cases.keys():\n",
    "        cv = RepeatedStratifiedKFold(n_splits=5, n_repeats=10, random_state=42)\n",
    "        model = LogisticRegression(random_state=42, solver='liblinear')\n",
    "        scoring = ['precision', 'recall', 'f1', 'balanced_accuracy', 'roc_auc']\n",
    "        X = data_proteomics\n",
    "        y = test_cases[test_case]['y']\n",
    "        in_both = y.index.intersection(X.index)\n",
    "        _X = X.loc[in_both]\n",
    "        _y = y.loc[in_both]\n",
    "        for n_features in range(1, n_features_max):\n",
    "            selected_features = mrmr_classif(_X, _y, K=n_features)\n",
    "            _X_mrmr = _X[selected_features]\n",
    "            scores = cross_validate(model, _X_mrmr, _y, scoring=scoring, cv=cv)\n",
    "            scores['n_features'] = n_features\n",
    "            scores['test_case'] = test_case\n",
    "            scores['n_observations'] = _X.shape[0]\n",
    "            results = pd.DataFrame(scores)\n",
    "            summary.append(results)\n",
    "    summary = [pd.DataFrame(_) for _ in summary]\n",
    "    summary_n_features = pd.concat(summary)\n",
    "    summary_n_features.to_pickle(RESULT_FEATURE_COMPARISON_MRMR)\n",
    "    return(summary_n_features)\n",
    "\n",
    "if not RECALCULATE_FEATURES:\n",
    "    try:\n",
    "        summary_n_features_mrmr = pd.read_pickle(RESULT_FEATURE_COMPARISON_MRMR)\n",
    "    except FileNotFoundError:\n",
    "        summary_n_features_mrmr = n_features_comparison_mrmr()\n",
    "else:\n",
    "    summary_n_features_mrmr = n_features_comparison_mrmr()\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "combined_mrmr = summary_n_features_mrmr.groupby(['test_case','n_features']).mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(5,5))\n",
    "sns.lineplot(x='n_features',y='test_roc_auc',hue='test_case', data=summary_n_features_mrmr)\n",
    "plt.ylim([0.5,1])\n",
    "plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.title('Number of Features vs. ROC-AUC', fontsize=12)\n",
    "plt.xlabel('Number of features', fontsize =12)\n",
    "plt.ylabel('ROC-AUC', fontsize=12)\n",
    "plt.xticks(fontsize=12)\n",
    "plt.yticks(fontsize=12)\n",
    "#plt.savefig('figures/proteomics_mrmr.png', bbox_inches='tight', dpi=120)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "summary_n_features_mrmr.to_csv('data/processed/source_data_figS10a-c.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig=plt.figure(figsize=(14.5, 4))\n",
    "metrics = ['test_f1', 'test_balanced_accuracy', 'test_roc_auc']\n",
    "for i in range(3):\n",
    "    plt.subplot(1, 3, i+1)\n",
    "    sns.lineplot(x='n_features',y=metrics[i],hue='test_case', data=summary_n_features_mrmr[summary_n_features_mrmr['test_case']!='F3'], palette=['darkblue', 'gray', 'darkred'])\n",
    "    plt.ylim(0, 1)\n",
    "    plt.ylabel(metrics[i], fontsize=14)\n",
    "    plt.xlabel('Number of features', fontsize=14)\n",
    "    plt.xticks(fontsize=14);\n",
    "    plt.yticks(fontsize=14);\n",
    "    plt.ylim(0.4, 1)\n",
    "    plt.title('{} vs.\\n number of features'.format(metrics[i]), fontsize=14)\n",
    "plt.savefig('figures/1_vs_panel.png', dpi=120, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Top k mrmr selected proteins for prediction "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Based on roc auc plot shown above, the minimal number of features to achieve a \"good enough\" performance for \n",
    "-- F2 is 9\n",
    "-- F3 is 5\n",
    "-- I2 is 6\n",
    "-- S1 is 12"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_cases = {}\n",
    "test_cases['F2'] = {'n_features': 9, 'y':kleiner_ge_2}\n",
    "test_cases['F3'] = {'n_features': 5, 'y':kleiner_ge_3}\n",
    "test_cases['I2'] = {'n_features': 6, 'y':inflamation_ge_2} \n",
    "test_cases['S1'] = {'n_features': 12, 'y':steatosis_ge_1}\n",
    "\n",
    "\n",
    "mrmr_markers = {}\n",
    "for test_case in test_cases.keys():\n",
    "    X = data_proteomics\n",
    "    y = test_cases[test_case]['y']\n",
    "    in_both = y.index.intersection(X.index)\n",
    "    _X = X.loc[in_both]\n",
    "    _y = y.loc[in_both]\n",
    "    markers = mrmr_classif(_X, _y, K=test_cases[test_case]['n_features'])\n",
    "    markers = pd.DataFrame({'Protein ID':markers, test_case: [IDmapping_UniprotID_to_Genename[i]for i in markers]}).set_index('Protein ID')\n",
    "    mrmr_markers[test_case]=markers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "FILE_MARKERS_UPGRADE = os.path.join(TABLEFOLDER, 'marker panels mrmr upgrade.xlsx')\n",
    "\n",
    "with pd.ExcelWriter(FILE_MARKERS_UPGRADE) as writer:\n",
    "    for test_case in ['F2', 'F3', 'I2', 'S1']:\n",
    "        mrmr_markers[test_case].to_excel(writer, sheet_name=test_case+'feature_mrmr')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "best_mrmr = combined_mrmr.loc[[('F2', 9),('I2', 6), ('S1', 12)]]\n",
    "best_mrmr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "single_protein_mrmr = combined_mrmr.loc[[('F2', 1), ('I2', 1), ('S1', 1)]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "one_vs_panel_mrmr = pd.concat([best_mrmr, single_protein_mrmr]).sort_values(by='test_case')\n",
    "one_vs_panel_mrmr.to_csv('tables/1_vs_panel_mrmr.csv')\n",
    "one_vs_panel_mrmr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Inspect marker distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df=data_plasma_filtered_log.rename_axis('Sample ID', axis=1).T.copy().reset_index()\n",
    "df2=data_cli.copy().reset_index()\n",
    "df2=df2[df2['kleiner']!=0.5]\n",
    "df3=df.merge(df2[['kleiner', 'nas_inflam','nas_steatosis_ordinal', 'Sample ID']], on='Sample ID').set_index('Sample ID')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "source_data_figs7a = df3[mrmr_markers['F2'].index.tolist() + ['kleiner']].melt(id_vars=['kleiner'], value_name='Log2 intensity', var_name='Protein ID').dropna(subset=['kleiner', 'Log2 intensity'])\n",
    "source_data_figs7b = df3[mrmr_markers['I2'].index.tolist() + ['nas_inflam']].melt(id_vars=['nas_inflam'], value_name='Log2 intensity', var_name='Protein ID').dropna(subset=['nas_inflam', 'Log2 intensity'])\n",
    "source_data_figs7c = df3[mrmr_markers['S1'].index.tolist() + ['nas_steatosis_ordinal']].melt(id_vars=['nas_steatosis_ordinal'], value_name='Log2 intensity', var_name='Protein ID').dropna(subset=['nas_steatosis_ordinal', 'Log2 intensity'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "source_data_figs7a['Gene name']=source_data_figs7a['Protein ID'].map(IDmapping_UniprotID_to_Genename)\n",
    "source_data_figs7b['Gene name']=source_data_figs7b['Protein ID'].map(IDmapping_UniprotID_to_Genename)\n",
    "source_data_figs7c['Gene name']=source_data_figs7c['Protein ID'].map(IDmapping_UniprotID_to_Genename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SOURCE_DATA_FIGS7 = os.path.join('data/processed', 'source_data_figS7a-c.xlsx')\n",
    "\n",
    "with pd.ExcelWriter(SOURCE_DATA_FIGS7) as writer:\n",
    "    source_data_figs7a.to_excel(writer, sheet_name='FigS7a')\n",
    "    source_data_figs7b.to_excel(writer, sheet_name='FigS7b')\n",
    "    source_data_figs7c.to_excel(writer, sheet_name='FigS7c')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "hist_scores = {\"targets\":{'F2':'kleiner', 'I2':'nas_inflam', 'S1':'nas_steatosis_ordinal'}, \n",
    "               \"labels\":{'F2':'kleiner score', 'I2':'inflammatory score', 'S1':'steatosis score'},\n",
    "              \"colors\":{'F2':'mediumblue', 'I2':'darkred', 'S1':'gray'}}\n",
    "for test_case in ['F2', 'I2', 'S1']:\n",
    "    for marker in mrmr_markers[test_case].index:\n",
    "    #for marker in markers_biorxiv[test_case].index:\n",
    "        fig=plt.subplots(figsize=(2.5, 2.5))\n",
    "        x=hist_scores['targets'][test_case]\n",
    "        y=marker\n",
    "        sns.boxplot(x=x,y=y, data=df3, \n",
    "                    color='white', width=0.6, medianprops=dict(color=\"white\", alpha=1))\n",
    "        sns.stripplot(x=x, y=y, data=df3, color=hist_scores['colors'][test_case], size=3)\n",
    "        gene_name = IDmapping_UniprotID_to_Genename[marker]\n",
    "        plt.title(gene_name, fontsize=14)\n",
    "        plt.xticks(fontsize=12)\n",
    "        plt.xlabel(hist_scores['labels'][test_case], fontsize=12)\n",
    "        plt.yticks(fontsize=12)\n",
    "        plt.ylabel('MS signal [Log2]', fontsize=12)\n",
    "        plt.savefig('figures/markers/{}_{}.pdf'.format(test_case, gene_name), dpi=120, bbox_inches='tight')\n",
    "        print('{}-{}-distribution{}'.format(x, y, df3[[x, y]].dropna().value_counts(x).sort_index()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Magnitude of the changes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df1=data_proteomics.copy()\n",
    "df2=data_cli.copy()\n",
    "ep_cols = ['kleiner', 'nas_inflam', 'nas_steatosis_ordinal']\n",
    "df3 = df1.join(df2[ep_cols])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_k = df3[mrmr_markers['F2'].index.tolist() + ['kleiner']].groupby('kleiner').median().T\n",
    "df_i = df3[mrmr_markers['I2'].index.tolist() + ['nas_inflam']].groupby('nas_inflam').median().T\n",
    "df_s = df3[mrmr_markers['S1'].index.tolist() + ['nas_steatosis_ordinal']].groupby('nas_steatosis_ordinal').median().T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for stage in [2, 3, 4]:\n",
    "    col_name = 'F{}/F0'.format(stage)\n",
    "    func = lambda x: 2**(x[stage]-x[0])\n",
    "    df_k[col_name]=df_k.apply(func, axis=1)\n",
    "    df_k['Gene name']=df_k.index.map(IDmapping_UniprotID_to_Genename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for stage in [2, 3, 4, 5]:\n",
    "    col_name = 'I{}/I0'.format(stage)\n",
    "    func = lambda x: 2**(x[stage]-x[0])\n",
    "    df_i[col_name]=df_i.apply(func, axis=1)\n",
    "    df_i['Gene name']=df_i.index.map(IDmapping_UniprotID_to_Genename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for stage in [1, 2, 3]:\n",
    "    col_name = 'S{}/S0'.format(stage)\n",
    "    func = lambda x: 2**(x[stage]-x[0])\n",
    "    df_s[col_name]=df_s.apply(func, axis=1)\n",
    "    df_s['Gene name']=df_s.index.map(IDmapping_UniprotID_to_Genename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Write results to excel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "FILE_RESULTS = os.path.join(TABLEFOLDER, 'Magnitudes_marker_panels_mrmr.xlsx')\n",
    "\n",
    "with pd.ExcelWriter(FILE_RESULTS) as writer:\n",
    "    df_k.to_excel(writer, sheet_name='F2')\n",
    "    df_i.to_excel(writer, sheet_name='I2')\n",
    "    df_s.to_excel(writer, sheet_name='S1')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Old marker panels in the bioRxiv version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "markers_bioRxiv = pd.read_csv('data/processed/Old_marker_panels_bioRxiv.csv')\n",
    "markers_biorxiv={'F2':markers_bioRxiv['F2'], 'I2':markers_bioRxiv['I2'], 'S1':markers_bioRxiv['S1']}\n",
    "for key in markers_biorxiv:\n",
    "    markers_biorxiv[key]=key_ProteinID.loc[data_proteomics.columns.intersection(markers_biorxiv[key])]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "protein = 'P10643'\n",
    "in_both_old = set(kleiner_ge_2.index) & set(data_proteomics[protein].index)\n",
    "roc_auc_score(y_true=kleiner_ge_2.loc[in_both_old], y_score=data_proteomics.loc[in_both_old][protein])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "### Executor \n",
    "\n",
    "corresponds to main function in a script. Allows changes over all endpoints simultaneously. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "from src.cross_validation import MainExecutorCV\n",
    "cv_executor = MainExecutorCV(proteomics_data=data_proteomics, clinical_data=data_cli, demographics=demographics, clf_sklearn=clf_sklearn, cutoffs_clinic=cutoffs_clinic)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "1. Models based on clinical marker thresholds (Clinical Reference Models) - defined by `cutoffs_clinic`\n",
    "2. Additional Models based on clinical markers (having no standard cutoffs defined or if data dependent cutoff is wanted) - defined by list `additional_markers`\n",
    "    - As this depends on the endpoint, it is \n",
    "3. Proteomics Models based on protein intensities\n",
    "\n",
    "Performance depends on the number of available features (varies!) \n",
    "\n",
    "Result tables:\n",
    "- metrics are reported for the test set\n",
    "- `N_obs` is the number of patients both in the training and testing data set which is split into 80% training and 20% testing splits."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ADD_DEMOGRAPHICS = False\n",
    "INTERACTION_DEGREE = 1\n",
    "VERBOSE = False"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "### Fibrosis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "#### F2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "kleiner_ge_2.value_counts(dropna=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "cutoffs_f2 = cutoffs_clinic['F2'].dropna().to_dict()\n",
    "cutoffs_f2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "y = kleiner_ge_2.astype(int)\n",
    "y.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "proteins_selected = mrmr_markers['F2']\n",
    "#proteins_selected = proteins_selected_f2\n",
    "f2_results, f2_auc_scores, f2_prc_scores = cv_executor.run_evaluation(y=kleiner_ge_2,\n",
    "                                                       endpoint='F2',\n",
    "                                                       additional_markers=['forns', 'p3np'],\n",
    "                                                       proteins_selected=proteins_selected,\n",
    "                                                       add_demographics=ADD_DEMOGRAPHICS,\n",
    "                                                       interactions_degree=1,\n",
    "                                                       cv=cv_train_test_indices,\n",
    "                                                       verbose=VERBOSE)\n",
    "result_table_f2 = _get_cv_means(f2_results).sort_values(('f1', 'mean'), ascending = False)\n",
    "result_table_f2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert result_table_f2.loc['F2_prot_Logistic', ('roc_auc', 'mean')].round(5) - 0.90237 < 0.0001 , \"Final results not reproduced.\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "#### F3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "kleiner_ge_3.value_counts(dropna=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "cutoffs_f3 = cutoffs_clinic['F3'].dropna().to_dict()\n",
    "cutoffs_f3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "Two clinical markers have no cutoff defined in the literature. Therefore we have to learn these"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "proteins_selected = mrmr_markers['F3']\n",
    "#proteins_selected = proteins_selected_f3\n",
    "f3_results, f3_auc_scores, f3_prc_scores = cv_executor.run_evaluation(y=kleiner_ge_3, endpoint='F3', \n",
    "                                                       additional_markers=['p3np'], \n",
    "                                                       proteins_selected=proteins_selected,\n",
    "                                                       add_demographics=ADD_DEMOGRAPHICS,\n",
    "                                                       interactions_degree=1,\n",
    "                                                       cv=cv_train_test_indices,\n",
    "                                                       verbose=VERBOSE)\n",
    "result_table_f3 = _get_cv_means(f3_results).sort_values(('f1', 'mean'), ascending = False)\n",
    "result_table_f3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "Using one of the models (or the ensemble), one could expect some predictions of fibrosis patients in the untested healthy patient (hp) cohort."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "### Inflamation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "inflamation_ge_2.value_counts(dropna=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "cutoffs_i2 = cutoffs_clinic['I2'].dropna().to_dict()\n",
    "cutoffs_i2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "proteins_selected = mrmr_markers['I2']\n",
    "#proteins_selected = proteins_selected_I2\n",
    "#proteins_selected = markers_biorxiv['I2']\n",
    "i2_results, i2_auc_scores, i2_prc_scores = cv_executor.run_evaluation(y=inflamation_ge_2, endpoint='I2', \n",
    "                                                       additional_markers=['m30', 'm65', 'alt', 'ast', 'm30m65_ratio'], \n",
    "                                                       proteins_selected=proteins_selected,\n",
    "                                                       add_demographics=ADD_DEMOGRAPHICS,\n",
    "                                                       interactions_degree=1,\n",
    "                                                       cv=cv_train_test_indices,\n",
    "                                                       verbose=VERBOSE)\n",
    "result_table_i2 = _get_cv_means(i2_results).sort_values(('f1', 'mean'), ascending = False)\n",
    "result_table_i2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "### Steatosis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "steatosis_ge_1.value_counts(dropna=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "cutoffs_s1 = cutoffs_clinic['S1'].dropna().to_dict()\n",
    "cutoffs_s1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "y = steatosis_ge_1.astype(int)\n",
    "y.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "proteins_selected = mrmr_markers['S1']\n",
    "#proteins_selected = proteins_selected_s1\n",
    "#proteins_selected = markers_biorxiv['S1']\n",
    "s1_results, s1_auc_scores, s1_prc_scores = cv_executor.run_evaluation(y=steatosis_ge_1, endpoint='S1', \n",
    "                                                       additional_markers=[], \n",
    "                                                       proteins_selected=proteins_selected,\n",
    "                                                       add_demographics=ADD_DEMOGRAPHICS,\n",
    "                                                       interactions_degree=1,\n",
    "                                                       cv=cv_train_test_indices,\n",
    "                                                       verbose=VERBOSE)\n",
    "result_table_s1 = _get_cv_means(s1_results).sort_values(('f1', 'mean'), ascending = False)\n",
    "result_table_s1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "### Write results to Excel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "FILE_RESULTS = os.path.join(TABLEFOLDER, 'CV_results_proteomics.xlsx')\n",
    "\n",
    "with pd.ExcelWriter(FILE_RESULTS) as writer:\n",
    "    result_table_f2.to_excel(writer, sheet_name='F2_featureOptim')\n",
    "    result_table_f3.to_excel(writer, sheet_name='F3_featureOptim')\n",
    "    result_table_i2.to_excel(writer, sheet_name='I2_featureOptim')\n",
    "    result_table_s1.to_excel(writer, sheet_name='S1_featureOptim')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "## Plot Results of Cross validation for three endpoints (F2, I2, S1)\n",
    "\n",
    "- create [enumeration of subplots](https://stackoverflow.com/a/25544329/9684872) starting at a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_names = pd.read_csv(os.path.join(FOLDER_DATA_RAW, 'naming_scheme.csv'), index_col='name_in_clinical_data')\n",
    "map_names = map_names['name_in_plot'].to_dict()\n",
    "map_names"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Custom function to transform index names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "display(result_table_s1)\n",
    "\n",
    "def _process_names(index, map_names=map_names):\n",
    "    \"\"\"Helper function for custom labeling of models.\n",
    "    This function is specific to any dataset and has to be rewritten.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    index: pandas.Index\n",
    "        Index to transform. Index names are composite word\n",
    "        combined with '_' here.\n",
    "    map_names: dict\n",
    "        Mapping of names to apply to words.\n",
    "    \"\"\"\n",
    "    names = list(index)\n",
    "    names = [x.split('_') for x in names]\n",
    "    endpoint = names[0][0]\n",
    "    \n",
    "    def _process_index_names(_l:list):\n",
    "        REMOVE = 'marker'\n",
    "        if REMOVE in _l:\n",
    "            _l.remove(REMOVE)\n",
    "        _l = [word if not word in map_names else map_names[word] for word in _l]\n",
    "        CHANGE = {'Logistic': 'Model',\n",
    "                  'prot': 'Proteomics'}\n",
    "        _l = [word if not word in CHANGE else CHANGE[word] for word in _l]\n",
    "        if CHANGE['Logistic'] not in _l:\n",
    "            _l.append('Test')\n",
    "        return _l\n",
    "    \n",
    "    for _l in names: assert endpoint == _l[0] , f\"Mixed endpoints: {endpoint} and {_l[0]}\"\n",
    "    names = [\" \".join(_process_index_names(_l[1:])) for _l in names]\n",
    "    return names\n",
    "\n",
    "_process_names(result_table_s1.index)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "### Performance Plots based on results DataFrame for a endpoint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "from src.plots import plot_performance\n",
    "fig, ax = plt.subplots(figsize=(10,10))\n",
    "plot_performance(ax, result_table_s1, 'balanced_accuracy', 'Steatosis', _process_index=_process_names)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Build bar plots with individual dots "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from src.plots import plot_performance_adddots\n",
    "fig, ax = plt.subplots(figsize=(10,10))\n",
    "plot_performance_adddots(ax, result_table_f2, f2_results, 'balanced_accuracy', 'Fibrosis', _process_index=_process_names)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "### AUC-ROC Curves based on CV result for an endpoint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "from src.plots import plot_roc_curve\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8, 8))      \n",
    "        \n",
    "plot_roc_curve(ax, roc_curve_results['F2_Logistic'], 'TARGET')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Precision Recall Curve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from src.plots import plot_prc_curve\n",
    "fig, ax = plt.subplots(figsize=(8, 8))\n",
    "plot_prc_curve(ax, precision_recall_results['F2_Logistic'], 'TARGET')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def collect_source_data(cross_validation_results, test_case):\n",
    "    df=cross_validation_results.copy()\n",
    "    df_source_data = pd.DataFrame.from_dict(df).loc[test_case].to_frame().reset_index()\n",
    "    df_source_data_split = pd.DataFrame(df_source_data[test_case].tolist(), columns=['CV_' + str(i) for i in range(50)])\n",
    "    df_source_data_combined=pd.concat([df_source_data[['index']], df_source_data_split], axis=1)\n",
    "    return df_source_data_combined"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "source_data_fig4a = pd.DataFrame(pd.DataFrame.from_dict(f2_auc_scores)['F2_prot_Logistic'].to_frame()['F2_prot_Logistic'].tolist(), columns=['fpr', 'tpr', 'thresholds'])\n",
    "source_data_fig4b = pd.DataFrame(pd.DataFrame.from_dict(i2_auc_scores)['I2_prot_Logistic'].to_frame()['I2_prot_Logistic'].tolist(), columns=['fpr', 'tpr', 'thresholds'])\n",
    "source_data_fig4c = pd.DataFrame(pd.DataFrame.from_dict(s1_auc_scores)['S1_prot_Logistic'].to_frame()['S1_prot_Logistic'].tolist(), columns=['fpr', 'tpr', 'thresholds'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SOURCE_DATA_FIG4 = os.path.join('data/processed/', 'source_data_fig4.xlsx')\n",
    "\n",
    "with pd.ExcelWriter(SOURCE_DATA_FIG4) as writer:\n",
    "    for score in ['f1', 'balanced_accuracy']:\n",
    "        collect_source_data(f2_results, score).to_excel(writer, sheet_name='F2_{}'.format(score))\n",
    "        collect_source_data(i2_results, score).to_excel(writer, sheet_name='I2_{}'.format(score))\n",
    "        collect_source_data(s1_results, score).to_excel(writer, sheet_name='S1_{}'.format(score))\n",
    "        source_data_fig4a.to_excel(writer, sheet_name='F2_AUC')\n",
    "        source_data_fig4b.to_excel(writer, sheet_name='I2_AUC')\n",
    "        source_data_fig4c.to_excel(writer, sheet_name='S1_AUC')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "### Build final figure for publication"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "import string\n",
    "fig, axs = plt.subplots(3,3,figsize=(20,20))\n",
    "\n",
    "n=0\n",
    "result_tuples = [\n",
    " (result_table_f2, f2_results, f2_auc_scores, 'Fibrosis F2-F4', 'F2_prot_Logistic'), \n",
    " (result_table_i2, i2_results, i2_auc_scores, 'NAS Inflamation $\\geq 2$', 'I2_prot_Logistic'), \n",
    " (result_table_s1, s1_results, s1_auc_scores, 'NAS Steatosis $\\geq 5$%', 'S1_prot_Logistic'), \n",
    "    \n",
    "]\n",
    "\n",
    "for col, (result_table, result_data, result_auc_scores, endpoint_title, auc_model_name) in enumerate(result_tuples):\n",
    "    \n",
    "    ax = axs[0,col]\n",
    "    plot_roc_curve(ax, result_auc_scores[auc_model_name], endpoint_title)\n",
    "    _ = ax.text(-0.5, 1.1, f\"{string.ascii_lowercase[n]})\", transform=ax.transAxes, \n",
    "                size=20, weight='bold')\n",
    "    n+=1\n",
    "    \n",
    "    ax = axs[1,col]\n",
    "    plot_performance_adddots(ax, result=result_table, result_data=result_data, metric='f1', title=endpoint_title,  _process_index=_process_names)   \n",
    "    _ = ax.text(-0.5, 1.1, f\"{string.ascii_lowercase[n]})\", transform=ax.transAxes, \n",
    "                size=20, weight='bold')\n",
    "    n+=1\n",
    "    \n",
    "    ax = axs[2,col]\n",
    "    _ = ax.text(-0.5, 1.1, f\"{string.ascii_lowercase[n]})\", transform=ax.transAxes, \n",
    "                size=20, weight='bold')\n",
    "    plot_performance_adddots(ax, result=result_table, result_data=result_data, metric='balanced_accuracy', title=endpoint_title, _process_index=_process_names)   \n",
    "\n",
    "    n+=1\n",
    "\n",
    "fig.tight_layout()\n",
    "fig.savefig(FIGURE_FOLDER / 'Figure_6_new.pdf', dpi=120, pad_inches=0.1, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Final model\n",
    "\n",
    "1. Either pick one of the models run during CV\n",
    "2. Aggregate metrics over all CV runs for obs in test set (~ mean of CV results)\n",
    "3. Perform new train/test split.\n",
    "\n",
    "\n",
    "Steps to implement\n",
    "\n",
    "1.  Select model with median aucroc performance\n",
    "2.  Report summary statistics (mean, median, min, max)\n",
    "3.  [DeLong](https://github.com/llniu/roc_comparison)\n",
    "4.  Target Scores for three endpoints of final prediction model\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Check Cross Validation results (for comparison)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Look at descriptive statistics of CV\n",
    "\n",
    "- depending on the split the performance varies. The `min`, `max`, `mean`, etc. are given per model, therefore it's not the given split. Selecting a split of the data which supports one's conclusion can be misleading."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = {\n",
    " \"F2\":f2_results,\n",
    " \"F3\":f3_results,\n",
    " \"I2\":i2_results,\n",
    " \"S1\":s1_results\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_combinded = {}\n",
    "for _results in results.values(): results_combinded.update(_results)\n",
    "results_combinded.keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def show_summary(results_dict, metric='f1', sort=True, save=False):\n",
    "    _df = pd.DataFrame(results_dict).loc[metric].apply(pd.Series).T.describe()\n",
    "    _df.index.name = metric\n",
    "    if sort:\n",
    "        _df = _df.sort_values(by='mean', axis=1, ascending=False)\n",
    "    if save:\n",
    "        fname = \"cv_stats_{}{}.xlsx\".format(metric, '_sorted' if save else '')\n",
    "        fname = os.path.join(TABLEFOLDER, fname)\n",
    "        _df.to_excel(fname)\n",
    "        print(f'Saved Table to: {fname}')\n",
    "    return _df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "out= widgets.interact(show_summary, metric=scoring, results_dict=widgets.fixed(results_combinded))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Select a model from the CV run (and compare it to others)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "l_models = list(results_combinded.keys())\n",
    "ref_model=l_models[14] # 'F2_prot_Logistic'\n",
    "metric=scoring[2]      # 'f1'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compare_models(ref_model:str, metric:str, data:dict, fct=np.median):\n",
    "    \"\"\"Select comparison. The first column is the referenc model.\n",
    "    If the metric of the summary statistic is present several times in the model, this value is returned.\"\"\"\n",
    "    ref_model_metric_values = np.array(data[ref_model][metric])\n",
    "    ref_model_metric = fct(ref_model_metric_values)\n",
    "    ref_model_metric_idx = (np.abs(ref_model_metric_values - ref_model_metric)).argmin() #ToDo: return both closest values for median.\n",
    "    ref_model_metric = ref_model_metric_values[ref_model_metric_idx]\n",
    "    matches = [index for index, item in enumerate(ref_model_metric_values) if item == ref_model_metric]\n",
    "\n",
    "    _results = pd.DataFrame(ref_model_metric_values[matches], columns=[ref_model], index=matches)\n",
    "    _results.index.name = 'run'\n",
    "    \n",
    "    _selected_results = {}\n",
    "    \n",
    "    \n",
    "    for _model, _result in data.items():\n",
    "        if not _model == ref_model:\n",
    "            _selected_results[_model] = [_result[metric][item] for item in matches]\n",
    "    #sorting of values over last result\n",
    "    _other = pd.DataFrame(_selected_results, index=matches)\n",
    "    _other = _other.sort_values(by=matches[-1], axis=1, ascending=False)\n",
    "    return _results.join(_other)\n",
    "\n",
    "# compare_models(ref_model=ref_model, metric='precision', data=f2_results, fct=lambda x: np.quantile(x, q=0.6))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "protein_model_name = {endpoint: f'{endpoint}_prot_Logistic'for endpoint in end_points}\n",
    "metrics_np_fct = {'median': np.median, '3rd quintile': lambda x: np.quantile(x, q=0.6), 'mean': np.mean, 'max': np.max, 'min': np.min }\n",
    "\n",
    "#needs global dictionaries: protein_model_name, metrics_np_fct, results\n",
    "def _caller_comp(metric, endpoint, selector):\n",
    "    \"\"\"Helper function to use with ipykernel\"\"\"\n",
    "    df = compare_models(ref_model=protein_model_name[endpoint], metric=metric, data=results[endpoint], fct=metrics_np_fct[selector])\n",
    "    df.columns.name = selector\n",
    "    return df\n",
    "# _caller_comp('f1', 'I2', 'median')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "out2= widgets.interact(_caller_comp, metric=scoring, endpoint=results.keys(), selector=metrics_np_fct.keys())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Train a new final model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Adjust `run_cv_binary`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train, X_test, y_train, y_test = sklearn.model_selection.train_test_split(data_cli_missing_table, stratifier, test_size=0.2, stratify=stratifier, random_state=42)\n",
    "logging.info(f\"N train: {len(y_train)}\")\n",
    "logging.info(f\"N test: {len(y_test)}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "combine_value_counts(pd.DataFrame({'train': y_train, 'test': y_test}))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_test_split_final = [(X_train.index, X_test.index)] # list of (train, test) indices. Does not throw a nice error message otherwise."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Adapted the executor to take as argument a different function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "arguments = {}\n",
    "arguments['F2'] = {'target': kleiner_ge_2, 'proteins': mrmr_markers['F2'], 'add_markers': ['forns', 'p3np']} \n",
    "arguments['F3'] = {'target': kleiner_ge_3, 'proteins': mrmr_markers['F3'], 'add_markers': ['p3np']}\n",
    "arguments['I2'] = {'target': inflamation_ge_2, 'proteins': mrmr_markers['I2'], 'add_markers': ['m30', 'm65', 'alt', 'ast', 'm30m65_ratio']}\n",
    "arguments['S1'] = {'target': steatosis_ge_1, 'proteins': mrmr_markers['S1'], 'add_markers': []}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ADD_DEMOGRAPHICS = False\n",
    "FOLDER_FINAL_SCORES = 'final_model_scores'\n",
    "\n",
    "import importlib; importlib.reload(src.cross_validation)\n",
    "from src.cross_validation import run_cv_binary\n",
    "prediction_folder = os.path.join(TABLEFOLDER, FOLDER_FINAL_SCORES)\n",
    "\n",
    "\n",
    "from functools import partial\n",
    "run_cv_binary_final = partial(run_cv_binary, save_predictions=True, folder=prediction_folder)\n",
    "\n",
    "\n",
    "results_final = {}\n",
    "\n",
    "for endpoint, args in arguments.items():\n",
    "    results_final[endpoint] = cv_executor.run_evaluation(y=args['target'],\n",
    "                                                       endpoint=endpoint,\n",
    "                                                       additional_markers=args['add_markers'],\n",
    "                                                       proteins_selected=args['proteins'],\n",
    "                                                       evaluator_fct=run_cv_binary_final,\n",
    "                                                       add_demographics=ADD_DEMOGRAPHICS,\n",
    "                                                       interactions_degree=1,\n",
    "                                                       cv=train_test_split_final,\n",
    "                                                       verbose=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def display_result(endpoint):\n",
    "    return pd.DataFrame(results_final[endpoint][0]).applymap(lambda x: x[0]).T.sort_values('f1', ascending=False)\n",
    "\n",
    "\n",
    "# pd.DataFrame(results_final['F2'][0]).T.sort_values('f1', ascending=False)\n",
    "out3= widgets.interact(display_result, endpoint=results_final.keys())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_final_df = {endpoint: display_result(endpoint) for endpoint in results_final.keys()}\n",
    "\n",
    "RESULTS_FINAL_MODEL = os.path.join(TABLEFOLDER ,'final_model_results.xlsx')\n",
    "\n",
    "with pd.ExcelWriter(RESULTS_FINAL_MODEL) as writer:\n",
    "    for sheet_name, _df in results_final_df.items():\n",
    "        _df.to_excel(writer, sheet_name=sheet_name)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### DeLong-Test on final model\n",
    "\n",
    "First a check of the [implementation](https://github.com/yandexdataschool/roc_comparison) by using a toy example derived in detail in Rachel Draelos' [blog-post](https://glassboxmedicine.com/2020/02/04/comparing-aucs-of-machine-learning-models-with-delongs-test/), where she also references the original DeLong paper from 1988 and a paper on a fast implementation described by Xu Sun and Weichao Xu from 2014."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from roc_comparison import compare_auc_delong_xu\n",
    "\n",
    "ground_trouth = np.array([0,0,1,1,1])\n",
    "pred_model_a  = np.array([0.1,0.2,0.6,0.7,0.8])\n",
    "pred_model_b  = np.array([0.3,0.6,0.2,0.7,0.9])\n",
    "\n",
    "log10_pvalue = compare_auc_delong_xu.delong_roc_test(ground_truth=ground_trouth,\n",
    "                                      predictions_one=pred_model_a,\n",
    "                                      predictions_two=pred_model_b\n",
    "                                 )\n",
    "assert np.round(10**log10_pvalue[0][0], 4) ==  0.3173 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The implementation currently returns an array of an array with a single float. This could be change to an normal array or a plain floating point number:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "log10_pvalue[0][0], log10_pvalue[0] # one of the two"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get list of dumped results from models."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "folder_final_scores = os.path.join(TABLEFOLDER, FOLDER_FINAL_SCORES)\n",
    "l_scores = [_csv for _csv in os.listdir(folder_final_scores) if 'Logistic.csv' in _csv]\n",
    "l_scores"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare all models between all endpoints and highlight the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from src.delong import calc_p_value_delong_xu\n",
    "\n",
    "model_1 = l_scores[0]\n",
    "model_2 = l_scores[5]\n",
    "\n",
    "print(\n",
    "f\"Delong-Test p-value between scores of model {model_1.split('.',1)[0]} and {model_2.split('.')[0]}: \"\n",
    "f\"{calc_p_value_delong_xu(model_1=model_1, model_2=model_2, folder_dumps=folder_final_scores, verbose=True):.4f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "auc_comp_dict = {}\n",
    "\n",
    "def highlight_significant(val):\n",
    "    \"\"\"\n",
    "    Takes a scalar and returns a string with\n",
    "    the css property `'color: red'` for floats smaller\n",
    "    equal to 0.05\n",
    "    \"\"\"\n",
    "    color = 'red' if val <= 0.05 else 'black'\n",
    "    return 'color: %s' % color\n",
    "\n",
    "for endpoint in end_points:\n",
    "    # models = [x for x in l_scores if endpoint in x]\n",
    "    # model_names = [\" \".join(x.split('.csv')[0].split('_')) for x in models]\n",
    "    model_names = pd.DataFrame(results_final[endpoint][0]).applymap(lambda x: x[0]).T.sort_values('f1', ascending=False).index.to_list()\n",
    "    model_names = [model for model in model_names if 'Logistic' in model]\n",
    "    models = [f'{name}.csv' for name in model_names]\n",
    "    _df = pd.DataFrame(0, index=models, columns=models)\n",
    "    \n",
    "    for i, model_1 in enumerate(models):\n",
    "        for model_2 in models[i:]:\n",
    "            if model_1 == model_2:\n",
    "                _df.loc[model_1, model_2] = 1.0\n",
    "            else:\n",
    "                _auc_p_value = calc_p_value_delong_xu(model_1, model_2, folder_dumps=folder_final_scores)\n",
    "                _df.loc[model_1, model_2] = _auc_p_value\n",
    "                _df.loc[model_2, model_1] = _auc_p_value\n",
    "    model_names = [\" \".join(x.split('.csv')[0].split('_')) for x in model_names]\n",
    "    _df.columns = model_names\n",
    "    _df.index   = model_names\n",
    "    _df =  _df.style.applymap(highlight_significant)\n",
    "    display(_df)\n",
    "    auc_comp_dict[endpoint] = _df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save results for supplementary materials:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "FILE_DELONG = Path(TABLEFOLDER) / 'compare_delong.xlsx'\n",
    "\n",
    "with pd.ExcelWriter(FILE_DELONG) as writer:\n",
    "    auc_comp_dict['F2'].to_excel(writer, sheet_name='F2')\n",
    "    auc_comp_dict['F3'].to_excel(writer, sheet_name='F3')\n",
    "    auc_comp_dict['I2'].to_excel(writer, sheet_name='I2')\n",
    "    auc_comp_dict['S1'].to_excel(writer, sheet_name='S1')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Final Model evaluation\n",
    "\n",
    "The healthy patients in the study - have been selected based on age and gender range of the non-healthy patients. These have an `NaN` assigned on the three endpoints.\n",
    "Another group has been previously excluded from the analysis and can be used here as an pseudo external dataset \n",
    "\n",
    "> Note: Pseudo as for now all proteomics pre-processing is done on the entire dataset. Either the procedure is done individuelly for each subset of patients, \n",
    "or the mean and std deviation from the training dataset are used to sample random values for missing protein intensities on the test data. The standardisation of the \n",
    "proteins intensities after imputation should ideally also be based on the statistics from the training data (which is assumed to the \"global\" value\").\n",
    "\n",
    "Patiens with an assigned fibrosis score of `0.5` are known to be heavy drinkers without being diagnosed with a fibrosis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_cli_full = pd.read_csv(f_data_clinic, index_col=COL_ID)\n",
    "# previous selection: data_cli[data_cli['kleiner']!=0.5] \n",
    "all_kleiner_score = data_cli_full.kleiner.value_counts(dropna=False)\n",
    "all_kleiner_score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#n_healty, n_at_risk = all_kleiner_score.loc[[np.nan, 0.5]]\n",
    "#One patient in the disease cohort had a 'NaN' kleiner score\n",
    "n_healthy = data_cli[data_cli['group']=='HP'].shape[0]\n",
    "n_at_risk = data_cli_full[data_cli_full['kleiner']==0.5].shape[0]\n",
    "print(f\"N Healthy (selected to match ill patients: {n_healthy}\")\n",
    "print(f\"N Unhealthy behaviour, but not sick: {n_at_risk}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# from src.pandas import combine_value_counts\n",
    "combine_value_counts(data_cli_full[TARGETS], dropna=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "combine_value_counts(data_cli_full.loc[data_cli_full.kleiner == 0.5, TARGETS], dropna=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_cli_full.loc[data_cli_full.kleiner == 0.5, TARGETS]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Load final models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from joblib import load\n",
    "endpoints = ['F2', 'F3', 'I2', 'S1']\n",
    "fname_final_model = '{}_prot_Logistic.joblib'\n",
    "final_prot_model = {}\n",
    "for endpoint in endpoints :\n",
    "    _fname = os.path.join('tables', FOLDER_FINAL_SCORES, fname_final_model.format(endpoint))\n",
    "    print(f\"Load model from : {_fname}\")\n",
    "    final_prot_model[endpoint] = load(_fname)\n",
    "    print(final_prot_model[endpoint].coef_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Protein by endpoint can be retrieved using the previously defined dictionary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from src.final_model import FinalPredictor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "final_predictor = FinalPredictor(data_clinic=data_cli_full, \n",
    "                                 data_proteomics=data_proteomics, \n",
    "                                 final_models=final_prot_model, \n",
    "                                 features_dict=arguments, \n",
    "                                 endpoints=endpoints)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Feature Importance of final models\n",
    "\n",
    "- how important are the single proteins, assessed by the model weights."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# protein panel. Some order?\n",
    "d_model_weights = {}\n",
    "d_model_intercepts = {}\n",
    "for _endpoint, _model in final_predictor.final_models.items():\n",
    "    d_model_weights[_endpoint] = dict(zip(final_predictor.features_dict[_endpoint]['proteins'].index.to_list(), _model.coef_[0]))\n",
    "    d_model_intercepts[_endpoint] = _model.intercept_[0]\n",
    "del _endpoint, _model\n",
    "\n",
    "model_weights = pd.DataFrame(d_model_weights)\n",
    "\n",
    "model_weights.index =  [f'{_index}_{_gene}' for _index, _gene \n",
    " in key_ProteinID.loc[model_weights.index].to_dict()['Gene names'].items()\n",
    "]\n",
    "\n",
    "model_weights = model_weights.append(pd.Series(d_model_intercepts, name='intercept'))\n",
    "del d_model_weights, d_model_intercepts\n",
    "model_weights.to_excel(Path(DATAFOLDER) / 'final_model_weights.xlsx') \n",
    "new_index = [i.split('_')[1] for i in model_weights.index[:-1]]\n",
    "new_index.append('intercept')\n",
    "model_weights.index=new_index\n",
    "model_weights=model_weights.sort_values(by='F2', ascending=False).iloc[::-1].dropna(how='all')\n",
    "model_weights_F2=model_weights[model_weights['F2'].notnull()][['F2']]\n",
    "model_weights_I2=model_weights[model_weights['I2'].notnull()][['I2']].sort_values(by='I2')\n",
    "model_weights_S1=model_weights[model_weights['S1'].notnull()][['S1']].sort_values(by='S1')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = model_weights[['F2', 'I2', 'S1']].plot(kind='barh', figsize=(3,8), color=['darkblue', 'darkred', 'gray'])\n",
    "plt.xticks(fontsize=14);\n",
    "plt.yticks(fontsize=12);\n",
    "plt.xlabel('final model weights', fontsize=14);\n",
    "ax.figure.savefig(FIGURE_FOLDER / 'final_model_weights.png', dpi=120, pad_inches=0.1, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = model_weights[['S1']].sort_values(by='S1').dropna().plot(kind='barh', figsize=(2, 7), color='gray')\n",
    "plt.xticks(fontsize=14);\n",
    "plt.yticks(fontsize=12);\n",
    "plt.xlabel('final model weights', fontsize=14);\n",
    "ax.figure.savefig(FIGURE_FOLDER / 'final_model_weights_S1.png', dpi=120, pad_inches=0.1, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Risk scores from the final predictor on all patients in GALA-ALD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ald_index = annotation_file[(annotation_file['Sample type']=='Plasma') & (annotation_file['Group2']=='ALD')]['Sample ID']\n",
    "risk_scores_ald = final_predictor.predict_score(indices=ald_index)\n",
    "risk_scores_ald['Participant ID']=risk_scores_ald.index.map(dict(zip(annotation_file['Sample ID'], annotation_file['Participant ID'])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert risk_scores_ald.loc['Plate1_A4', 'F2'].round(4) == 0.6181, 'Fail to reproduce final risk score'\n",
    "print('Final results reproduced')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "risk_scores_ald['Participant ID']=risk_scores_ald['Participant ID'].str.replace('SIPHON', 'ALD')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "risk_scores_ald=risk_scores_ald.set_index('Participant ID')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "risk_scores_ald.to_csv('tables/final_model_scores/tables_to_Maja/risk_scores_10.16.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### At risk population\n",
    "\n",
    "How many patients at risk have been predicted to have a certain disease stage?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "at_risk_score = final_predictor.predict_score(indices=data_cli_full.kleiner[data_cli_full.kleiner == 0.5].index)\n",
    "at_risk_pred = final_predictor.predict(indices=data_cli_full.kleiner[data_cli_full.kleiner == 0.5].index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# minimal check to catch big errors\n",
    "from numpy.testing import assert_array_almost_equal\n",
    "assert_array_almost_equal(x=at_risk_score.loc['Plate4_D9'].values,  \n",
    "                          y=[0.58544614, 0.0241852 , 0.44440634, 0.66941346], decimal=4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def highlight_surpassed_score(val):\n",
    "    \"\"\"\n",
    "    Takes a scalar and returns a string with\n",
    "    the css property `'color: red'` for floats smaller\n",
    "    equal to 0.05\n",
    "    \"\"\"\n",
    "    color = 'yellow' if val >= 0.5 else None\n",
    "    return f'background-color: {color}'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "at_risk_score = at_risk_score.sort_values(by=list(at_risk_score.columns), ascending=False).style.applymap(highlight_surpassed_score)\n",
    "at_risk_pred = at_risk_pred.loc[at_risk_score.index]\n",
    "at_risk_pred.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def build_results_table(predictions:pd.DataFrame, dropna=True):\n",
    "    \"\"\"Convert a table of predicitons (for several columns) into a result table.\"\"\"\n",
    "    predictions_tab = combine_value_counts(predictions)\n",
    "    predictions_tab = predictions_tab.join(combine_value_counts(predictions, dropna=dropna) / len(predictions), rsuffix='_freq')\n",
    "    predictions_tab.loc['Total'] = predictions_tab.sum()\n",
    "    return predictions_tab.convert_dtypes()\n",
    "\n",
    "at_risk_results_tab  = build_results_table(at_risk_pred)\n",
    "at_risk_results_tab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ToDo: move to src (replace index building in cross_validation._get_cv_means)\n",
    "def build_two_level_index(initial_columns, present_key='freq', added_key='prop'):\n",
    "    \"\"\"Build a custom multi-index object for data. The initial columns to build the \n",
    "    data object are passed along the key describing the intial data and added information by \n",
    "    a second key. Could be generalized to work with any number of keys.\"\"\"\n",
    "    column_map = []\n",
    "    for x in initial_columns:\n",
    "        column_map += [x, x + '_freq']\n",
    "\n",
    "    levels = [initial_columns, [present_key, added_key]]\n",
    "    multi_index = pd.MultiIndex.from_product(\n",
    "        levels, names=['variable', 'statistics'])\n",
    "    return column_map, multi_index\n",
    "\n",
    "column_map, multi_index = build_two_level_index(initial_columns=at_risk_pred.columns)\n",
    "at_risk_results_tab = at_risk_results_tab[column_map]\n",
    "at_risk_results_tab.columns = multi_index\n",
    "at_risk_results_tab"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "### Healthy\n",
    "\n",
    "How many healthy patients would be predicted to have fibrosis by the final model?\n",
    "\n",
    "Here we check for formally \"healthy\" patients (they have not been diagnosed at the initial time of the data collection) and see who would be predicted to \n",
    "have a fibrosis score of 2 and above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will predict for patients without a fibriosis score. We load the **final** proteomics model from the previous step, lookup the blood plasma protein intensities for these patients on the selected proteins for the classification model, and finally predict their outcome."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_cli.fibrosis_class.value_counts(dropna=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_cli[data_cli['group']=='HP'].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "healthy_cohort_mask = data_cli[data_cli['group']=='HP']\n",
    "healthy_cohort_indices = healthy_cohort_mask.index\n",
    "healthy_pred = final_predictor.predict(indices=healthy_cohort_indices)\n",
    "healthy_risk_score = final_predictor.predict_score(indices=healthy_cohort_indices)\n",
    "healthy_risk_score = healthy_risk_score.sort_values(by=list(healthy_risk_score.columns), ascending=False).style.applymap(highlight_surpassed_score)\n",
    "healthy_pred = healthy_pred.loc[healthy_risk_score.index]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "healthy_results_tab  = build_results_table(healthy_pred)\n",
    "column_map, multi_index = build_two_level_index(initial_columns=healthy_pred.columns)\n",
    "healthy_results_tab = healthy_results_tab[column_map]\n",
    "healthy_results_tab.columns = multi_index\n",
    "healthy_results_tab"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### build figure for manuscript"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "negative_assignments = {'healthy_≥'+i:healthy_results_tab[i]['prop'].loc[0].round(2) for i in ['F2', 'F3', 'I2', 'S1']}\n",
    "for i in ['F2', 'F3']:\n",
    "    negative_assignments['ALD subset_≥'+i]=at_risk_results_tab[i]['prop'].loc[0].round(2)\n",
    "neg_ass = pd.DataFrame({'cohort_model':negative_assignments.keys(),'proportion of negative assignments':negative_assignments.values()})\n",
    "neg_ass['model']=neg_ass['cohort_model'].str.split('_').str[1]\n",
    "neg_ass['cohort'] = neg_ass['cohort_model'].str.split('_').str[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax =plt.subplots(figsize=(4,4))\n",
    "plt.bar(x=neg_ass['cohort_model'], height=neg_ass['proportion of negative assignments'], width=0.6)\n",
    "plt.xticks(ticks=range(len(neg_ass['cohort_model'])), labels=neg_ass['model'], fontsize=16);\n",
    "plt.ylabel('Proportion of\\nnegative assignments', fontsize=16)\n",
    "plt.yticks(fontsize=16);\n",
    "#plt.savefig('figures/negative_assignments.pdf', dpi=120, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Final model performance on test set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "targets_final_test = pd.DataFrame(targets_dict).loc[train_test_split_final[0][1]]\n",
    "pred_final_test = final_predictor.predict(indices=train_test_split_final[0][1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import confusion_matrix\n",
    "def confusion_matrix_scorer(y, y_pred):\n",
    "    cm = confusion_matrix(y, y_pred)\n",
    "    return {'tn': cm[0, 0], 'fp': cm[0, 1],\n",
    "              'fn': cm[1, 0], 'tp': cm[1, 1]}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_cases = ['F2','F3', 'I2', 'S1']\n",
    "cm_scores_final_test = {}\n",
    "for test_case in test_cases:\n",
    "    overlap_ids = set(pred_final_test[test_case].dropna().index) & set(targets_final_test[test_case].dropna().index)\n",
    "    cm_scores_final_test[test_case]=confusion_matrix_scorer(targets_final_test.loc[overlap_ids][test_case], pred_final_test.loc[overlap_ids][test_case]) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def performance_score_metric(confusion_matrix):\n",
    "    \"\"\"\n",
    "    confusion_matrix as a dictionary with keys of 'tn', 'fp', 'fn', 'tp'.\n",
    "    \"\"\"\n",
    "    dict_cm =confusion_matrix\n",
    "    tps = dict_cm['tp']\n",
    "    fns = dict_cm['fn']\n",
    "    fps = dict_cm['fp']\n",
    "    tns = dict_cm['tn']\n",
    "    PPV = tps/(tps + fps)\n",
    "    NPV = tns/(tns + fns)\n",
    "    recall = tps/(tps + fns)\n",
    "    specificity = tns/(tns + fps)\n",
    "    balanced_accuracy = (recall + specificity)/2\n",
    "    f1_score = 2*PPV*recall/(PPV + recall)\n",
    "    scores = {'PPV': PPV.round(2), 'NPV': NPV.round(2), 'recall':recall.round(2), \n",
    "              'specificity':specificity.round(2), 'F1 score':f1_score.round(2), 'balanced accuracy':balanced_accuracy.round(2)}\n",
    "    return(scores)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "performance_score_metric(cm_scores_final_test['F2'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Independent ALD validation cohort"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_cli_validation=pd.read_csv('data/raw/validation_clinical_data_aldupgrade.csv')\n",
    "Y_val=data_cli_validation[['Sample ID', 'kleiner', 'nas_steatosis_ordinal', 'nas_inflam']].set_index('Sample ID')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kleiner_ge_2_val     = create_dichotome(Y_val['kleiner'], 2)\n",
    "kleiner_ge_3_val     = create_dichotome(Y_val['kleiner'], 3)\n",
    "steatosis_ge_1_val   = create_dichotome(Y_val['nas_steatosis_ordinal'], 1)\n",
    "inflamation_ge_2_val = create_dichotome(Y_val['nas_inflam'], 2)\n",
    "\n",
    "end_points_val = ['F2', 'F3', 'S1', 'I2']\n",
    "dichotomies_val = [kleiner_ge_2_val, kleiner_ge_3_val, steatosis_ge_1_val, inflamation_ge_2_val]\n",
    "targets_dict_val  = {k: v for k, v in zip(end_points_val, dichotomies_val)}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "freq_targets_val = pd.DataFrame(\n",
    "    {'kleiner>=2': kleiner_ge_2_val.value_counts(dropna=False, sort=False),\n",
    "     'kleiner>=3': kleiner_ge_3_val.value_counts(dropna=False, sort=False),\n",
    "     'steatosis>=1' : steatosis_ge_1_val.value_counts(dropna=False, sort=False),\n",
    "     'inflamation>=2':inflamation_ge_2_val.value_counts(dropna=False, sort=False)\n",
    "    })\n",
    "freq_targets_val.loc['total'] = freq_targets_val.sum()\n",
    "freq_targets_val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_cli_validation['nas_steatosis_ordinal'].value_counts(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "targets_val = pd.DataFrame(targets_dict_val).rename_axis('Sample ID', axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pred_val = final_predictor.predict(indices=targets_val.index)\n",
    "score_val = final_predictor.predict_score(indices=targets_val.index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "source_data_fig5b=pd.concat([pred_val, score_val], axis=1).reset_index().drop('Sample ID', axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_cases = ['F2','F3', 'I2', 'S1']\n",
    "cm_scores = {}\n",
    "for test_case in test_cases:\n",
    "    overlap_ids = set(pred_val[test_case].dropna().index) & set(targets_val[test_case].dropna().index)\n",
    "    cm_scores[test_case]=confusion_matrix_scorer(targets_val.loc[overlap_ids][test_case], pred_val.loc[overlap_ids][test_case]) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "final_model_val = pd.DataFrame({i:performance_score_metric(cm_scores[i]) for i in ['F2', 'F3', 'I2', 'S1']})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "roc_auc_scores = {}\n",
    "roc_curves = {}\n",
    "for test_case in test_cases:\n",
    "    overlap_ids = set(score_val[test_case].dropna().index) & set(targets_val[test_case].dropna().index)\n",
    "    roc_auc_scores[test_case]=roc_auc_score(targets_val.loc[overlap_ids][test_case], score_val.loc[overlap_ids][test_case]).round(2)\n",
    "    fpr, tpr, thresholds = roc_curve(targets_val.loc[overlap_ids][test_case], score_val.loc[overlap_ids][test_case])\n",
    "    roc_curves[test_case]= (fpr, tpr, thresholds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "final_model_val.loc['ROC-AUC'] = list(roc_auc_scores.values())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "final_model_val.T.to_csv('tables/final_model_validation_cohort.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(4,4))\n",
    "lw = 2\n",
    "colors=['navy', 'darkred', 'gray']\n",
    "for test_case, color in zip(roc_curves.keys(), colors):\n",
    "    fpr = roc_curves[test_case][0]\n",
    "    tpr = roc_curves[test_case][1]\n",
    "    roc_auc = roc_auc_scores[test_case]\n",
    "    plt.plot(fpr, tpr, color=color, lw=lw, \n",
    "            label='ROC curve of {} (area = {})'.format(test_case, roc_auc_scores[test_case]))\n",
    "plt.plot([0, 1], [0, 1], color='gray', lw=lw, linestyle='--')\n",
    "plt.xlim([0.0, 1.0])\n",
    "plt.ylim([0.0, 1.05])\n",
    "plt.xlabel('False Positive Rate')\n",
    "plt.ylabel('True Positive Rate')\n",
    "plt.title('Validation cohort\\nReceiver operating characteristic', fontsize=12)\n",
    "plt.legend(loc=\"lower right\")\n",
    "plt.savefig('figures/ROC curve_validation.pdf', dpi=120, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Assess roc-auc of other tests"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = data_cli_validation[['Sample ID', 'te_1', 'swe', 'elf_b', 'cap_1', 'p3np_b', 'fib4', 'apri', 'forns', 'aar', 'alt', 'ast']]\n",
    "df1 = df.merge(targets_val.reset_index(), on='Sample ID').set_index('Sample ID')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "comparators = pd.read_csv('data/raw/comparators_validation_cohort.csv', index_col=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "roc_aucs_comparators = []\n",
    "for test, test_case in zip(comparators['comparators'], comparators['test_case']):\n",
    "    df = df1[[test,test_case]].dropna()\n",
    "    auroc = roc_auc_score(df[test_case], df[test]).round(2)\n",
    "    roc_aucs_comparators.append(auroc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "comparators['roc_auc'] = roc_aucs_comparators"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "proteomics_validation = pd.DataFrame({'comparators':roc_auc_scores.keys(), \n",
    "                                      'test_case':roc_auc_scores.keys(),\n",
    "                                      'roc_auc':roc_auc_scores.values()})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "func_sort = lambda x: x.sort_values(by='roc_auc', ascending = False)\n",
    "df_plot=comparators.append(proteomics_validation).groupby('test_case').apply(func_sort)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "source_data_fig5c=df_plot.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig=plt.subplots(figsize=(2.8,5))\n",
    "plt.barh(y=df_plot['comparators'], width=df_plot['roc_auc'], height=0.6)\n",
    "plt.gca().invert_yaxis()\n",
    "plt.xlim(0, 1)\n",
    "plt.savefig('figures/validation_roc_auc_comparison.pdf', dpi=120, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Write results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "excel_sheets = {'at_risk_results': at_risk_results_tab,\n",
    "                'at_risk_pred':at_risk_pred,\n",
    "                'at_risk_score': at_risk_score,\n",
    "                'healthy_results': healthy_results_tab, \n",
    "                'healthy_pred': healthy_pred,\n",
    "                'heatlhy_score': healthy_risk_score}\n",
    "FILE_FINAL_MODEL_EVALUATION = os.path.join(TABLEFOLDER ,'final_model_evaluation.xlsx')\n",
    "\n",
    "with pd.ExcelWriter(FILE_FINAL_MODEL_EVALUATION) as writer:\n",
    "    for sheet_name, _df in excel_sheets.items():\n",
    "        _df.to_excel(writer, sheet_name=sheet_name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SOURCE_DATA_FIG5 = os.path.join('data/processed/', 'source_data_fig5.xlsx')\n",
    "\n",
    "with pd.ExcelWriter(SOURCE_DATA_FIG5) as writer:\n",
    "    source_data_fig5b.to_excel(writer, sheet_name='fig5b')\n",
    "    source_data_fig5c.to_excel(writer, sheet_name='fig5c')\n",
    "    at_risk_results_tab.to_excel(writer, sheet_name='fib5a1')\n",
    "    healthy_results_tab.to_excel(writer, sheet_name='fib5a2')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### NRI calclulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_pred_val = pred_val.rename({i: i+'proteomics_pred' for i in pred_val.columns}, axis=1)\n",
    "cols_to_keep = [['te_1', 'swe', 'elf_b', 'cap_1', 'fib4', 'apri', 'aar','F2',\n",
    "                 'F2proteomics_pred', 'I2', 'I2proteomics_pred', 'S1proteomics_pred', 'S1']]\n",
    "df_nri_val = pd.concat([df1, df_pred_val], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_nri_val = df_nri_val.rename({'te_1':'te', 'elf_b':'elf', 'cap_1':'cap'}, axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_ald_pred = pred_final_test.copy()\n",
    "df_ald_pred = df_ald_pred.rename({i: i+'proteomics_pred' for i in df_ald_pred.columns}, axis=1)\n",
    "df_ald_pred2 = df_ald_pred.merge(targets_final_test, left_on=df_ald_pred.index, right_on=targets_final_test.index, left_index=True).drop('key_0', axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cols_to_keep = ['te', 'swe', 'elf', 'fib4', 'apri', 'aar', 'cap']\n",
    "df_nri_ald = pd.concat([df_ald_pred2, data_cli[cols_to_keep]], axis=1).loc[df_ald_pred2.index]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_nri_combined = pd.concat([df_nri_val, df_nri_ald])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# cutoffs in col_pairs >, cutoffs in col_pairs>=\n",
    "col_pairs = [('te_pred', 'te', 7.0), \n",
    "            ('swe_pred', 'swe', 8.6),\n",
    "            ('elf_pred', 'elf', 7.7),\n",
    "            ('fib4_pred', 'fib4', 1.45),\n",
    "            ('apri_pred', 'apri', 0.5),\n",
    "            ('aar_pred', 'aar', 2),\n",
    "            ('cap_pred', 'cap', 290),]\n",
    "col_pairs2 = [('te_pred_8kpa', 'te', 8.0),\n",
    "             ('fib4_pred_1.3', 'fib4', 1.3)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i, j, k in col_pairs:\n",
    "    sub_df = df_nri_combined[[j]].dropna() > k\n",
    "    new_col = sub_df[j].replace({True:1, False:0})\n",
    "    df_nri_combined[i] = new_col\n",
    "for i,j,k in col_pairs2:\n",
    "    sub_df = df_nri_combined[[j]].dropna()>= k\n",
    "    new_col = sub_df[j].replace({True:1, False:0})\n",
    "    df_nri_combined[i] = new_col"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_nri(data_frame, col_oldmodel, col_newmodel, col_label):\n",
    "    df=data_frame.copy()\n",
    "    df=df[[col_oldmodel, col_newmodel, col_label]].dropna()\n",
    "    event_count=df[df[col_label]==1].shape[0]\n",
    "    nonevent_count=df[df[col_label]==0].shape[0]\n",
    "    up_event=df[(df[col_label]==1) & (df[col_oldmodel]==0) & (df[col_newmodel]==1)].shape[0]\n",
    "    down_event=df[(df[col_label]==1) & (df[col_oldmodel]==1) & (df[col_newmodel]==0)].shape[0]\n",
    "    down_nonevent=df[(df[col_label]==0) & (df[col_oldmodel]==1) & (df[col_newmodel]==0)].shape[0]\n",
    "    up_nonevent=df[(df[col_label]==0) & (df[col_oldmodel]==0) & (df[col_newmodel]==1)].shape[0]\n",
    "    nri_e = (up_event-down_event)/event_count\n",
    "    nri_ne = (down_nonevent-up_nonevent)/nonevent_count\n",
    "    nri = nri_e + nri_ne\n",
    "    nr_obs = df.shape[0]\n",
    "    results = [up_event, down_event, down_nonevent, up_nonevent, event_count, nonevent_count, nr_obs, nri_e, nri_ne, nri]\n",
    "    index = ['up_event', 'down_event', 'down_nonevent', 'up_nonevent', 'event', 'nonevent', 'nr_observations', 'nri_e', 'nri_ne', 'nri']\n",
    "    results=pd.DataFrame(results, index=index, columns=[col_oldmodel])\n",
    "    return (results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_nri_f2 = []\n",
    "for baseline_model in ['te_pred','te_pred_8kpa', 'swe_pred', 'elf_pred', 'fib4_pred', 'fib4_pred_1.3','apri_pred']:\n",
    "    result= calculate_nri(data_frame=df_nri_combined, col_oldmodel=baseline_model, col_newmodel='F2proteomics_pred', col_label='F2')\n",
    "    results_nri_f2.append(result)\n",
    "results_nri_f2 = pd.concat(results_nri_f2, axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_nri_i2 = calculate_nri(data_frame=df_nri_combined, col_oldmodel='aar_pred', col_newmodel='I2proteomics_pred', col_label='I2')\n",
    "results_nri_s1 = calculate_nri(data_frame=df_nri_combined, col_oldmodel='cap_pred', col_newmodel='S1proteomics_pred', col_label='S1')\n",
    "results_nri = pd.concat([results_nri_f2, results_nri_i2, results_nri_s1], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_nri.to_csv('tables/nri_values.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_nri"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extract predicted positive cases and concordance with other non-invasive markers\n",
    "- F2: swe , te. Cut-offs for ≥F2: 8.6kPa for swe and 7kPa for te\n",
    "- I2: no equivalent \n",
    "- S1: cap. Cut-offs for ≥S1: 290\n",
    "- Values above cut-offs are color-coded"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "    doubleID=pd.read_csv('data/raw/DoubleIDkey.csv', index_col=False)\n",
    "    display(doubleID)\n",
    "except:\n",
    "    doubleID=None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### At risk population"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "at_risk_pred_pos_all = at_risk_pred[at_risk_pred.max(axis=1)==1]\n",
    "data_cli_atrisk = data_cli_full[data_cli_full['kleiner']==0.5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "at_risk_pos = at_risk_pred_pos_all.join(data_cli_atrisk[['swe', 'te', 'cap']], how='left')\n",
    "if not doubleID.empty:\n",
    "    at_risk_pos = at_risk_pos.join(doubleID.set_index('Sample ID'), how='left')\n",
    "m=at_risk_pos.style.apply(lambda x: ['background: yellow' if v>0 else \"\" for v in x],\n",
    "                        subset=['F2', 'I2', 'S1'], axis=1)\n",
    "m.apply(lambda x:['background: pink' if v>290 else\"\" for v in x],\n",
    "             subset=['cap'], axis=0)\n",
    "os.makedirs(os.path.join('tables', 'validation'), exist_ok=True)\n",
    "m.data.to_csv('tables/validation/predicted_scores_at_risk_postive.csv')\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Healthy population"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "healthy_pred_pos_all = healthy_pred[healthy_pred.max(axis=1)==1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "healthy_pos = healthy_pred_pos_all.join(data_cli[['swe', 'te', 'cap']], how='left')\n",
    "if not doubleID.empty:\n",
    "    healthy_pos = healthy_pos.join(doubleID.set_index('Sample ID'), how='left')\n",
    "l=healthy_pos.style.apply(lambda x: ['background: yellow' if v>0 else \"\" for v in x],\n",
    "                        subset=['F2', 'I2', 'S1'], axis=1)\n",
    "l.apply(lambda x:['background: pink' if v>290 else\"\" for v in x],\n",
    "             subset=['cap'], axis=0)\n",
    "l.data.to_csv('tables/validation/predicted_scores_healthy_postive.csv')\n",
    "l"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot predicted score distribution from the final model in three populations\n",
    "- healthy population (N=136).\n",
    "- at risk population (N=98). Heavy alcohol drinkers but benigh liver\n",
    "- high risk population (N=72). Healvy alcohol drinkers with various forms of liver injuries "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "final_model_test = final_predictor.predict(indices=train_test_split_final[0][1])\n",
    "final_model_test_score = final_predictor.predict_score(indices=train_test_split_final[0][1])\n",
    "final_model_test_score = final_model_test_score.sort_values(by=list(final_model_test_score.columns), ascending=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_naming = {'F2':'Fibrosis F2-F4', 'I2':'NAS Inflamation $\\geq 2$', 'S1':'NAS Steatosis $\\geq 5$%'}\n",
    "def plot_score_distribution(ax, endpoint):\n",
    "    ax = sns.kdeplot(data=healthy_risk_score.data[endpoint], label='healthy', color='steelblue', shade=True, ax=ax, cut=0)\n",
    "    ax = sns.kdeplot(data=at_risk_score.data[endpoint], label='at risk', color='darkorange', shade=True, ax=ax, cut=0)\n",
    "    ax = sns.kdeplot(data=final_model_test_score[endpoint], label='high risk', color='red', shade=True, ax=ax, cut=0)\n",
    "    ax.set_xlabel('Risk score predicted by\\n proteomic model')\n",
    "    ax.set_xlim(-0.0, 1.0)\n",
    "    ax.set_title(model_naming[endpoint])\n",
    "    return ax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(ncols=3,figsize=(12,4))\n",
    "for i, endpoint in enumerate(['F2', 'I2', 'S1']):\n",
    "    plot_score_distribution(axes[i], endpoint)\n",
    "fig.savefig(FIGURE_FOLDER / 'Score_kde_all', dpi=120, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for endpoint in ['F2', 'I2', 'S1']:\n",
    "    fig, ax = plt.subplots(figsize=(4,4))\n",
    "    plot_score_distribution(ax, endpoint)\n",
    "    fig.savefig(FIGURE_FOLDER / 'Score_kde_{}'.format(endpoint), dpi=120, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_score_hist(ax, data, label, color, bins=10, kde=False, rug=True):\n",
    "    ax = sns.distplot(data, ax=ax, label=label, color=color, bins=bins, kde=kde, rug=rug)\n",
    "    ax.set_xlabel('Risk score predicted by\\n proteomic model')\n",
    "    ax.set_xlim(-0.0, 1.0)\n",
    "    ax.set_title(label)\n",
    "    return \n",
    "\n",
    "for endpoint in ['F2', 'I2', 'S1']:\n",
    "    fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(12,4))\n",
    "    plot_score_hist(ax=axes[0], data=healthy_risk_score.data[endpoint], bins=10, label='healthy',   color='steelblue', kde=False, rug=True)\n",
    "    plot_score_hist(ax=axes[1], data=at_risk_score.data[endpoint],      bins=10, label='at risk',   color='darkorange',kde=False, rug=True)\n",
    "    plot_score_hist(ax=axes[2], data=final_model_test_score[endpoint],  bins=10, label='high risk', color='red',       kde=False, rug=True)\n",
    "    fig.suptitle(model_naming[endpoint])\n",
    "    fig.savefig(FIGURE_FOLDER / 'Score_hist_{}'.format(endpoint), dpi=120, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The rugs (small vertical lines at the botoom of the figure) indicate a single prediction. The scores are devided in 10 bins which means that each bin is decimal range (0.0 to 0.1, 0.1 to 0.2 and so forth)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(nrows=3, ncols=3, figsize=(12,12))\n",
    "for i, endpoint in enumerate(['F2', 'I2', 'S1']):\n",
    "    sns.distplot(healthy_risk_score.data[endpoint], ax=axes[i, 0], bins=10, label='healthy',   axlabel=False, color='steelblue', kde=False, rug=True)\n",
    "    sns.distplot(at_risk_score.data[endpoint],      ax=axes[i, 1], bins=10, label='at risk',   axlabel=False, color='darkorange',kde=False, rug=True)\n",
    "    sns.distplot(final_model_test_score[endpoint],  ax=axes[i, 2], bins=10, label='high risk', axlabel=False, color='red',       kde=False, rug=True)\n",
    "\n",
    "pad = 5 # in point\n",
    "for i, (_title, _endpoint) in enumerate(zip(['healthy', 'at risk', 'high risk'], ['F2', 'I2', 'S1'])):\n",
    "    axes[-1, i].set_xlabel('Risk score predicted by\\n proteomic model')\n",
    "    axes[0, i].set_title(_title)\n",
    "    _ax = axes[i, 0]\n",
    "    _ax.set_ylabel('frequency')\n",
    "    _ax.annotate(_endpoint, xy=(0, 0.5), \n",
    "                 xytext=(-_ax.yaxis.labelpad - pad, 0),\n",
    "                xycoords=_ax.yaxis.label, textcoords='offset points',\n",
    "                size='large', ha='right', va='center')\n",
    "        \n",
    "_ = fig.suptitle('Histograms by endpoint', y=.93, fontsize=16 )\n",
    "fig.savefig(FIGURE_FOLDER / 'Score_hist_all', dpi=120, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Prognostication"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_prog=pd.read_csv('data/raw/Prognostication_performance_20211018.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "title_name = {'LRE': 'Liver-related events', 'ACM': 'All-cause mortality'}\n",
    "for i in ['LRE', 'ACM']:\n",
    "    fig = plt.subplots(figsize=(2.5,4.5))\n",
    "    results_prog = results_prog.sort_values(by=i+ \" Harrell's C\", ascending =False)\n",
    "    bar_list = plt.barh(width=results_prog[i + \" Harrell's C\"], y=results_prog['Risk model'], \n",
    "             height=0.6, yerr=results_prog[i + ' SD'], facecolor='white', edgecolor='navy')\n",
    "    plt.gca().invert_yaxis()\n",
    "    plt.xlabel(\"Harrell's C\", fontsize=12)\n",
    "    plt.yticks(fontsize=12);\n",
    "    plt.xticks(fontsize=12);\n",
    "    plt.title(title_name[i], fontsize=14);\n",
    "    bar_list[0].set_color('steelblue')\n",
    "    plt.savefig(os.path.join(FIGURE_FOLDER, 'prognostication_{}.pdf'.format(i)), dpi=120, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "title_name = {'3 year LRE': '3-year liver-related events',\n",
    "              '5 year LRE': '5-year liver-related events',\n",
    "              '3 year ACM': '3-year all-cause mortality',\n",
    "             '5 year ACM': '5-year all-cause mortality'}\n",
    "for i in ['3 year LRE', '5 year LRE', '3 year ACM', '5 year ACM']:\n",
    "    fig = plt.subplots(figsize=(2.5,4.5))\n",
    "    results_prog = results_prog.sort_values(by=i+ \" AUC\", ascending =False)\n",
    "    bar_list = plt.barh(width=results_prog[i + \" AUC\"], y=results_prog['Risk model'], \n",
    "             height=0.6, yerr=results_prog[i + ' SD'], facecolor='white', edgecolor='navy')\n",
    "    plt.gca().invert_yaxis()\n",
    "    plt.xlabel(\"ROC-AUC\", fontsize=12)\n",
    "    plt.yticks(fontsize=12);\n",
    "    plt.xticks(fontsize=12);\n",
    "    plt.title(title_name[i], fontsize=14);\n",
    "    bar_list[0].set_color('steelblue')\n",
    "    plt.savefig(os.path.join(FIGURE_FOLDER, 'prognostication_{}.pdf'.format(i)), dpi=120, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "# Python Package Versions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "!pip list | grep pandas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "!pip list | grep scikit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Find the packages in the `requirements.txt` or `environment.yml`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %load environment.yml"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.10"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  },
  "toc-autonumbering": true,
  "toc-showcode": false,
  "toc-showmarkdowntxt": false,
  "toc-showtags": true
 },
 "nbformat": 4,
 "nbformat_minor": 4
}