Code14_Obesity-classifier_Arivale-gut-microbiome.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multiomics BMI Paper — Gut Microbiome-based Obesity Classifier in the Arivale Cohort"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***Analyzed by Tomasz Wilmanski originally, and modified by Kengo Watanabe***  \n",
    "\n",
    "This Jupyter Notebook (with Python 3 kernel) generated the random forest models for classifying participants into normal vs. obese class (based on either BMI or MetBMI) from the Arivale baseline gut microbiome dataset, and calculated the testing (hold-out) set-derived class (label and probability) predictions for the Arivale cohort. Because DeLong's test did not seem available in Python yet, DeLong's test was performed in another sub-notebook with R kernel.  \n",
    "\n",
    "Input files:  \n",
    "* Arivale baseline gut microbiome taxon abundances: 220902_Multiomics-BMI-NatMed1stRevision_Microbiome-DataCleaning_AlphaDiversity-and-TaxonAbundance_final.tsv  \n",
    "* Arivale baseline BMI and MetBMI: 220803_Multiomics-BMI-NatMed1stRevision_DeltaBMI-misclassification_biologicalBMI-baseline-summary-BothSex.tsv  \n",
    "* ROC curve of the classifiers: 221010_Multiomics-BMI-NatMed1stRevision_Microbiome-RFclassifier-DeLong-ver5_Arivale-wenceslaus_\\[BMI/MetBMI\\]-ROC-curve.tsv  \n",
    "* DeLong's test result: 221010_Multiomics-BMI-NatMed1stRevision_Microbiome-RFclassifier-DeLong-ver5_Arivale-wenceslaus_result-summary.tsv  \n",
    "\n",
    "Output figures and tables:  \n",
    "* Figure 4c, 4d  \n",
    "* Tables for Supplementary Data 1, 10  \n",
    "* Intermediate tables for the sub-notebook (class predictions)  \n",
    "\n",
    "Original notebook (memo for my future tracing):  \n",
    "* wenceslaus:\\[JupyterLab HOME\\]/220621_Multiomics-BMI-NatMedRevision/221010_Multiomics-BMI-NatMed1stRevision_Microbiome-RFclassifier-ver5_Arivale-wenceslaus.ipynb  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# packages in environment at /opt/conda/envs/arivale-py3:\n",
      "#\n",
      "# Name                    Version                   Build  Channel\n",
      "_libgcc_mutex             0.1                 conda_forge    conda-forge\n",
      "_openmp_mutex             4.5                       1_gnu    conda-forge\n",
      "analytics                 0.1                      pypi_0    pypi\n",
      "argon2-cffi               21.1.0           py39h3811e60_0    conda-forge\n",
      "arivale-data-interface    0.1.0                    pypi_0    pypi\n",
      "async_generator           1.10                       py_0    conda-forge\n",
      "atk-1.0                   2.36.0               h3371d22_4    conda-forge\n",
      "attrs                     21.2.0             pyhd8ed1ab_0    conda-forge\n",
      "backcall                  0.2.0              pyh9f0ad1d_0    conda-forge\n",
      "backports                 1.0                        py_2    conda-forge\n",
      "backports.functools_lru_cache 1.6.4              pyhd8ed1ab_0    conda-forge\n",
      "biopython                 1.79             py39h3811e60_0    conda-forge\n",
      "bleach                    4.1.0              pyhd8ed1ab_0    conda-forge\n",
      "bokeh                     2.4.1            py39hf3d152e_1    conda-forge\n",
      "boto                      2.49.0                     py_0    conda-forge\n",
      "boto3                     1.19.2                   pypi_0    pypi\n",
      "botocore                  1.22.2                   pypi_0    pypi\n",
      "brotlipy                  0.7.0           py39h3811e60_1001    conda-forge\n",
      "bzip2                     1.0.8                h7f98852_4    conda-forge\n",
      "c-ares                    1.17.2               h7f98852_0    conda-forge\n",
      "ca-certificates           2022.6.15.1          ha878542_0    conda-forge\n",
      "cachecontrol              0.12.6                     py_0    conda-forge\n",
      "cairo                     1.16.0            h6cf1ce9_1008    conda-forge\n",
      "certifi                   2022.6.15.1        pyhd8ed1ab_0    conda-forge\n",
      "cffi                      1.14.6           py39h4bc2ebd_1    conda-forge\n",
      "chardet                   4.0.0            py39hf3d152e_1    conda-forge\n",
      "charset-normalizer        2.0.0              pyhd8ed1ab_0    conda-forge\n",
      "click                     8.0.3            py39hf3d152e_0    conda-forge\n",
      "colorama                  0.4.4              pyh9f0ad1d_0    conda-forge\n",
      "cryptography              35.0.0           py39h95dcef6_1    conda-forge\n",
      "cycler                    0.10.0                     py_2    conda-forge\n",
      "cython                    0.29.24          py39he80948d_0    conda-forge\n",
      "debugpy                   1.4.1            py39he80948d_0    conda-forge\n",
      "decorator                 4.4.2                      py_0    conda-forge\n",
      "defusedxml                0.7.1              pyhd8ed1ab_0    conda-forge\n",
      "english                   2020.7.0           pyhd8ed1ab_0    conda-forge\n",
      "entrypoints               0.3             py39hde42818_1002    conda-forge\n",
      "et_xmlfile                1.0.1                   py_1001    conda-forge\n",
      "expat                     2.4.1                h9c3ff4c_0    conda-forge\n",
      "font-ttf-dejavu-sans-mono 2.37                 hab24e00_0    conda-forge\n",
      "font-ttf-inconsolata      3.000                h77eed37_0    conda-forge\n",
      "font-ttf-source-code-pro  2.038                h77eed37_0    conda-forge\n",
      "font-ttf-ubuntu           0.83                 hab24e00_0    conda-forge\n",
      "fontconfig                2.13.1            hba837de_1005    conda-forge\n",
      "fonts-conda-ecosystem     1                             0    conda-forge\n",
      "fonts-conda-forge         1                             0    conda-forge\n",
      "freetype                  2.10.4               h0708190_1    conda-forge\n",
      "fribidi                   1.0.10               h36c2ea0_0    conda-forge\n",
      "future                    0.18.2           py39hf3d152e_3    conda-forge\n",
      "gdk-pixbuf                2.42.6               h04a7f16_0    conda-forge\n",
      "gettext                   0.19.8.1          h73d1719_1008    conda-forge\n",
      "giflib                    5.2.1                h36c2ea0_2    conda-forge\n",
      "graphite2                 1.3.13            h58526e2_1001    conda-forge\n",
      "graphviz                  2.49.1               h85b4f2f_0    conda-forge\n",
      "gtk2                      2.24.33              h539f30e_1    conda-forge\n",
      "gts                       0.7.6                h64030ff_2    conda-forge\n",
      "harfbuzz                  3.0.0                h83ec7ef_1    conda-forge\n",
      "hdmedians                 0.14.2           py39hce5d2b2_0    conda-forge\n",
      "icu                       68.1                 h58526e2_0    conda-forge\n",
      "idna                      3.1                pyhd3deb0d_0    conda-forge\n",
      "importlib-metadata        4.8.1            py39hf3d152e_0    conda-forge\n",
      "iniconfig                 1.1.1              pyh9f0ad1d_0    conda-forge\n",
      "interchange               2021.0.4           pyhd8ed1ab_0    conda-forge\n",
      "ipykernel                 6.4.2            py39hef51801_0    conda-forge\n",
      "ipython                   7.28.0           py39hef51801_0    conda-forge\n",
      "ipython_genutils          0.2.0                      py_1    conda-forge\n",
      "ipywidgets                7.6.5              pyhd8ed1ab_0    conda-forge\n",
      "jbig                      2.1               h7f98852_2003    conda-forge\n",
      "jedi                      0.18.0           py39hf3d152e_2    conda-forge\n",
      "jinja2                    3.0.2              pyhd8ed1ab_0    conda-forge\n",
      "jmespath                  0.10.0                   pypi_0    pypi\n",
      "joblib                    1.1.0              pyhd8ed1ab_0    conda-forge\n",
      "jpeg                      9d                   h36c2ea0_0    conda-forge\n",
      "jsonschema                4.1.2              pyhd8ed1ab_0    conda-forge\n",
      "jupyter_client            7.0.6              pyhd8ed1ab_0    conda-forge\n",
      "jupyter_core              4.8.1            py39hf3d152e_0    conda-forge\n",
      "jupyterlab_pygments       0.1.2              pyh9f0ad1d_0    conda-forge\n",
      "jupyterlab_widgets        1.0.2              pyhd8ed1ab_0    conda-forge\n",
      "kiwisolver                1.3.2            py39h1a9c180_0    conda-forge\n",
      "krb5                      1.19.2               hcc1bbae_2    conda-forge\n",
      "lcms2                     2.12                 hddcbb42_0    conda-forge\n",
      "ld_impl_linux-64          2.36.1               hea4e1c9_2    conda-forge\n",
      "lerc                      2.2.1                h9c3ff4c_0    conda-forge\n",
      "libblas                   3.9.0           12_linux64_openblas    conda-forge\n",
      "libcblas                  3.9.0           12_linux64_openblas    conda-forge\n",
      "libcurl                   7.79.1               h2574ce0_1    conda-forge\n",
      "libdeflate                1.7                  h7f98852_5    conda-forge\n",
      "libedit                   3.1.20191231         he28a2e2_2    conda-forge\n",
      "libev                     4.33                 h516909a_1    conda-forge\n",
      "libffi                    3.4.2                h9c3ff4c_4    conda-forge\n",
      "libgcc-ng                 11.2.0              h1d223b6_11    conda-forge\n",
      "libgd                     2.3.3                h6ad9fb6_0    conda-forge\n",
      "libgfortran-ng            11.2.0              h69a702a_11    conda-forge\n",
      "libgfortran5              11.2.0              h5c6108e_11    conda-forge\n",
      "libglib                   2.70.0               h174f98d_1    conda-forge\n",
      "libgomp                   11.2.0              h1d223b6_11    conda-forge\n",
      "libiconv                  1.16                 h516909a_0    conda-forge\n",
      "liblapack                 3.9.0           12_linux64_openblas    conda-forge\n",
      "libnghttp2                1.43.0               h812cca2_1    conda-forge\n",
      "libopenblas               0.3.18          pthreads_h8fe5266_0    conda-forge\n",
      "libpng                    1.6.37               h21135ba_2    conda-forge\n",
      "libpq                     13.3                 hd57d9b9_1    conda-forge\n",
      "librsvg                   2.52.2               hc3c00ef_0    conda-forge\n",
      "libsodium                 1.0.18               h36c2ea0_1    conda-forge\n",
      "libssh2                   1.10.0               ha56f1ee_2    conda-forge\n",
      "libstdcxx-ng              11.2.0              he4da1e4_11    conda-forge\n",
      "libtiff                   4.3.0                hf544144_1    conda-forge\n",
      "libtool                   2.4.6             h9c3ff4c_1008    conda-forge\n",
      "libuuid                   2.32.1            h7f98852_1000    conda-forge\n",
      "libwebp                   1.2.1                h3452ae3_0    conda-forge\n",
      "libwebp-base              1.2.1                h7f98852_0    conda-forge\n",
      "libxcb                    1.13              h7f98852_1003    conda-forge\n",
      "libxml2                   2.9.12               h72842e0_0    conda-forge\n",
      "libxslt                   1.1.33               h15afd5d_2    conda-forge\n",
      "libzlib                   1.2.11            h36c2ea0_1013    conda-forge\n",
      "lockfile                  0.12.2                     py_1    conda-forge\n",
      "lxml                      4.6.3            py39h107f48f_0    conda-forge\n",
      "lz4-c                     1.9.3                h9c3ff4c_1    conda-forge\n",
      "markupsafe                2.0.1            py39h3811e60_0    conda-forge\n",
      "matplotlib-base           3.4.3            py39h2fa2bec_1    conda-forge\n",
      "matplotlib-inline         0.1.3              pyhd8ed1ab_0    conda-forge\n",
      "missingpy                 0.2.0                    pypi_0    pypi\n",
      "mistune                   0.8.4           py39h3811e60_1004    conda-forge\n",
      "more-itertools            8.10.0             pyhd8ed1ab_0    conda-forge\n",
      "mscorefonts               0.0.1                         3    conda-forge\n",
      "msgpack-python            1.0.2            py39h1a9c180_1    conda-forge\n",
      "natsort                   7.1.1              pyhd8ed1ab_0    conda-forge\n",
      "nbclient                  0.5.4              pyhd8ed1ab_0    conda-forge\n",
      "nbconvert                 6.2.0            py39hf3d152e_0    conda-forge\n",
      "nbformat                  5.1.3              pyhd8ed1ab_0    conda-forge\n",
      "ncurses                   6.2                  h58526e2_4    conda-forge\n",
      "neo4j-python-driver       4.3.7              pyhd8ed1ab_0    conda-forge\n",
      "neobolt                   1.7.17           py39h3811e60_2    conda-forge\n",
      "neotime                   1.7.4                      py_0    conda-forge\n",
      "nest-asyncio              1.5.1              pyhd8ed1ab_0    conda-forge\n",
      "networkx                  2.6.3              pyhd8ed1ab_0    conda-forge\n",
      "notebook                  6.4.5              pyha770c72_0    conda-forge\n",
      "numpy                     1.21.3           py39hdbf815f_0    conda-forge\n",
      "olefile                   0.46               pyh9f0ad1d_1    conda-forge\n",
      "openjpeg                  2.4.0                hb52868f_1    conda-forge\n",
      "openpyxl                  3.0.9              pyhd8ed1ab_0    conda-forge\n",
      "openssl                   1.1.1l               h7f98852_0    conda-forge\n",
      "packaging                 21.0               pyhd8ed1ab_0    conda-forge\n",
      "pandas                    1.3.4            py39hde0f152_0    conda-forge\n",
      "pandoc                    2.15                 h7f98852_0    conda-forge\n",
      "pandocfilters             1.5.0              pyhd8ed1ab_0    conda-forge\n",
      "pango                     1.48.10              h54213e6_2    conda-forge\n",
      "parso                     0.8.2              pyhd8ed1ab_0    conda-forge\n",
      "patsy                     0.5.2              pyhd8ed1ab_0    conda-forge\n",
      "pcre                      8.45                 h9c3ff4c_0    conda-forge\n",
      "pexpect                   4.8.0              pyh9f0ad1d_2    conda-forge\n",
      "pickleshare               0.7.5           py39hde42818_1002    conda-forge\n",
      "pillow                    8.3.2            py39ha612740_0    conda-forge\n",
      "pip                       21.3.1             pyhd8ed1ab_0    conda-forge\n",
      "pixman                    0.40.0               h36c2ea0_0    conda-forge\n",
      "plotly                    5.3.1              pyhd8ed1ab_0    conda-forge\n",
      "pluggy                    1.0.0            py39hf3d152e_1    conda-forge\n",
      "prometheus_client         0.11.0             pyhd8ed1ab_0    conda-forge\n",
      "prompt-toolkit            3.0.21             pyha770c72_0    conda-forge\n",
      "prompt_toolkit            3.0.21               hd8ed1ab_0    conda-forge\n",
      "psycopg2                  2.9.1            py39h3811e60_0    conda-forge\n",
      "pthread-stubs             0.4               h36c2ea0_1001    conda-forge\n",
      "ptyprocess                0.7.0              pyhd3deb0d_0    conda-forge\n",
      "py                        1.10.0             pyhd3deb0d_0    conda-forge\n",
      "py2neo                    2021.2.3           pyhd8ed1ab_0    conda-forge\n",
      "pycparser                 2.20               pyh9f0ad1d_2    conda-forge\n",
      "pygments                  2.10.0             pyhd8ed1ab_0    conda-forge\n",
      "pygraphviz                1.7              py39h78163bd_0    conda-forge\n",
      "pyopenssl                 21.0.0             pyhd8ed1ab_0    conda-forge\n",
      "pyparsing                 2.4.7              pyh9f0ad1d_0    conda-forge\n",
      "pyrsistent                0.17.3           py39h3811e60_2    conda-forge\n",
      "pysam                     0.17.0           py39h051187c_0    bioconda\n",
      "pysocks                   1.7.1            py39hf3d152e_3    conda-forge\n",
      "pytabix                   0.1              py39h98c8e45_1    bioconda\n",
      "pytest                    6.2.5            py39hf3d152e_0    conda-forge\n",
      "python                    3.9.7           hb7a2778_3_cpython    conda-forge\n",
      "python-dateutil           2.8.2              pyhd8ed1ab_0    conda-forge\n",
      "python-graphviz           0.17               pyhaef67bd_0    conda-forge\n",
      "python_abi                3.9                      2_cp39    conda-forge\n",
      "pytz                      2021.3             pyhd8ed1ab_0    conda-forge\n",
      "pyvcf                     0.6.8           py39hde42818_1002    conda-forge\n",
      "pyyaml                    6.0              py39h3811e60_0    conda-forge\n",
      "pyzmq                     22.3.0           py39h37b5a0c_0    conda-forge\n",
      "readline                  8.1                  h46c0cb4_0    conda-forge\n",
      "requests                  2.26.0             pyhd8ed1ab_0    conda-forge\n",
      "s3transfer                0.5.0                    pypi_0    pypi\n",
      "scikit-bio                0.5.6            py39h16ac069_4    conda-forge\n",
      "scikit-learn              1.0.1            py39h7c5d8c9_0    conda-forge\n",
      "scipy                     1.7.1            py39hee8e79c_0    conda-forge\n",
      "seaborn                   0.11.2               hd8ed1ab_0    conda-forge\n",
      "seaborn-base              0.11.2             pyhd8ed1ab_0    conda-forge\n",
      "send2trash                1.8.0              pyhd8ed1ab_0    conda-forge\n",
      "setuptools                58.2.0           py39hf3d152e_0    conda-forge\n",
      "six                       1.16.0             pyh6c4a22f_0    conda-forge\n",
      "sqlite                    3.36.0               h9cd32fc_2    conda-forge\n",
      "statsmodels               0.13.0           py39hce5d2b2_0    conda-forge\n",
      "tenacity                  8.0.1              pyhd8ed1ab_0    conda-forge\n",
      "terminado                 0.12.1           py39hf3d152e_0    conda-forge\n",
      "testpath                  0.5.0              pyhd8ed1ab_0    conda-forge\n",
      "threadpoolctl             3.0.0              pyh8a188c0_0    conda-forge\n",
      "tk                        8.6.11               h27826a3_1    conda-forge\n",
      "toml                      0.10.2             pyhd8ed1ab_0    conda-forge\n",
      "tornado                   6.1              py39h3811e60_1    conda-forge\n",
      "traitlets                 5.1.0              pyhd8ed1ab_0    conda-forge\n",
      "typing_extensions         3.10.0.2           pyha770c72_0    conda-forge\n",
      "tzdata                    2021e                he74cb21_0    conda-forge\n",
      "urllib3                   1.26.7             pyhd8ed1ab_0    conda-forge\n",
      "wcwidth                   0.2.5              pyh9f0ad1d_2    conda-forge\n",
      "webencodings              0.5.1                      py_1    conda-forge\n",
      "wheel                     0.37.0             pyhd8ed1ab_1    conda-forge\n",
      "widgetsnbextension        3.5.1            py39hf3d152e_4    conda-forge\n",
      "wordcloud                 1.8.1                    pypi_0    pypi\n",
      "xlrd                      2.0.1              pyhd8ed1ab_3    conda-forge\n",
      "xorg-kbproto              1.0.7             h7f98852_1002    conda-forge\n",
      "xorg-libice               1.0.10               h7f98852_0    conda-forge\n",
      "xorg-libsm                1.2.3             hd9c2040_1000    conda-forge\n",
      "xorg-libx11               1.7.2                h7f98852_0    conda-forge\n",
      "xorg-libxau               1.0.9                h7f98852_0    conda-forge\n",
      "xorg-libxdmcp             1.1.3                h7f98852_0    conda-forge\n",
      "xorg-libxext              1.3.4                h7f98852_1    conda-forge\n",
      "xorg-libxrender           0.9.10            h7f98852_1003    conda-forge\n",
      "xorg-renderproto          0.11.1            h7f98852_1002    conda-forge\n",
      "xorg-xextproto            7.3.0             h7f98852_1002    conda-forge\n",
      "xorg-xproto               7.0.31            h7f98852_1007    conda-forge\n",
      "xz                        5.2.5                h516909a_1    conda-forge\n",
      "yaml                      0.2.5                h516909a_0    conda-forge\n",
      "zeromq                    4.3.4                h9c3ff4c_1    conda-forge\n",
      "zipp                      3.6.0              pyhd8ed1ab_0    conda-forge\n",
      "zlib                      1.2.11            h36c2ea0_1013    conda-forge\n",
      "zstd                      1.5.0                ha95c52a_0    conda-forge\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scipy.stats as stats\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import seaborn as sns\n",
    "#For Arial font\n",
    "#!conda install -c conda-forge -y mscorefonts\n",
    "##-> The below was also needed in matplotlib 3.4.2\n",
    "#import shutil\n",
    "#import matplotlib\n",
    "#shutil.rmtree(matplotlib.get_cachedir())\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "from IPython.display import display\n",
    "import time\n",
    "\n",
    "import random\n",
    "import statsmodels.formula.api as smf\n",
    "from statsmodels.stats.anova import anova_lm\n",
    "from statsmodels.stats import multitest as multi\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.decomposition import PCA\n",
    "from decimal import Decimal, ROUND_HALF_UP\n",
    "from sklearn.model_selection import GridSearchCV\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.metrics import confusion_matrix\n",
    "from sklearn.metrics import roc_curve, precision_recall_curve, auc\n",
    "from sklearn.metrics import roc_auc_score, average_precision_score\n",
    "from statsmodels.stats import weightstats\n",
    "\n",
    "!conda list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Data preparation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> The necessary files were copied from the dalek server in advance.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1-1. Gut microbiome taxon abundance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Import cleaned table for baseline gut microbiome data\n",
    "fileDir = './ImportData/'\n",
    "ipynbName = '220902_Multiomics-BMI-NatMed1stRevision_Microbiome-DataCleaning_'\n",
    "fileName = 'AlphaDiversity-and-TaxonAbundance_final.tsv'\n",
    "tempDF = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t', dtype={'public_client_id': str})\n",
    "tempDF = tempDF.set_index('public_client_id')\n",
    "print('Before:', tempDF.shape)\n",
    "\n",
    "#Take only the taxon abundances\n",
    "tempL1 = ['Phylum', 'Class', 'Order', 'Family', 'Genus', 'Species']\n",
    "tempL2 = []\n",
    "for taxon_rank in tempL1:\n",
    "    tempL = tempDF.loc[:, tempDF.columns.str.contains(taxon_rank+':')].columns.tolist()\n",
    "    tempL = sorted(tempL)\n",
    "    for taxon in tempL:\n",
    "        tempL2.append(taxon)\n",
    "tempDF = tempDF[tempL2]\n",
    "print('After taking the target taxonomic ranks:', tempDF.shape)\n",
    "tempDF1 = tempDF.columns.to_series().str.split(pat=':', expand=True)\n",
    "display(tempDF1[0].value_counts())\n",
    "\n",
    "display(tempDF)\n",
    "\n",
    "biomeDF = tempDF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1-2. BMI and omics-inferred BMI classes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Import cleaned table for baseline measured and biological BMIs\n",
    "fileDir = './ImportData/'\n",
    "ipynbName = '220803_Multiomics-BMI-NatMed1stRevision_DeltaBMI-misclassification_'\n",
    "fileName = 'biologicalBMI-baseline-summary-BothSex.tsv'\n",
    "tempDF = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t', dtype={'public_client_id': str})\n",
    "tempDF = tempDF.set_index('public_client_id')\n",
    "print('Original:', len(tempDF))\n",
    "\n",
    "#Take the participants having gut microbiome data\n",
    "tempDF = tempDF.loc[tempDF.index.isin(biomeDF.index.tolist())]\n",
    "print(' -> with gut microbiome data:', len(tempDF))\n",
    "\n",
    "#Clean to handle easier in this notebook\n",
    "tempDF.columns = tempDF.columns.str.replace('Base', '')\n",
    "\n",
    "#Select the BMI class and covariates (just for the display in Jupyter notebook)\n",
    "tempL1 = tempDF.loc[:, tempDF.columns.str.contains('_class')].columns.tolist()\n",
    "tempL2 = ['BMI', 'Sex', 'Age', 'PC1', 'PC2', 'PC3', 'PC4', 'PC5']\n",
    "tempL = [col_n for sublist in [tempL1, tempL2] for col_n in sublist]\n",
    "tempDF = tempDF[tempL]\n",
    "\n",
    "display(tempDF)\n",
    "tempL = []\n",
    "for bmi_class in tempL1:\n",
    "    tempL.append(tempDF[bmi_class].value_counts())\n",
    "tempDF1 = pd.concat(tempL, axis=1)\n",
    "tempDF1 = tempDF1.sort_index(ascending=True)\n",
    "display(tempDF1)\n",
    "\n",
    "bmiDF = tempDF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1-3. Split the cohort into 5 sets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> In this case, the split sets are different between BMI and MetBMI classes. Hence, the randomization step is added to reduce the bias of common participant existence, which would happen if simply using the row order.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempD1 = {'BothSex':bmiDF}\n",
    "yvarL = ['BMI', 'MetBMI']#Classifier target in this study; Fix through this notebook\n",
    "classD = {'Normal':0, 'Obese':1}#Target class label and its code in this study; Fix through this notebook\n",
    "nmodels = 5#Fix through this notebook\n",
    "random.seed(123)#For reproducibility (while checking the following independency)\n",
    "\n",
    "tempD2 = {}\n",
    "for sex in tempD1.keys():\n",
    "    tempD3 = {}\n",
    "    for yvar in yvarL:\n",
    "        #Select target classes for the classifier\n",
    "        tempDF = tempD1[sex]\n",
    "        tempDF = tempDF.loc[tempDF[yvar+'_class'].isin(list(classD.keys()))]\n",
    "        #Randomize the row order\n",
    "        tempL = tempDF.index.tolist()\n",
    "        tempL = random.sample(tempL, len(tempL))\n",
    "        tempDF = tempDF.loc[tempL]\n",
    "        #Split cohort to define the training and testing (hold-out) sets\n",
    "        tempL = np.array_split(tempDF, nmodels)#List of DFs\n",
    "        tempD = {}\n",
    "        for model_k in range(nmodels):\n",
    "            tempDF1 = tempL[model_k]\n",
    "            model_n = 'Model_'+str(model_k+1).zfill(2)\n",
    "            tempS = pd.Series(np.repeat(model_n, len(tempDF1)),\n",
    "                              index=tempDF1.index, name='Testing_'+yvar)\n",
    "            tempD[model_k] = tempS\n",
    "        tempS = pd.concat(list(tempD.values()), axis=0)\n",
    "        tempD3[yvar] = tempS\n",
    "    tempDF = pd.concat(list(tempD3.values()), axis=1)#NaN for out-of-target participant\n",
    "    #Add the info to bmiDF\n",
    "    tempDF1 = tempD1[sex]\n",
    "    tempDF = pd.merge(tempDF1, tempDF, left_index=True, right_index=True, how='left')\n",
    "    tempD2[sex] = tempDF\n",
    "    \n",
    "    print(sex)\n",
    "    display(tempDF)\n",
    "    for yvar in yvarL:\n",
    "        tempDF1 = tempDF.loc[~tempDF['Testing_'+yvar].isnull()]\n",
    "        print(' - '+yvar)\n",
    "        display(tempDF1.describe(include='all'))\n",
    "        display(tempDF1['Testing_'+yvar].value_counts())\n",
    "    print('')\n",
    "#Update\n",
    "bmiDF = tempD2['BothSex']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1-4. Check independency among the set"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Perform Pearson's chi-squared test for categorical variables (using scipy library) and ANOVA for numeric variables (using statsmodels library). Note that scipy.stats.f_oneway() doesn't report degrees of freedom.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Compare the sets per classifier target (i.e., 5 sets per comparison)\n",
    "tempD1 = {'BothSex':bmiDF}\n",
    "for sex in tempD1.keys():\n",
    "    tempD2 = {}\n",
    "    for yvar in yvarL:\n",
    "        tempDF = tempD1[sex]\n",
    "        tempDF = tempDF.loc[~tempDF['Testing_'+yvar].isnull()]\n",
    "        tempL1 = []#For variable name\n",
    "        tempL2 = []#Test name\n",
    "        tempL3 = []#For degrees of freedom\n",
    "        tempL4 = []#For test statistic\n",
    "        tempL5 = []#For P-value\n",
    "        #Categorical variables\n",
    "        #tempL = tempDF.select_dtypes(exclude=[np.number]).columns.tolist()\n",
    "        tempL = [yvar+'_class', 'Sex']\n",
    "        for col_n in tempL:\n",
    "            if (col_n=='Sex')&(sex in ['Female', 'Male']):\n",
    "                continue\n",
    "            else:\n",
    "                tempDF1 = pd.crosstab(tempDF[col_n], tempDF['Testing_'+yvar])\n",
    "                #Pearson's chi-squared test\n",
    "                chi2, pval, dof, tempA = stats.chi2_contingency(tempDF1, correction=False)\n",
    "                tempL1.append(col_n)\n",
    "                tempL2.append('Pearson\\'s chi-squared test')\n",
    "                tempL3.append(dof)\n",
    "                tempL4.append(chi2)\n",
    "                tempL5.append(pval)\n",
    "        #Numeric variables\n",
    "        tempL = tempDF.select_dtypes(include=[np.number]).columns.tolist()\n",
    "        for col_n in tempL:\n",
    "            #ANOVA\n",
    "            formula = col_n+' ~ C(Testing_'+yvar+')'\n",
    "            model = smf.ols(formula, data=tempDF).fit()\n",
    "            tempDF1 = anova_lm(model, typ=1)#ANOVA type doesn't matter this case\n",
    "            tempL1.append(col_n)\n",
    "            tempL2.append('ANOVA')\n",
    "            dof1 = tempDF1['df'].astype('int64').loc['C(Testing_'+yvar+')']#Between-groups\n",
    "            dof2 = tempDF1['df'].astype('int64').loc['Residual']#Within-groups\n",
    "            tempL3.append((dof1, dof2))\n",
    "            tempL4.append(tempDF1['F'].loc['C(Testing_'+yvar+')'])\n",
    "            tempL5.append(tempDF1['PR(>F)'].loc['C(Testing_'+yvar+')'])\n",
    "        tempDF = pd.DataFrame({'Variable':tempL1,\n",
    "                               'StatisticalTest':tempL2,\n",
    "                               'N':len(tempDF),\n",
    "                               'DoF':tempL3,\n",
    "                               'Statistic':tempL4,\n",
    "                               'Pval':tempL5})\n",
    "        #P-value adjustment (within classifier target) by using Benjamini–Hochberg method\n",
    "        tempDF['AdjPval_within'] = multi.multipletests(tempDF['Pval'], alpha=0.05, method='fdr_bh',\n",
    "                                                       is_sorted=False, returnsorted=False)[1]\n",
    "        tempDF['Classifier'] = yvar+' class'\n",
    "        tempD2[yvar] = tempDF\n",
    "    tempDF = pd.concat(list(tempD2.values()), axis=0)\n",
    "    #P-value adjustment (across all tests) by using Benjamini–Hochberg method\n",
    "    tempDF['AdjPval_all'] = multi.multipletests(tempDF['Pval'], alpha=0.05, method='fdr_bh',\n",
    "                                                is_sorted=False, returnsorted=False)[1]\n",
    "    tempDF = tempDF.set_index(['Classifier', 'Variable'])\n",
    "    print(sex)\n",
    "    display(tempDF)\n",
    "    #Save\n",
    "    fileDir = './ExportData/'\n",
    "    ipynbName = '221010_Multiomics-BMI-NatMed1stRevision_Microbiome-RFclassifier-ver5_Arivale-wenceslaus_'\n",
    "    fileName = 'split-sets-independency_within-classifier-'+sex+'.tsv'\n",
    "    tempDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Compare the sets across classifier targets per BMI class (i.e., 10 sets per comparison)\n",
    "tempD1 = {'BothSex':bmiDF}\n",
    "for sex in tempD1.keys():\n",
    "    tempD2 = {}\n",
    "    for bmi_class in classD.keys():\n",
    "        tempD = {}\n",
    "        for yvar in yvarL:\n",
    "            tempDF = tempD1[sex]\n",
    "            tempDF = tempDF.loc[tempDF[yvar+'_class']==bmi_class]\n",
    "            tempDF['Testing'] = yvar+'_'+tempDF['Testing_'+yvar]\n",
    "            tempDF = tempDF.loc[:, ~tempDF.columns.str.contains('Testing_')]\n",
    "            tempD[yvar] = tempDF\n",
    "        tempDF = pd.concat(list(tempD.values()), axis=0)\n",
    "        tempL1 = []#For variable name\n",
    "        tempL2 = []#Test name\n",
    "        tempL3 = []#For degrees of freedom\n",
    "        tempL4 = []#For test statistic\n",
    "        tempL5 = []#For P-value\n",
    "        #Categorical variables\n",
    "        #tempL = tempDF.select_dtypes(exclude=[np.number]).columns.tolist()\n",
    "        tempL = ['Sex']\n",
    "        for col_n in tempL:\n",
    "            if (col_n=='Sex')&(sex in ['Female', 'Male']):\n",
    "                continue\n",
    "            else:\n",
    "                tempDF1 = pd.crosstab(tempDF[col_n], tempDF['Testing'])\n",
    "                #Pearson's chi-squared test\n",
    "                chi2, pval, dof, tempA = stats.chi2_contingency(tempDF1, correction=False)\n",
    "                tempL1.append(col_n)\n",
    "                tempL2.append('Pearson\\'s chi-squared test')\n",
    "                tempL3.append(dof)\n",
    "                tempL4.append(chi2)\n",
    "                tempL5.append(pval)\n",
    "        #Numeric variables\n",
    "        tempL = tempDF.select_dtypes(include=[np.number]).columns.tolist()\n",
    "        for col_n in tempL:\n",
    "            #ANOVA\n",
    "            formula = col_n+' ~ C(Testing)'\n",
    "            model = smf.ols(formula, data=tempDF).fit()\n",
    "            tempDF1 = anova_lm(model, typ=1)#ANOVA type doesn't matter this case\n",
    "            tempL1.append(col_n)\n",
    "            tempL2.append('ANOVA')\n",
    "            dof1 = tempDF1['df'].astype('int64').loc['C(Testing)']#Between-groups\n",
    "            dof2 = tempDF1['df'].astype('int64').loc['Residual']#Within-groups\n",
    "            tempL3.append((dof1, dof2))\n",
    "            tempL4.append(tempDF1['F'].loc['C(Testing)'])\n",
    "            tempL5.append(tempDF1['PR(>F)'].loc['C(Testing)'])\n",
    "        tempDF = pd.DataFrame({'Variable':tempL1,\n",
    "                               'StatisticalTest':tempL2,\n",
    "                               'N':len(tempDF),\n",
    "                               'DoF':tempL3,\n",
    "                               'Statistic':tempL4,\n",
    "                               'Pval':tempL5})\n",
    "        #P-value adjustment (within BMI class) by using Benjamini–Hochberg method\n",
    "        tempDF['AdjPval_within'] = multi.multipletests(tempDF['Pval'], alpha=0.05, method='fdr_bh',\n",
    "                                                       is_sorted=False, returnsorted=False)[1]\n",
    "        tempDF['BMIorMetBMIclass'] = bmi_class\n",
    "        tempD2[bmi_class] = tempDF\n",
    "    tempDF = pd.concat(list(tempD2.values()), axis=0)\n",
    "    #P-value adjustment (across all tests) by using Benjamini–Hochberg method\n",
    "    tempDF['AdjPval_all'] = multi.multipletests(tempDF['Pval'], alpha=0.05, method='fdr_bh',\n",
    "                                                is_sorted=False, returnsorted=False)[1]\n",
    "    tempDF = tempDF.set_index(['BMIorMetBMIclass', 'Variable'])\n",
    "    print(sex)\n",
    "    display(tempDF)\n",
    "    #Save\n",
    "    fileDir = './ExportData/'\n",
    "    ipynbName = '221010_Multiomics-BMI-NatMed1stRevision_Microbiome-RFclassifier-ver5_Arivale-wenceslaus_'\n",
    "    fileName = 'split-sets-independency_across-classifiers-'+sex+'.tsv'\n",
    "    tempDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Preprocessing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2-1. Standardization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Although standardization is unnecessary for RF classifiers in general, certain standardization would be required to harmonize the Arivale and TwinsUK datasets.  \n",
    "> –> In this version, the harmonization is not needed because classifiers are generated per cohort. However, standardizaiton is performed to apply PCA. Of note, log-transformation is applied before standardization since it was not done in the data cleaning notebook.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF = biomeDF\n",
    "tempD1 = {'BothSex':bmiDF}\n",
    "\n",
    "#Standardization\n",
    "tempD2 = {}\n",
    "for sex in tempD1.keys():\n",
    "    tempDF1 = tempD1[sex]\n",
    "    tempDF1 = tempDF.loc[tempDF1.index.tolist()]\n",
    "    \n",
    "    #Check\n",
    "    print(sex+':', tempDF1.shape)\n",
    "    print(' - Negative values in DF:', (tempDF<0).to_numpy().sum(axis=None))\n",
    "    print(' - Before:')\n",
    "    tempS = pd.Series(stats.skew(tempDF1), index=tempDF1.columns, name='Skewness')\n",
    "    display(tempS.describe())\n",
    "    \n",
    "    #log-transformation\n",
    "    tempDF1 = np.log(tempDF1 + 1)\n",
    "    \n",
    "    #Z-score transformation\n",
    "    scaler = StandardScaler(copy=True, with_mean=True, with_std=True)\n",
    "    tempA = scaler.fit_transform(tempDF1)#Column direction\n",
    "    tempDF1 = pd.DataFrame(data=tempA, index=tempDF1.index, columns=tempDF1.columns)\n",
    "    \n",
    "    tempD2[sex] = tempDF1\n",
    "    \n",
    "    #Confirmation\n",
    "    print(' - After:')\n",
    "    tempS = pd.Series(stats.skew(tempDF1), index=tempDF1.columns, name='Skewness')\n",
    "    display(tempS.describe())\n",
    "    tempDF2 = tempDF1.describe()\n",
    "    tempDF2.loc['Skewness'] = tempS\n",
    "    display(tempDF2)\n",
    "    sns.set(style='ticks', font='Arial', context='notebook')\n",
    "    plt.figure(figsize=(4, 3))\n",
    "    for col_i in range(3):\n",
    "        sns.distplot(tempDF1.iloc[:, col_i], label=tempDF1.columns[col_i])\n",
    "    sns.despine()\n",
    "    plt.xlabel(r'$Z$'+'-score')\n",
    "    plt.ylabel('Density')\n",
    "    plt.legend(bbox_to_anchor=(1, 0.5), loc='center left', borderaxespad=1)\n",
    "    plt.show()\n",
    "\n",
    "#Update\n",
    "biomeDF = tempD2['BothSex']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2-2. PCA for feature selection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Because of many features vs. small sample size, the input taxon abundances are applied to PCA.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempD1 = {'BothSex':biomeDF}\n",
    "#pc_cutoff = 0.1#Cutoff of the PC's explained variance [%]\n",
    "pc_topX = 50#Set the maximum number of input features\n",
    "nplot = 5\n",
    "tempD2 = {'BothSex':bmiDF}\n",
    "tempD3 = {'Underweight':'blue', 'Normal':'green', 'Overweight':'orange', 'Obese':'red'}\n",
    "\n",
    "#PCA for feature selection\n",
    "tempD4 = {}\n",
    "for sex in tempD1.keys():\n",
    "    print(sex)\n",
    "    tempDF = tempD1[sex]#Already standardized\n",
    "    \n",
    "    #Perform PCA\n",
    "    nPCs = np.min(tempDF.shape)\n",
    "    model = PCA(n_components=nPCs, svd_solver='randomized', iterated_power='auto', random_state=123)\n",
    "    model.fit(tempDF)\n",
    "    \n",
    "    #Explained variance\n",
    "    tempS = pd.Series(data=model.explained_variance_ratio_*100,\n",
    "                      index=['PC'+str(i+1) for i in range(nPCs)], name='ExplainedVariance')\n",
    "    print(' - Percentage of variance explained by each component:')\n",
    "    display(tempS)\n",
    "    ##Retrieve cutoff value at top X\n",
    "    pc_cutoff = tempS.iloc[(pc_topX-1)]\n",
    "    ##Scree plot\n",
    "    tempDF1 = tempS.reset_index()\n",
    "    tempDF1['PC'] = [i+1 for i in range(nPCs)]\n",
    "    sns.set(style='ticks', font='Arial', context='talk')\n",
    "    plt.figure(figsize=(4, 3))\n",
    "    p = sns.lineplot(data=tempDF1, x='PC', y='ExplainedVariance', color='k')\n",
    "    sns.despine()\n",
    "    p.set(xlim=(0.5, nPCs+0.5))\n",
    "    p.axhline(y=pc_cutoff, linestyle=\"--\", color='crimson', zorder=0)\n",
    "    plt.ylabel('Explained varaince [%]')\n",
    "    plt.xlabel('Principal component number')\n",
    "    plt.show()\n",
    "    ##Scree plot (log-scale)\n",
    "    tempDF1['ExplainedVariance_log10'] = np.log10(tempDF1['ExplainedVariance'])\n",
    "    sns.set(style='ticks', font='Arial', context='talk')\n",
    "    plt.figure(figsize=(4, 3))\n",
    "    p = sns.lineplot(data=tempDF1, x='PC', y='ExplainedVariance_log10', color='k')\n",
    "    sns.despine()\n",
    "    p.set(xlim=(0.5, nPCs+0.5))\n",
    "    p.axhline(y=np.log10(pc_cutoff), linestyle=\"--\", color='crimson', zorder=0)\n",
    "    plt.ylabel('Explained varaince [%]\\n(log-scale)')\n",
    "    plt.xlabel('Principal component number')\n",
    "    plt.show()\n",
    "    ##Integrate the explained variance into PC label\n",
    "    tempL = []\n",
    "    for i in range(nPCs):\n",
    "        round_value = Decimal(str(tempS['PC'+str(i+1)])).quantize(Decimal('0.1'), rounding=ROUND_HALF_UP)\n",
    "        tempL.append('PC'+str(i+1)+' ('+str(round_value)+'%)')\n",
    "    \n",
    "    #Select PCs as input features\n",
    "    tempS = tempS.loc[tempS>=pc_cutoff]\n",
    "    print(' - PCs for input features:')\n",
    "    print('   - Total explained variance:', tempS.sum())\n",
    "    display(tempS.describe())\n",
    "    \n",
    "    #Projection\n",
    "    tempDF1 = pd.DataFrame(data=model.transform(tempDF), index=tempDF.index, columns=tempL)\n",
    "    tempDF1 = tempDF1.iloc[:, :len(tempS)]\n",
    "    print(' - Projection DF:', tempDF1.shape)\n",
    "    display(tempDF1)\n",
    "    \n",
    "    #PC component\n",
    "    tempDF2 = pd.DataFrame(data=model.components_, index=tempL, columns=tempDF.columns)\n",
    "    tempDF2 = tempDF2.iloc[:len(tempS), :]\n",
    "    print('PC component DF:', tempDF2.shape)\n",
    "    display(tempDF2)\n",
    "    \n",
    "    #Save\n",
    "    fileDir = './ExportData/'\n",
    "    ipynbName = '221010_Multiomics-BMI-NatMed1stRevision_Microbiome-RFclassifier-ver5_Arivale-wenceslaus_'\n",
    "    fileName = 'TaxonAbundance-PCA-projection-'+sex+'.tsv'\n",
    "    tempDF1.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "    fileName = 'TaxonAbundance-PCA-component-'+sex+'.tsv'\n",
    "    tempDF2.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "    \n",
    "    tempD4[sex] = tempDF1\n",
    "    \n",
    "    #Visualize sample distribution in the projected spaces\n",
    "    tempDF = tempD2[sex]\n",
    "    tempS = tempDF['BMI_class']\n",
    "    tempDF = pd.merge(tempDF1.iloc[:, :(nplot+1)], tempS, left_index=True, right_index=True, how='left')\n",
    "    sns.set(style='ticks', font='Arial', context='talk')\n",
    "    p = sns.PairGrid(data=tempDF, hue='BMI_class', hue_order=tempD3.keys(), palette=tempD3)\n",
    "    p.map_lower(sns.scatterplot, edgecolor='0.3', alpha=0.5, s=15)\n",
    "    for i, j in zip(*np.triu_indices_from(p.axes, 0)):\n",
    "        p.axes[i, j].set_visible(False)\n",
    "    p.add_legend(bbox_to_anchor=(0.6, 0.6), loc='center right', frameon=True)\n",
    "    plt.show()\n",
    "    \n",
    "    print('')\n",
    "\n",
    "#Update\n",
    "biomeDF = tempD4['BothSex']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2-3. Class label encoding"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> In this case, only two labels: Normal vs. Obese. Hence, manually encode.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempD1 = {'BothSex':bmiDF}\n",
    "\n",
    "tempD2 = {}\n",
    "for sex in tempD1.keys():\n",
    "    tempDF = tempD1[sex]\n",
    "    #Encoding\n",
    "    for yvar in yvarL:\n",
    "        tempDF[yvar+'_class_code'] = tempDF[yvar+'_class'].map(classD)\n",
    "    tempD2[sex] = tempDF\n",
    "    \n",
    "    #Confirmation\n",
    "    print(sex+':', tempDF.shape)\n",
    "    display(tempDF.describe(include='all'))\n",
    "    print('')\n",
    "\n",
    "#Update\n",
    "bmiDF = tempD2['BothSex']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. BMI classifier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "yvar = 'BMI'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3-1. Select the target class participants"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "tempD1 = {'BothSex':bmiDF}\n",
    "tempD2 = {'BothSex':biomeDF}\n",
    "\n",
    "tempD3 = {}\n",
    "tempD4 = {}\n",
    "for sex in tempD1.keys():\n",
    "    #Select target classes for the classifier\n",
    "    tempDF1 = tempD1[sex]\n",
    "    tempDF1 = tempDF1.loc[tempDF1[yvar+'_class'].isin(list(classD.keys()))]\n",
    "    tempDF2 = tempD2[sex]\n",
    "    tempDF2 = tempDF2.loc[tempDF1.index.tolist()]\n",
    "    \n",
    "    tempD3[sex] = tempDF1\n",
    "    tempD4[sex] = tempDF2\n",
    "    \n",
    "    print(sex+':')\n",
    "    display(tempDF1.describe(include='all'))\n",
    "    display(tempDF2.describe(include='all'))\n",
    "    print('')\n",
    "\n",
    "yDF_B = tempD3['BothSex']\n",
    "xDF_B = tempD4['BothSex']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3-2. Random forest with cross-validation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Gut microbiome-based RF classifiers are generated using the 5-fold iteration scheme (with 5-fold cross-validation).  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Both sex model\n",
    "tempDF1 = xDF_B#Unstandardized input variables\n",
    "tempDF2 = yDF_B#Encoded true class label and testing set label\n",
    "ncvs = 5\n",
    "parameters = {'n_estimators':[int(value) for value in np.geomspace(100, 1000, num=10)],\n",
    "              'max_features':[value for sublist in [['log2', 'sqrt'], np.linspace(0.05, 1.0, num=20)] for value in sublist]}\n",
    "model = RandomForestClassifier(\n",
    "    #n_estimators=100,\n",
    "    criterion='entropy',#Cheged from default (gini) according to Wilmanski, T. et al. Nat. Biotechnol 2019\n",
    "    max_depth=None,\n",
    "    min_samples_split=2,\n",
    "    min_samples_leaf=1,\n",
    "    min_weight_fraction_leaf=0.0,\n",
    "    #max_features='sqrt',\n",
    "    max_leaf_nodes=None,\n",
    "    min_impurity_decrease=0.0,\n",
    "    bootstrap=True,\n",
    "    oob_score=False,#Remain default, because manually calculate scores\n",
    "    n_jobs=-1,#Use all processors = Need to care about the other jobs\n",
    "    random_state=123,#To maintain reproducibility\n",
    "    verbose=0,\n",
    "    warm_start=False,\n",
    "    class_weight='balanced_subsample',#Automatically adjust class weights based on the bootstrap sample\n",
    "    ccp_alpha=0.0,\n",
    "    max_samples=None)\n",
    "modelcv = GridSearchCV(model, parameters, refit=True, cv=ncvs)\n",
    "\n",
    "#Perform random forest\n",
    "featureDF = pd.DataFrame(index=tempDF1.columns).astype('float64')#For feature importance\n",
    "featureDF.index.rename('Variable', inplace=True)\n",
    "predictS_label = pd.Series(name=yvar+'_class_predicted')#For predictions of class label\n",
    "predictS_proba = pd.Series(name=yvar+'_class_predicted-probability')#For predictions of class probability\n",
    "t_start = time.time()\n",
    "for model_k in range(nmodels):\n",
    "    #Prepare training and testing (hold-out) datasets in model k\n",
    "    model_n = 'Model_'+str(model_k+1).zfill(2)\n",
    "    yDF_train = tempDF2.loc[tempDF2['Testing_'+yvar]!=model_n]\n",
    "    yDF_test = tempDF2.loc[tempDF2['Testing_'+yvar]==model_n]\n",
    "    xDF_train = tempDF1.loc[yDF_train.index.tolist()]\n",
    "    xDF_test = tempDF1.loc[yDF_test.index.tolist()]\n",
    "    #Retrieve the true class label\n",
    "    yDF_train = pd.DataFrame(yDF_train[yvar+'_class_code'])#Not Series but DF\n",
    "    #yDF_test = pd.DataFrame(yDF_test[yvar+'_class_code'])#Performance score is calculated later\n",
    "    \n",
    "    #Fitting model with cross-validation using training dataset (i.e., internal training/validation datasets)\n",
    "    modelcv.fit(xDF_train, yDF_train, sample_weight=None)#Weight was set by class_weight parameter\n",
    "    #Check the best hyperparameter set decided by cross validation\n",
    "    print(model_n+':', modelcv.best_params_)\n",
    "    #Extract the best estimator\n",
    "    model_best = modelcv.best_estimator_\n",
    "    #Save feature importance\n",
    "    featureDF[model_n] = model_best.feature_importances_#Impurity-based feature importances\n",
    "    \n",
    "    #Prediction for testing dataset using the fitted model k\n",
    "    ##Label\n",
    "    tempS = pd.Series(model_best.predict(xDF_test),\n",
    "                      index=xDF_test.index, name=predictS_label.name)\n",
    "    predictS_label = pd.concat([predictS_label, tempS], axis=0)\n",
    "    ##Probability\n",
    "    tempS = pd.Series(model_best.predict_proba(xDF_test)[:, 1],#Take only the values for class 1\n",
    "                      index=xDF_test.index, name=predictS_proba.name)\n",
    "    predictS_proba = pd.concat([predictS_proba, tempS], axis=0)\n",
    "t_elapsed = time.time() - t_start\n",
    "print('Elapsed time for '+str(nmodels)+' models of '+str(ncvs)+'-fold CV RF:',\n",
    "      round(t_elapsed//60), 'min', round(t_elapsed%60, 1), 'sec')\n",
    "\n",
    "#Clean prediction DF\n",
    "tempDF = pd.merge(tempDF2[['Testing_'+yvar, yvar+'_class', yvar+'_class_code']], predictS_label,\n",
    "                  left_index=True, right_index=True, how='left')\n",
    "tempDF = pd.merge(tempDF, predictS_proba, left_index=True, right_index=True, how='left')\n",
    "display(tempDF)\n",
    "\n",
    "#Save the cleaned prediction DF\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '221010_Multiomics-BMI-NatMed1stRevision_Microbiome-RFclassifier-ver5_Arivale-wenceslaus_'\n",
    "fileName = yvar+'class-BothSex.tsv'\n",
    "tempDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "measBMI_B = tempDF\n",
    "featureDF_B = featureDF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3-3. Model performance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Evaluation with testing (hold-out) dataset\n",
    "tempD1 = {'BothSex':measBMI_B}\n",
    "tempD2 = {}\n",
    "for sex in tempD1.keys():\n",
    "    tempDF = tempD1[sex]\n",
    "    tempDF1 = pd.DataFrame(index=pd.Index(['Sensitivity', 'Specificity', 'Precision', 'AUC-ROC', 'AUC-PR', 'AP'],\n",
    "                                          name='Metric'))\n",
    "    for model_k in range(nmodels):\n",
    "        model_n = 'Model_'+str(model_k+1).zfill(2)\n",
    "        tempS1 = tempDF[yvar+'_class_code'].loc[tempDF['Testing_'+yvar]==model_n]\n",
    "        tempS2 = tempDF[yvar+'_class_predicted'].loc[tempDF['Testing_'+yvar]==model_n]\n",
    "        tempS3 = tempDF[yvar+'_class_predicted-probability'].loc[tempDF['Testing_'+yvar]==model_n]\n",
    "        \n",
    "        tn, fp, fn, tp = confusion_matrix(tempS1, tempS2).ravel()\n",
    "        sensitivity = tp/(tp+fn)#a.k.a. recall, TPR\n",
    "        specificity = tn/(tn+fp)#a.k.a. 1-FPR\n",
    "        precision = tp/(tp+fp)\n",
    "        \n",
    "        #ROC curve\n",
    "        fprA, tprA, thresholdA = roc_curve(tempS1, tempS3, pos_label=1)\n",
    "        auc_roc = auc(fprA, tprA)\n",
    "        #auc_roc = roc_auc_score(tempS1, tempS3)#Same with the above\n",
    "        \n",
    "        #Precision-recall curve\n",
    "        precisionA, recallA, thresholdA = precision_recall_curve(tempS1, tempS3, pos_label=1)\n",
    "        auc_pr = auc(recallA, precisionA)#Linear interpolation\n",
    "        ap = average_precision_score(tempS1, tempS3)#Without interpolation\n",
    "        \n",
    "        tempDF1[model_n] = [sensitivity, specificity, precision, auc_roc, auc_pr, ap]\n",
    "    #Summarize\n",
    "    tempS1 = tempDF1.mean(axis=1)\n",
    "    tempS1.name = 'Mean'\n",
    "    tempS2 = tempDF1.std(axis=1, ddof=1)#Sample standard deviation\n",
    "    tempS2.name = 'SD'\n",
    "    tempS3 = tempDF1.std(axis=1, ddof=1)/np.sqrt(len(tempDF1.columns))\n",
    "    tempS3.name = 'SEM'\n",
    "    tempDF = pd.concat([tempDF1, tempS1, tempS2, tempS3], axis=1)\n",
    "    \n",
    "    #Save the cleaned DF\n",
    "    fileDir = './ExportData/'\n",
    "    ipynbName = '221010_Multiomics-BMI-NatMed1stRevision_Microbiome-RFclassifier-ver5_Arivale-wenceslaus_'\n",
    "    fileName = yvar+'class-'+sex+'-performance.tsv'\n",
    "    tempDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "    \n",
    "    tempD2[sex] = tempDF\n",
    "    print(sex)\n",
    "    display(tempDF)\n",
    "    print('')\n",
    "\n",
    "measBMI_B_metrics = tempD2['BothSex']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3-4. Clean feature importance dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempD1 = {'BothSex':featureDF_B}\n",
    "for sex in tempD1.keys():\n",
    "    tempDF = tempD1[sex]\n",
    "    #Summarize\n",
    "    tempL1 = []\n",
    "    tempL2 = []\n",
    "    for row_n in tempDF.index.tolist():\n",
    "        tempL1.append(tempDF.loc[row_n].mean())\n",
    "        tempL2.append(tempDF.loc[row_n].std(ddof=1))#Sample standard deviation\n",
    "    tempDF['Mean'] = tempL1\n",
    "    tempDF['SD'] = tempL2\n",
    "    #Save the cleaned DF\n",
    "    fileDir = './ExportData/'\n",
    "    ipynbName = '221010_Multiomics-BMI-NatMed1stRevision_Microbiome-RFclassifier-ver5_Arivale-wenceslaus_'\n",
    "    fileName = yvar+'class-'+sex+'-feature-importance.tsv'\n",
    "    tempDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "    \n",
    "    #Check\n",
    "    print(sex+':')\n",
    "    print(' - Variables:', len(tempDF))\n",
    "    tempDF1 = tempDF.loc[tempDF['Mean']>0.01]\n",
    "    print(' - Variables with the mean of feature importances > 0.01:', len(tempDF1),\n",
    "          '(', len(tempDF1)/len(tempDF)*100, '%)')\n",
    "    tempDF1 = tempDF.loc[tempDF['Mean']>0.05]\n",
    "    print(' - Variables with the mean of feature importances > 0.05:', len(tempDF1),\n",
    "          '(', len(tempDF1)/len(tempDF)*100, '%)')\n",
    "    tempDF1 = tempDF.sort_values(by='Mean', ascending=False)\n",
    "    display(tempDF1.iloc[:30])\n",
    "    print('')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. MetBMI classifier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "yvar = 'MetBMI'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4-1. Select the target class participants"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "tempD1 = {'BothSex':bmiDF}\n",
    "tempD2 = {'BothSex':biomeDF}\n",
    "\n",
    "tempD3 = {}\n",
    "tempD4 = {}\n",
    "for sex in tempD1.keys():\n",
    "    #Select target classes for the classifier\n",
    "    tempDF1 = tempD1[sex]\n",
    "    tempDF1 = tempDF1.loc[tempDF1[yvar+'_class'].isin(list(classD.keys()))]\n",
    "    tempDF2 = tempD2[sex]\n",
    "    tempDF2 = tempDF2.loc[tempDF1.index.tolist()]\n",
    "    \n",
    "    tempD3[sex] = tempDF1\n",
    "    tempD4[sex] = tempDF2\n",
    "    \n",
    "    print(sex+':')\n",
    "    display(tempDF1.describe(include='all'))\n",
    "    display(tempDF2.describe(include='all'))\n",
    "    print('')\n",
    "\n",
    "yDF_B = tempD3['BothSex']\n",
    "xDF_B = tempD4['BothSex']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4-2. Random forest with cross-validation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Gut microbiome-based RF classifiers are generated using the 5-fold iteration scheme (with 5-fold cross-validation).  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Both sex model\n",
    "tempDF1 = xDF_B#Unstandardized input variables\n",
    "tempDF2 = yDF_B#Encoded true class label and testing set label\n",
    "ncvs = 5\n",
    "parameters = {'n_estimators':[int(value) for value in np.geomspace(100, 1000, num=10)],\n",
    "              'max_features':[value for sublist in [['log2', 'sqrt'], np.linspace(0.05, 1.0, num=20)] for value in sublist]}\n",
    "model = RandomForestClassifier(\n",
    "    #n_estimators=100,\n",
    "    criterion='entropy',#Cheged from default (gini) according to Wilmanski, T. et al. Nat. Biotechnol 2019\n",
    "    max_depth=None,\n",
    "    min_samples_split=2,\n",
    "    min_samples_leaf=1,\n",
    "    min_weight_fraction_leaf=0.0,\n",
    "    #max_features='sqrt',\n",
    "    max_leaf_nodes=None,\n",
    "    min_impurity_decrease=0.0,\n",
    "    bootstrap=True,\n",
    "    oob_score=False,#Remain default, because manually calculate scores\n",
    "    n_jobs=-1,#Use all processors = Need to care about the other jobs\n",
    "    random_state=123,#To maintain reproducibility\n",
    "    verbose=0,\n",
    "    warm_start=False,\n",
    "    class_weight='balanced_subsample',#Automatically adjust class weights based on the bootstrap sample\n",
    "    ccp_alpha=0.0,\n",
    "    max_samples=None)\n",
    "modelcv = GridSearchCV(model, parameters, refit=True, cv=ncvs)\n",
    "\n",
    "#Perform random forest\n",
    "featureDF = pd.DataFrame(index=tempDF1.columns).astype('float64')#For feature importance\n",
    "featureDF.index.rename('Variable', inplace=True)\n",
    "predictS_label = pd.Series(name=yvar+'_class_predicted')#For predictions of class label\n",
    "predictS_proba = pd.Series(name=yvar+'_class_predicted-probability')#For predictions of class probability\n",
    "t_start = time.time()\n",
    "for model_k in range(nmodels):\n",
    "    #Prepare training and testing (hold-out) datasets in model k\n",
    "    model_n = 'Model_'+str(model_k+1).zfill(2)\n",
    "    yDF_train = tempDF2.loc[tempDF2['Testing_'+yvar]!=model_n]\n",
    "    yDF_test = tempDF2.loc[tempDF2['Testing_'+yvar]==model_n]\n",
    "    xDF_train = tempDF1.loc[yDF_train.index.tolist()]\n",
    "    xDF_test = tempDF1.loc[yDF_test.index.tolist()]\n",
    "    #Retrieve the true class label\n",
    "    yDF_train = pd.DataFrame(yDF_train[yvar+'_class_code'])#Not Series but DF\n",
    "    #yDF_test = pd.DataFrame(yDF_test[yvar+'_class_code'])#Performance score is calculated later\n",
    "    \n",
    "    #Fitting model with cross-validation using training dataset (i.e., internal training/validation datasets)\n",
    "    modelcv.fit(xDF_train, yDF_train, sample_weight=None)#Weight was set by class_weight parameter\n",
    "    #Check the best hyperparameter set decided by cross validation\n",
    "    print(model_n+':', modelcv.best_params_)\n",
    "    #Extract the best estimator\n",
    "    model_best = modelcv.best_estimator_\n",
    "    #Save feature importance\n",
    "    featureDF[model_n] = model_best.feature_importances_#Impurity-based feature importances\n",
    "    \n",
    "    #Prediction for testing dataset using the fitted model k\n",
    "    ##Label\n",
    "    tempS = pd.Series(model_best.predict(xDF_test),\n",
    "                      index=xDF_test.index, name=predictS_label.name)\n",
    "    predictS_label = pd.concat([predictS_label, tempS], axis=0)\n",
    "    ##Probability\n",
    "    tempS = pd.Series(model_best.predict_proba(xDF_test)[:, 1],#Take only the values for class 1\n",
    "                      index=xDF_test.index, name=predictS_proba.name)\n",
    "    predictS_proba = pd.concat([predictS_proba, tempS], axis=0)\n",
    "t_elapsed = time.time() - t_start\n",
    "print('Elapsed time for '+str(nmodels)+' models of '+str(ncvs)+'-fold CV RF:',\n",
    "      round(t_elapsed//60), 'min', round(t_elapsed%60, 1), 'sec')\n",
    "\n",
    "#Clean prediction DF\n",
    "tempDF = pd.merge(tempDF2[['Testing_'+yvar, yvar+'_class', yvar+'_class_code']], predictS_label,\n",
    "                  left_index=True, right_index=True, how='left')\n",
    "tempDF = pd.merge(tempDF, predictS_proba, left_index=True, right_index=True, how='left')\n",
    "display(tempDF)\n",
    "\n",
    "#Save the cleaned prediction DF\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '221010_Multiomics-BMI-NatMed1stRevision_Microbiome-RFclassifier-ver5_Arivale-wenceslaus_'\n",
    "fileName = yvar+'class-BothSex.tsv'\n",
    "tempDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "metBMI_B = tempDF\n",
    "featureDF_B = featureDF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4-3. Model performance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Evaluation with testing (hold-out) dataset\n",
    "tempD1 = {'BothSex':metBMI_B}\n",
    "tempD2 = {}\n",
    "for sex in tempD1.keys():\n",
    "    tempDF = tempD1[sex]\n",
    "    tempDF1 = pd.DataFrame(index=pd.Index(['Sensitivity', 'Specificity', 'Precision', 'AUC-ROC', 'AUC-PR', 'AP'],\n",
    "                                          name='Metric'))\n",
    "    for model_k in range(nmodels):\n",
    "        model_n = 'Model_'+str(model_k+1).zfill(2)\n",
    "        tempS1 = tempDF[yvar+'_class_code'].loc[tempDF['Testing_'+yvar]==model_n]\n",
    "        tempS2 = tempDF[yvar+'_class_predicted'].loc[tempDF['Testing_'+yvar]==model_n]\n",
    "        tempS3 = tempDF[yvar+'_class_predicted-probability'].loc[tempDF['Testing_'+yvar]==model_n]\n",
    "        \n",
    "        tn, fp, fn, tp = confusion_matrix(tempS1, tempS2).ravel()\n",
    "        sensitivity = tp/(tp+fn)#a.k.a. recall, TPR\n",
    "        specificity = tn/(tn+fp)#a.k.a. 1-FPR\n",
    "        precision = tp/(tp+fp)\n",
    "        \n",
    "        #ROC curve\n",
    "        fprA, tprA, thresholdA = roc_curve(tempS1, tempS3, pos_label=1)\n",
    "        auc_roc = auc(fprA, tprA)\n",
    "        #auc_roc = roc_auc_score(tempS1, tempS3)#Same with the above\n",
    "        \n",
    "        #Precision-recall curve\n",
    "        precisionA, recallA, thresholdA = precision_recall_curve(tempS1, tempS3, pos_label=1)\n",
    "        auc_pr = auc(recallA, precisionA)#Linear interpolation\n",
    "        ap = average_precision_score(tempS1, tempS3)#Without interpolation\n",
    "        \n",
    "        tempDF1[model_n] = [sensitivity, specificity, precision, auc_roc, auc_pr, ap]\n",
    "    #Summarize\n",
    "    tempS1 = tempDF1.mean(axis=1)\n",
    "    tempS1.name = 'Mean'\n",
    "    tempS2 = tempDF1.std(axis=1, ddof=1)#Sample standard deviation\n",
    "    tempS2.name = 'SD'\n",
    "    tempS3 = tempDF1.std(axis=1, ddof=1)/np.sqrt(len(tempDF1.columns))\n",
    "    tempS3.name = 'SEM'\n",
    "    tempDF = pd.concat([tempDF1, tempS1, tempS2, tempS3], axis=1)\n",
    "    \n",
    "    #Save the cleaned DF\n",
    "    fileDir = './ExportData/'\n",
    "    ipynbName = '221010_Multiomics-BMI-NatMed1stRevision_Microbiome-RFclassifier-ver5_Arivale-wenceslaus_'\n",
    "    fileName = yvar+'class-'+sex+'-performance.tsv'\n",
    "    tempDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "    \n",
    "    tempD2[sex] = tempDF\n",
    "    print(sex)\n",
    "    display(tempDF)\n",
    "    print('')\n",
    "\n",
    "metBMI_B_metrics = tempD2['BothSex']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4-4. Clean feature importance dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempD1 = {'BothSex':featureDF_B}\n",
    "for sex in tempD1.keys():\n",
    "    tempDF = tempD1[sex]\n",
    "    #Summarize\n",
    "    tempL1 = []\n",
    "    tempL2 = []\n",
    "    for row_n in tempDF.index.tolist():\n",
    "        tempL1.append(tempDF.loc[row_n].mean())\n",
    "        tempL2.append(tempDF.loc[row_n].std(ddof=1))#Sample standard deviation\n",
    "    tempDF['Mean'] = tempL1\n",
    "    tempDF['SD'] = tempL2\n",
    "    #Save the cleaned DF\n",
    "    fileDir = './ExportData/'\n",
    "    ipynbName = '221010_Multiomics-BMI-NatMed1stRevision_Microbiome-RFclassifier-ver5_Arivale-wenceslaus_'\n",
    "    fileName = yvar+'class-'+sex+'-feature-importance.tsv'\n",
    "    tempDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "    \n",
    "    #Check\n",
    "    print(sex+':')\n",
    "    print(' - Variables:', len(tempDF))\n",
    "    tempDF1 = tempDF.loc[tempDF['Mean']>0.01]\n",
    "    print(' - Variables with the mean of feature importances > 0.01:', len(tempDF1),\n",
    "          '(', len(tempDF1)/len(tempDF)*100, '%)')\n",
    "    tempDF1 = tempDF.loc[tempDF['Mean']>0.05]\n",
    "    print(' - Variables with the mean of feature importances > 0.05:', len(tempDF1),\n",
    "          '(', len(tempDF1)/len(tempDF)*100, '%)')\n",
    "    tempDF1 = tempDF.sort_values(by='Mean', ascending=False)\n",
    "    display(tempDF1.iloc[:30])\n",
    "    print('')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# — Move to the R sub-notebook —"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Comparison b/w classifiers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5-1. ROC curve"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5-1-1. Per model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Refer to Examples of scikit-learn website.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempD1 = {'BMI':measBMI_B, 'MetBMI':metBMI_B}\n",
    "tempD2 = {'BMI':measBMI_B_metrics, 'MetBMI':metBMI_B_metrics}\n",
    "tempD3 = {'BMI':'k', 'MetBMI':'b'}\n",
    "\n",
    "#Prepare each ROC curve\n",
    "mean_fpr = np.linspace(0, 1, 100)#X-coordinate\n",
    "tempD = {}\n",
    "for yvar in tempD1.keys():\n",
    "    tempDF = tempD1[yvar]\n",
    "    tempL = tempDF['Testing_'+yvar].unique().tolist()\n",
    "    tprs = []\n",
    "    for model_n in tempL:\n",
    "        tempS1 = tempDF[yvar+'_class_code'].loc[tempDF['Testing_'+yvar]==model_n]\n",
    "        tempS2 = tempDF[yvar+'_class_predicted-probability'].loc[tempDF['Testing_'+yvar]==model_n]\n",
    "        fprA, tprA, thresholdA = roc_curve(tempS1, tempS2, pos_label=1)\n",
    "        interp_tpr = np.interp(mean_fpr, fprA, tprA)\n",
    "        interp_tpr[0] = 0.0#Left border\n",
    "        interp_tpr[-1] = 1.0#Right border\n",
    "        tprs.append(interp_tpr)\n",
    "    tempD[yvar] = tprs\n",
    "\n",
    "#Visualization\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(5, 5))\n",
    "plt.setp(ax, xlim=(0.0, 1.01), xticks=np.arange(0, 1.1, 0.2))\n",
    "plt.setp(ax, ylim=(0.0, 1.01), yticks=np.arange(0, 1.1, 0.2))\n",
    "for yvar in tempD.keys():\n",
    "    tprs = tempD[yvar]\n",
    "    print(yvar+': n =', len(tprs), 'models')#Check length for the following SEM calculation (just in case)\n",
    "    #Prepare mean and SEM of AUC\n",
    "    tempS = tempD2[yvar].loc['AUC-ROC']\n",
    "    mean = tempS.loc['Mean']\n",
    "    mean_text = str(Decimal(mean).quantize(Decimal('0.01'), rounding=ROUND_HALF_UP))\n",
    "    sem = tempS.loc['SEM']\n",
    "    sem_text = str(Decimal(sem).quantize(Decimal('0.01'), rounding=ROUND_HALF_UP))\n",
    "    display(tempS)\n",
    "    #Mean line\n",
    "    mean_tpr = np.mean(tprs, axis=0)\n",
    "    print(' -> Cf. AUC of the mean ROC curve in plot:', auc(mean_fpr, mean_tpr))\n",
    "    ax.plot(mean_fpr, mean_tpr, color=tempD3[yvar], lw=2,\n",
    "            label=yvar+' class\\n(AUC = '+mean_text+' ± '+sem_text+')')\n",
    "    #SEM range\n",
    "    sem_tpr = np.std(tprs, axis=0, ddof=1)/np.sqrt(len(tprs))\n",
    "    ax.fill_between(mean_fpr, mean_tpr-sem_tpr, mean_tpr+sem_tpr,\n",
    "                    color=tempD3[yvar], alpha=0.2)\n",
    "##Random classification line\n",
    "ax.plot([0, 1], [0, 1], linestyle=\"--\", lw=2, color='r', zorder=0)\n",
    "sns.despine()\n",
    "plt.ylabel('True positive rate')\n",
    "plt.xlabel('False postive rate')\n",
    "plt.legend(fontsize='small', title='Obesity classifier', title_fontsize='medium',\n",
    "           bbox_to_anchor=(1, 0), loc='lower right', borderaxespad=0.25,\n",
    "           handlelength=1.5, handletextpad=0.5)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5-1-2. Overall population"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempD1 = {'BMI':measBMI_B, 'MetBMI':metBMI_B}\n",
    "tempD2 = {'BMI':measBMI_B_metrics, 'MetBMI':metBMI_B_metrics}\n",
    "tempD3 = {'BMI':'k', 'MetBMI':'b'}\n",
    "\n",
    "#Prepare overall ROC curve\n",
    "fprs = np.linspace(0, 1, 100)#X-coordinate\n",
    "tempD = {}\n",
    "for yvar in tempD1.keys():\n",
    "    tempDF = tempD1[yvar]\n",
    "    tempS1 = tempDF[yvar+'_class_code']\n",
    "    tempS2 = tempDF[yvar+'_class_predicted-probability']\n",
    "    fprA, tprA, thresholdA = roc_curve(tempS1, tempS2, pos_label=1)\n",
    "    tprs = np.interp(fprs, fprA, tprA)\n",
    "    tprs[0] = 0.0#Left border\n",
    "    tprs[-1] = 1.0#Right border\n",
    "    tempD[yvar] = tprs\n",
    "\n",
    "#Visualization\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(5, 5))\n",
    "plt.setp(ax, xlim=(0.0, 1.01), xticks=np.arange(0, 1.1, 0.2))\n",
    "plt.setp(ax, ylim=(0.0, 1.01), yticks=np.arange(0, 1.1, 0.2))\n",
    "for yvar in tempD.keys():\n",
    "    tprs = tempD[yvar]\n",
    "    print(yvar+':', len(tprs), 'interpolated points')#Check just in case\n",
    "    #Prepare AUC\n",
    "    auc_roc = auc(fprs, tprs)\n",
    "    auc_text = str(Decimal(auc_roc).quantize(Decimal('0.001'), rounding=ROUND_HALF_UP))\n",
    "    print(' - AUC:', auc_roc)\n",
    "    tempS = tempD2[yvar].loc['AUC-ROC']\n",
    "    display(tempS)\n",
    "    #ROC curve\n",
    "    ax.plot(fprs, tprs, color=tempD3[yvar], lw=2,\n",
    "            label=yvar+' class\\n(AUC = '+auc_text+')')\n",
    "##Random classification line\n",
    "ax.plot([0, 1], [0, 1], linestyle=\"--\", lw=2, color='r', zorder=0)\n",
    "sns.despine()\n",
    "plt.ylabel('True positive rate')\n",
    "plt.xlabel('False postive rate')\n",
    "plt.legend(fontsize='small', title='Obesity classifier', title_fontsize='medium',\n",
    "           bbox_to_anchor=(1, 0), loc='lower right', borderaxespad=0.25,\n",
    "           handlelength=1.5, handletextpad=0.5)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5-1-3. Overall population (pROC package)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Because the above interpolated ROC curve is not completely same with the one used in DeLong's test, the used ROC curve is imported.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempD1 = {'BMI':'k', 'MetBMI':'b'}\n",
    "cohort = 'Arivale'\n",
    "\n",
    "#Prepare the overall ROC curve used in DeLong test\n",
    "tempD2 = {}\n",
    "for yvar in tempD1.keys():\n",
    "    #Import\n",
    "    fileDir = './ExportData/'\n",
    "    ipynbName = '221010_Multiomics-BMI-NatMed1stRevision_Microbiome-RFclassifier-DeLong-ver5_Arivale-wenceslaus_'\n",
    "    fileName = yvar+'-ROC-curve.tsv'\n",
    "    tempDF = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t')\n",
    "    #Calculate TPR and FPR\n",
    "    tempDF['TPR'] = tempDF['Sensitivity']\n",
    "    tempDF['FPR'] = 1 - tempDF['Specificity']\n",
    "    print(yvar)\n",
    "    display(tempDF)\n",
    "    tempD2[yvar] = tempDF\n",
    "\n",
    "#Prepare test summary\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '221010_Multiomics-BMI-NatMed1stRevision_Microbiome-RFclassifier-DeLong-ver5_Arivale-wenceslaus_'\n",
    "fileName = 'result-summary.tsv'\n",
    "tempDF = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t')\n",
    "tempDF = tempDF.set_index('Variable')\n",
    "display(tempDF)\n",
    "tempD3 = {}\n",
    "for yvar in tempD1.keys():\n",
    "    tempD3[yvar] = tempDF['Estimate_'+yvar].iloc[0]\n",
    "pval = tempDF['Pval'].iloc[0]\n",
    "\n",
    "#Visualization\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(5, 5))\n",
    "plt.setp(ax, xlim=(0.0, 1.01), xticks=np.arange(0, 1.1, 0.2))\n",
    "plt.setp(ax, ylim=(0.0, 1.01), yticks=np.arange(0, 1.1, 0.2))\n",
    "for yvar in tempD1.keys():\n",
    "    tempDF = tempD2[yvar]\n",
    "    #Prepare AUC\n",
    "    auc_roc = tempD3[yvar]\n",
    "    auc_text = str(Decimal(auc_roc).quantize(Decimal('0.001'), rounding=ROUND_HALF_UP))\n",
    "    #ROC curve\n",
    "    ax.plot(tempDF['FPR'], tempDF['TPR'], color=tempD1[yvar], lw=2,\n",
    "            label=yvar+' class\\n(AUC = '+auc_text+')')\n",
    "##Random classification line\n",
    "ax.plot([0, 1], [0, 1], linestyle=\"--\", lw=2, color='r', zorder=0)\n",
    "##Add annotation line\n",
    "xcoord_0 = 0.4#Manually adjusted\n",
    "xcoord_1 = 0.425#Manually adjusted\n",
    "ycoord_0 = 0.090#Manually adjusted\n",
    "ycoord_1 = 0.230#Manually adjusted\n",
    "ax.plot([xcoord_1, xcoord_0, xcoord_0, xcoord_1],\n",
    "        [ycoord_0, ycoord_0, ycoord_1, ycoord_1],\n",
    "        lw=1.5, c='k')\n",
    "##Add P-value annotation\n",
    "if pval<0.001:\n",
    "    label = '***'\n",
    "elif pval<0.01:\n",
    "    label = '**'\n",
    "elif pval<0.05:\n",
    "    label = '*'\n",
    "else:\n",
    "    pval_text = str(Decimal(pval).quantize(Decimal('0.01'), rounding=ROUND_HALF_UP))\n",
    "    label = r'$P$'+' = '+pval_text\n",
    "if label in ['***', '**', '*']:\n",
    "    text_offset = -0.015\n",
    "    text_size = 'medium'\n",
    "else:\n",
    "    text_offset = 0.0\n",
    "    text_size = 'x-small'\n",
    "##Add axis title\n",
    "ax.set_title(cohort+' cohort', {'fontsize':'large'})\n",
    "ax.annotate(label, xy=(2*xcoord_0-xcoord_1, (ycoord_0+ycoord_1)/2+text_offset),\n",
    "            horizontalalignment='right', verticalalignment='center',\n",
    "            fontsize=text_size, color='k')\n",
    "sns.despine()\n",
    "plt.ylabel('True positive rate')\n",
    "plt.xlabel('False postive rate')\n",
    "plt.legend(fontsize='small', title='Obesity classifier', title_fontsize='medium',\n",
    "           bbox_to_anchor=(1, 0), loc='lower right', borderaxespad=0.25,\n",
    "           handlelength=1.5, handletextpad=0.5)\n",
    "##Save\n",
    "fileDir = './ExportFigures/'\n",
    "ipynbName = '221010_Multiomics-BMI-NatMed1stRevision_Microbiome-RFclassifier-ver5_Arivale-wenceslaus_'\n",
    "fileName = 'ROC-curve.tif'\n",
    "plt.gcf().savefig(fileDir+ipynbName+fileName, dpi=300, bbox_inches='tight', pad_inches=0.04,\n",
    "                  pil_kwargs={'compression':'tiff_lzw'})\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5-2. Precision–recall curve"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5-2-1. Per model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Refer to Examples of scikit-learn website.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempD1 = {'BMI':measBMI_B, 'MetBMI':metBMI_B}\n",
    "tempD2 = {'BMI':measBMI_B_metrics, 'MetBMI':metBMI_B_metrics}\n",
    "tempD3 = {'BMI':'k', 'MetBMI':'b'}\n",
    "\n",
    "#Prepare each PR curve\n",
    "mean_recall = np.linspace(0, 1, 100)#X-coordinate\n",
    "tempD = {}\n",
    "for yvar in tempD1.keys():\n",
    "    tempDF = tempD1[yvar]\n",
    "    tempL = tempDF['Testing_'+yvar].unique().tolist()\n",
    "    precisions = []\n",
    "    for model_n in tempL:\n",
    "        tempS1 = tempDF[yvar+'_class_code'].loc[tempDF['Testing_'+yvar]==model_n]\n",
    "        tempS2 = tempDF[yvar+'_class_predicted-probability'].loc[tempDF['Testing_'+yvar]==model_n]\n",
    "        precisionA, recallA, thresholdA = precision_recall_curve(tempS1, tempS2, pos_label=1)\n",
    "        precisionA = np.flip(precisionA)#To make the first element 1\n",
    "        recallA = np.flip(recallA)#To make the first element 0\n",
    "        interp_precision = np.interp(mean_recall, recallA, precisionA)\n",
    "        interp_precision[0] = 1.0#Left border\n",
    "        random_precision = len(tempS1.loc[tempS1==1])/len(tempS1)#No-skill classifier precision\n",
    "        interp_precision[-1] = random_precision#Right border\n",
    "        precisions.append(interp_precision)\n",
    "    tempD[yvar] = precisions\n",
    "\n",
    "#Visualization\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(5, 5))\n",
    "plt.setp(ax, xlim=(0.0, 1.01), xticks=np.arange(0, 1.1, 0.2))\n",
    "plt.setp(ax, ylim=(0.0, 1.01), yticks=np.arange(0, 1.1, 0.2))\n",
    "for yvar in tempD.keys():\n",
    "    precisions = tempD[yvar]\n",
    "    print(yvar+': n =', len(precisions), 'models')#Check length for the following SEM calculation (just in case)\n",
    "    #Prepare mean and SEM of AUC\n",
    "    tempS = tempD2[yvar].loc['AUC-PR']\n",
    "    mean = tempS.loc['Mean']\n",
    "    mean_text = str(Decimal(mean).quantize(Decimal('0.01'), rounding=ROUND_HALF_UP))\n",
    "    sem = tempS.loc['SEM']\n",
    "    sem_text = str(Decimal(sem).quantize(Decimal('0.01'), rounding=ROUND_HALF_UP))\n",
    "    display(tempS)\n",
    "    #Mean line\n",
    "    mean_precision = np.mean(precisions, axis=0)\n",
    "    print(' -> Cf. AUC of the mean PR curve in plot:', auc(mean_recall, mean_precision))\n",
    "    ax.plot(mean_recall, mean_precision, color=tempD3[yvar], lw=2,\n",
    "            label=yvar+' class\\n(AUC = '+mean_text+' ± '+sem_text+')')\n",
    "    #SEM range\n",
    "    sem_precision = np.std(precisions, axis=0, ddof=1)/np.sqrt(len(precisions))\n",
    "    ax.fill_between(mean_recall, mean_precision-sem_precision, mean_precision+sem_precision,\n",
    "                    color=tempD3[yvar], alpha=0.2)\n",
    "    #Random classification line\n",
    "    tempDF = tempD1[yvar]\n",
    "    tempS = tempDF[yvar+'_class_code']\n",
    "    random_precision = len(tempS.loc[tempS==1])/len(tempS)#No-skill classifier precision\n",
    "    print(' - No-skill classifier precision (overall):', random_precision)\n",
    "    print('   -> Cf. the mean among testing sets:', mean_precision[-1])\n",
    "    ax.axhline(y=random_precision, linestyle=\"--\", lw=2, color=tempD3[yvar], zorder=0)\n",
    "sns.despine()\n",
    "plt.ylabel('Precision')\n",
    "plt.xlabel('Recall')\n",
    "plt.legend(fontsize='small', title='Obesity classifier', title_fontsize='medium',\n",
    "           bbox_to_anchor=(1, 0), loc='lower right', borderaxespad=0.25,\n",
    "           handlelength=1.5, handletextpad=0.5)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 5-2-2. Overall population"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempD1 = {'BMI':measBMI_B, 'MetBMI':metBMI_B}\n",
    "tempD2 = {'BMI':measBMI_B_metrics, 'MetBMI':metBMI_B_metrics}\n",
    "tempD3 = {'BMI':'k', 'MetBMI':'b'}\n",
    "\n",
    "#Prepare overall PR curve\n",
    "recalls = np.linspace(0, 1, 100)#X-coordinate\n",
    "tempD = {}\n",
    "for yvar in tempD1.keys():\n",
    "    tempDF = tempD1[yvar]\n",
    "    tempS1 = tempDF[yvar+'_class_code']\n",
    "    tempS2 = tempDF[yvar+'_class_predicted-probability']\n",
    "    precisionA, recallA, thresholdA = precision_recall_curve(tempS1, tempS2, pos_label=1)\n",
    "    precisionA = np.flip(precisionA)#To make the first element 1\n",
    "    recallA = np.flip(recallA)#To make the first element 0\n",
    "    precisions = np.interp(recalls, recallA, precisionA)\n",
    "    precisions[0] = 1.0#Left border\n",
    "    random_precision = len(tempS1.loc[tempS1==1])/len(tempS1)#No-skill classifier precision\n",
    "    precisions[-1] = random_precision#Right border\n",
    "    tempD[yvar] = precisions\n",
    "\n",
    "#Visualization\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(5, 5))\n",
    "plt.setp(ax, xlim=(0.0, 1.01), xticks=np.arange(0, 1.1, 0.2))\n",
    "plt.setp(ax, ylim=(0.0, 1.01), yticks=np.arange(0, 1.1, 0.2))\n",
    "for yvar in tempD.keys():\n",
    "    precisions = tempD[yvar]\n",
    "    print(yvar+':', len(precisions), 'interpolated points')#Check just in case\n",
    "    #Prepare AUC\n",
    "    auc_pr = auc(recalls, precisions)\n",
    "    auc_text = str(Decimal(auc_pr).quantize(Decimal('0.001'), rounding=ROUND_HALF_UP))\n",
    "    print(' - AUC:', auc_pr)\n",
    "    tempS = tempD2[yvar].loc['AUC-PR']\n",
    "    display(tempS)\n",
    "    #PR curve\n",
    "    ax.plot(recalls, precisions, color=tempD3[yvar], lw=2,\n",
    "            label=yvar+' class\\n(AUC = '+auc_text+')')\n",
    "    #Random classification line\n",
    "    tempDF = tempD1[yvar]\n",
    "    tempS = tempDF[yvar+'_class_code']\n",
    "    random_precision = len(tempS.loc[tempS==1])/len(tempS)#No-skill classifier precision\n",
    "    print(' - No-skill classifier precision (overall):', random_precision)\n",
    "    ax.axhline(y=random_precision, linestyle=\"--\", lw=2, color=tempD3[yvar], zorder=0)\n",
    "sns.despine()\n",
    "plt.ylabel('Precision')\n",
    "plt.xlabel('Recall')\n",
    "plt.legend(fontsize='small', title='Obesity classifier', title_fontsize='medium',\n",
    "           bbox_to_anchor=(1, 0), loc='lower right', borderaxespad=0.25,\n",
    "           handlelength=1.5, handletextpad=0.5)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5-3. Prediction metrics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cohort = 'Arivale'\n",
    "tempD1 = {'BMI class':measBMI_B_metrics, 'MetBMI class':metBMI_B_metrics}\n",
    "tempD2 = {'BMI class':'0.5', 'MetBMI class':'b'}\n",
    "tempD3 = {'AUC-ROC':'AUC (ROC)', 'Sensitivity':'Sensitivity',\n",
    "          'Specificity':'Specificity', 'AP':'Precision'}\n",
    "\n",
    "#Prepare DF\n",
    "tempL = []\n",
    "for classifier in tempD1.keys():\n",
    "    tempDF = tempD1[classifier].T\n",
    "    tempDF = tempDF.loc[tempDF.index.str.contains('Model_')]\n",
    "    tempDF.index.rename('Model', inplace=True)\n",
    "    tempDF = tempDF[list(tempD3.keys())]\n",
    "    tempDF.columns = tempDF.columns.map(tempD3)\n",
    "    tempDF['Classifier'] = classifier\n",
    "    tempL.append(tempDF)\n",
    "    \n",
    "    print(classifier)\n",
    "    display(tempDF.describe(include='all'))\n",
    "tempDF = pd.concat(tempL, axis=0)\n",
    "\n",
    "#Statistical tests\n",
    "control = list(tempD1.keys())[0]\n",
    "contrast = list(tempD1.keys())[1]\n",
    "tempDF1 = pd.DataFrame(columns=['Control', 'Contrast', 'Control_N', 'Contrast_N', 'DoF', 'tStat', 'Pval'])\n",
    "for metric in tempD3.values():\n",
    "    tempS1 = tempDF[metric].loc[tempDF['Classifier']==control]\n",
    "    tempS2 = tempDF[metric].loc[tempDF['Classifier']==contrast]\n",
    "    #Two-sided Welch's t-test\n",
    "    tstat, pval, dof = weightstats.ttest_ind(tempS1, tempS2,\n",
    "                                             alternative='two-sided', usevar='unequal')\n",
    "    size1 = len(tempS1)\n",
    "    size2 = len(tempS2)\n",
    "    tempDF1.loc[metric] = [control, contrast, size1, size2, dof, tstat, pval]\n",
    "tempDF1.index.rename('Metric', inplace=True)\n",
    "display(tempDF1)\n",
    "##Save\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '221010_Multiomics-BMI-NatMed1stRevision_Microbiome-RFclassifier-ver5_Arivale-wenceslaus_'\n",
    "fileName = 'performance-comparison.tsv'\n",
    "tempDF1.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "\n",
    "#Plot\n",
    "tempDF = tempDF.reset_index().melt(var_name='Metric', value_name='Value',\n",
    "                                   id_vars=['Classifier', 'Model'], value_vars=list(tempD3.values()))\n",
    "axis_ymin = 0.0\n",
    "axis_ymax = 1.1\n",
    "ymin = 0.0\n",
    "ymax = 1.0\n",
    "yinter = 0.2\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "plt.figure(figsize=(3.5, 3.5))\n",
    "sns.barplot(data=tempDF, y='Value', x='Metric', order=tempD3.values(),\n",
    "            hue='Classifier', hue_order=tempD2.keys(), dodge=True, palette=tempD2,\n",
    "            ci=95, capsize=0.4/2, errwidth=1.5, errcolor='black', edgecolor='black')\n",
    "p = sns.stripplot(data=tempDF, y='Value', x='Metric', order=tempD3.values(),\n",
    "                  hue='Classifier', hue_order=tempD2.keys(), dodge=True, jitter=0.3/2,\n",
    "                  size=5, edgecolor='black', color='gray', linewidth=1, alpha=0.4)\n",
    "sns.despine()\n",
    "p.set(ylim=(axis_ymin, axis_ymax), yticks=np.arange(ymin, ymax+yinter/10, yinter))\n",
    "plt.setp(p.get_xticklabels(), rotation=70,\n",
    "         horizontalalignment='right', verticalalignment='center', rotation_mode='anchor')\n",
    "##Random classification line\n",
    "p.axhline(y=0.5, linestyle=\"--\", lw=2, color='r', zorder=0)\n",
    "##P-value annotation\n",
    "for x_i, metric in enumerate(tempD3.values()):\n",
    "    #Control\n",
    "    rect_0 = p.patches[x_i]\n",
    "    xcoord_0 = rect_0.get_x() + rect_0.get_width()/2\n",
    "    ycoord_0 = rect_0.get_height()\n",
    "    #Contrast\n",
    "    rect_1 = p.patches[x_i+len(tempD3)]\n",
    "    xcoord_1 = rect_1.get_x() + rect_1.get_width()/2\n",
    "    ycoord_1 = rect_1.get_height()\n",
    "    #Standard point of marker\n",
    "    xcoord = (xcoord_0+xcoord_1)/2\n",
    "    ycoord = tempDF['Value'].loc[tempDF['Metric']==metric].max()\n",
    "    #Add annotation lines\n",
    "    aline_offset = yinter/10\n",
    "    aline_length = yinter/10 + aline_offset\n",
    "    p.plot([xcoord_0, xcoord_0, xcoord_1, xcoord_1],\n",
    "           [ycoord+aline_offset, ycoord+aline_length, ycoord+aline_length, ycoord+aline_offset],\n",
    "           lw=1.5, c='k')\n",
    "    #Retrieve P-value\n",
    "    pval = tempDF1['Pval'].loc[metric]\n",
    "    if pval<0.001:\n",
    "        label = '***'\n",
    "    elif pval<0.01:\n",
    "        label = '**'\n",
    "    elif pval<0.05:\n",
    "        label = '*'\n",
    "    else:\n",
    "        pval_text = str(Decimal(pval).quantize(Decimal('0.01'), rounding=ROUND_HALF_UP))\n",
    "        label = r'$P$'+' = '+pval_text\n",
    "    #Add annotation text\n",
    "    if label in ['***', '**', '*']:\n",
    "        text_offset = yinter/25\n",
    "        text_size = 'medium'\n",
    "    else:\n",
    "        text_offset = yinter/5\n",
    "        text_size = 'x-small'\n",
    "    p.annotate(label, xy=(xcoord, ycoord+text_offset),\n",
    "               horizontalalignment='center', verticalalignment='bottom',\n",
    "               fontsize=text_size, color='k')\n",
    "##Add axis title\n",
    "p.set_title(cohort+' cohort', {'fontsize':'large'})\n",
    "plt.ylabel('Out-of-sample metric')\n",
    "plt.xlabel('')\n",
    "##Modify legend\n",
    "h, l = p.get_legend_handles_labels()\n",
    "h = h[2:]#Remove sns.stripplot legend (legend=False didn't work...)\n",
    "l = l[2:]#Remove sns.stripplot legend (legend=False didn't work...)\n",
    "plt.legend(h, l, fontsize='medium', title='Obesity classifier', title_fontsize='medium',\n",
    "           bbox_to_anchor=(0.5, -0.45), loc='upper center', borderaxespad=0.25,\n",
    "           handlelength=1.5, handletextpad=0.5)\n",
    "##Save\n",
    "fileDir = './ExportFigures/'\n",
    "ipynbName = '221010_Multiomics-BMI-NatMed1stRevision_Microbiome-RFclassifier-ver5_Arivale-wenceslaus_'\n",
    "fileName = 'performance-comparison.tif'\n",
    "plt.gcf().savefig(fileDir+ipynbName+fileName, dpi=300, bbox_inches='tight', pad_inches=0.04,\n",
    "                  pil_kwargs={'compression':'tiff_lzw'})\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# — End of this notebook —"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "arivale-py3 - Python",
   "language": "python",
   "name": "conda-env-arivale-py3-py"
  },
  "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.9.7"
  },
  "toc-showcode": false,
  "toc-showmarkdowntxt": false,
  "toc-showtags": false,
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}