Code11_DeltaBMI_clinical-misclassification.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multiomics BMI Paper — Relationships between ∆BMI and Clinical Definition-based Misclassification"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***by Kengo Watanabe***  \n",
    "\n",
    "This Jupyter Notebook (with Python 3 kernel) assessed differences in ∆BMI (1) between metabolically healthy normal-weight (MHNW) and metabolically unhealthy normal-weight (MUNW) and (2) between metabolically healthy obese (MHO) and metabolically unhealthy obese (MUO), based on a clinical definition (in the baseline Arivale cohort).  \n",
    "\n",
    "Input files:  \n",
    "* Arivale baseline biological BMIs and covariates: 220803_Multiomics-BMI-NatMed1stRevision_DeltaBMI-misclassification_biologicalBMI-baseline-summary-BothSex.tsv  \n",
    "* Arivale baseline metabolic health condition: 220720_Multiomics-BMI-NatMedRevision_Misclassification_metabolic-health-summary.tsv  \n",
    "\n",
    "Output figures and tables:  \n",
    "* Figure 3b  \n",
    "* Table for Supplementary Data 10  \n",
    "\n",
    "Original notebook (memo for my future tracing):  \n",
    "* dalek:\\[JupyterLab HOME\\]/220621_Multiomics-BMI-NatMedRevision/220804_Multiomics-BMI-NatMed1stRevision_DeltaBMI-misclassification_ClinicalDefinition-ver2.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            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        py39hf3d152e_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",
      "matplotlib-venn           0.11.7             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.1o               h166bdaf_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",
      "venn                      0.1.3                    pypi_0    pypi\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",
    "from sklearn.preprocessing import StandardScaler\n",
    "import statsmodels.formula.api as smf\n",
    "from statsmodels.stats import multitest as multi\n",
    "from decimal import Decimal, ROUND_HALF_UP\n",
    "\n",
    "!conda list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Prepare datasets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1-1. ∆BMI and covariates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Import cleaned table for baseline measured and biological BMIs\n",
    "fileDir = './ExportData/'\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",
    "\n",
    "#Clean to handle easier in this notebook\n",
    "tempDF.columns = tempDF.columns.str.replace('Base', '')\n",
    "\n",
    "#Calculate the rate of difference\n",
    "tempL = ['MetBMI', 'ProtBMI', 'ChemBMI', 'CombiBMI']\n",
    "for bbmi in tempL:\n",
    "    tempDF['Delta'+bbmi] = (tempDF[bbmi] - tempDF['BMI']) / tempDF['BMI'] * 100\n",
    "\n",
    "#Select the misclassification and covariates (just for the display in Jupyter notebook)\n",
    "tempL1 = tempDF.loc[:, tempDF.columns.str.contains('Delta')].columns.tolist()\n",
    "tempL2 = tempDF.loc[:, tempDF.columns.str.contains('_class')].columns.tolist()\n",
    "tempL3 = ['BMI', 'Sex', 'Age', 'PC1', 'PC2', 'PC3', 'PC4', 'PC5']\n",
    "tempL = [col_n for sublist in [tempL1, tempL2, tempL3] for col_n in sublist]\n",
    "tempDF = tempDF[tempL]\n",
    "\n",
    "display(tempDF)\n",
    "display(tempDF.describe(include='all'))\n",
    "print('NaN in DF:', tempDF.isnull().to_numpy().sum(axis=None))\n",
    "\n",
    "bmiDF = tempDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF = bmiDF\n",
    "\n",
    "#Check skewness\n",
    "tempDF1 = tempDF.select_dtypes(include=[np.number])\n",
    "tempDF2 = tempDF1.describe()\n",
    "tempDF2.loc['Skewness'] = stats.skew(tempDF1)\n",
    "display(tempDF2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> –> ∆BMI can be surly assumed as normal distribution.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1-2. Metabolic health condition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Import cleaned table for metabolic health condition\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220720_Multiomics-BMI-NatMedRevision_Misclassification_'\n",
    "fileName = 'metabolic-health-summary.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",
    "\n",
    "display(tempDF)\n",
    "display(tempDF.describe(include='all'))\n",
    "print('NaN in DF:', tempDF.isnull().to_numpy().sum(axis=None))\n",
    "\n",
    "metabDF = tempDF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> –> Six participants had NaN for the metabolic health condition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF1 = bmiDF\n",
    "tempDF2 = metabDF\n",
    "\n",
    "#Check\n",
    "tempDF = pd.merge(tempDF1, tempDF2['Metabolically'], left_index=True, right_index=True, how='inner')\n",
    "for bmi_class in ['Normal', 'Obese']:\n",
    "    tempDF3 = tempDF.loc[tempDF['BMI_class']==bmi_class]\n",
    "    print(bmi_class+': n =', len(tempDF3))\n",
    "    display(tempDF3.groupby('Metabolically')['BMI'].describe())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> –> It would be safer to adjust the baseline BMI in the statistical tests.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Regression analysis for ∆BMI"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Of note, because ∆BMI values are assumed as normal distribution, OLS linear regression (i.e., GLM with Gaussian family) can be used simply.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2-1. Perform OLS linear regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Model: ∆BMI ~ b0 + b1\\*C(MetabolicCondition) + b2\\*BMI + b3\\*C(Sex) + b4\\*Age + b5\\*AncestryPCs  \n",
    "> Main aim: Assess the difference in each ∆BMI between the metabolically healthy and unhealthy groups.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF1 = bmiDF\n",
    "tempDF2 = metabDF\n",
    "tempL1 = ['Normal', 'Obese']\n",
    "tempL2 = ['MetBMI', 'ProtBMI', 'ChemBMI', 'CombiBMI']\n",
    "\n",
    "t_start = time.time()\n",
    "tempD1 = {}\n",
    "for bmi_class in tempL1:\n",
    "    #Processing for OLS linear regression\n",
    "    ##Gather all necessary variables into a single DF\n",
    "    tempS = tempDF2['Metabolically']\n",
    "    tempDF = pd.merge(tempDF1, tempS, left_index=True, right_index=True, how='left')\n",
    "    ##Select the target participants\n",
    "    tempDF = tempDF.loc[tempDF['BMI_class']==bmi_class]\n",
    "    ##Drop NaN in the metabolic health condition\n",
    "    tempDF = tempDF.dropna()\n",
    "    ##Z-score transformation\n",
    "    tempDF3 = tempDF.select_dtypes(include=[np.number])\n",
    "    scaler = StandardScaler(copy=True, with_mean=True, with_std=True)\n",
    "    tempA = scaler.fit_transform(tempDF3)#Column direction\n",
    "    tempDF3 = pd.DataFrame(data=tempA, index=tempDF3.index, columns=tempDF3.columns)\n",
    "    ###Recover the categorical variables\n",
    "    tempDF4 = tempDF.select_dtypes(exclude=[np.number])\n",
    "    tempDF = pd.merge(tempDF3, tempDF4, left_index=True, right_index=True, how='left')\n",
    "    ##Add a constant for the intercept\n",
    "    ###–> In statsmodels, a constant is automatically added as well as R!\n",
    "    ##Sort to make bcoef = 0 and 1 for Healthy and Unhealthy\n",
    "    tempDF = tempDF.sort_values(by='Metabolically', ascending=True)\n",
    "    ##One-hot encoding for categorical covariates\n",
    "    ###–> In statsmodels, categorical variables are automatically recognized!\n",
    "    \n",
    "    tempD2 = {}\n",
    "    for bbmi in tempL2:\n",
    "        #OLS linear regression\n",
    "        ##Fit univariate model\n",
    "        formula = 'Delta'+bbmi+' ~ C(Metabolically)'\n",
    "        fit_res1 = smf.ols(formula, data=tempDF).fit()\n",
    "        ##Fit full model\n",
    "        formula = 'Delta'+bbmi+' ~ C(Metabolically)'\\\n",
    "            '+ BMI + C(Sex) + Age + PC1 + PC2 + PC3 + PC4 + PC5'\n",
    "        fit_res2 = smf.ols(formula, data=tempDF).fit()\n",
    "        \n",
    "        #Summarize the result\n",
    "        tempDF3 = pd.DataFrame({'DeltaBMI':[bbmi]})\n",
    "        ##Save the sample size for each group\n",
    "        tempDF3['N'] = len(tempDF)\n",
    "        tempDF3['nHealthy'] = len(tempDF.loc[tempDF['Metabolically']=='Healthy'])\n",
    "        tempDF3['nUnhealthy'] = len(tempDF.loc[tempDF['Metabolically']=='Unhealthy'])\n",
    "        ##Save R2 [%]\n",
    "        tempDF3['UnivarR2'] = fit_res1.rsquared*100\n",
    "        tempDF3['R2'] = fit_res2.rsquared*100\n",
    "        ##Save beta-coefficient of the target variable\n",
    "        tempDF3['Bcoef'] = fit_res2.params['C(Metabolically)[T.Unhealthy]']\n",
    "        tempDF3['BcoefSE'] = fit_res2.bse['C(Metabolically)[T.Unhealthy]']\n",
    "        ##Save t-statistic of the target variable\n",
    "        tempDF3['tStat'] = fit_res2.tvalues['C(Metabolically)[T.Unhealthy]']\n",
    "        ##Save residual degrees of freedom\n",
    "        tempDF3['DoF'] = int(fit_res2.df_resid)\n",
    "        ##Save P-value of the target variable\n",
    "        tempDF3['Pval'] = fit_res2.pvalues['C(Metabolically)[T.Unhealthy]']\n",
    "        \n",
    "        tempD2[bbmi] = tempDF3\n",
    "    \n",
    "    #Clean the results (pd.DataFrame) across bBMIs\n",
    "    tempDF = pd.concat(list(tempD2.values()), axis=0)\n",
    "    tempDF['BMIclass'] = bmi_class\n",
    "    \n",
    "    tempD1[bmi_class] = tempDF\n",
    "t_elapsed = time.time() - t_start\n",
    "print('Elapsed time for',\n",
    "      len(tempL1)*len(tempL2), 'OLS linear regressions (',\n",
    "      len(tempL1), 'BMI classes x',\n",
    "      len(tempL2), 'bBMIs):',\n",
    "      round(t_elapsed//60), 'min', round(t_elapsed%60, 1), 'sec')\n",
    "\n",
    "#Clean the results (pd.DataFrame) across BMI classes\n",
    "tempDF = pd.concat(list(tempD1.values()), axis=0)\n",
    "##Clean the column order by setting index\n",
    "tempDF = tempDF.set_index(['BMIclass', 'DeltaBMI'])\n",
    "\n",
    "#P-value adjustment (across BMI classes within each bBMI) by using Benjamini–Hochberg method\n",
    "tempD = {}\n",
    "for bbmi in tempL2:\n",
    "    tempL = [(bmi_class, bbmi) for bmi_class in tempL1]\n",
    "    tempS = tempDF['Pval'].loc[tempL]\n",
    "    tempA = multi.multipletests(tempS, alpha=0.05, method='fdr_bh',\n",
    "                                is_sorted=False, returnsorted=False)[1]\n",
    "    tempS = pd.Series(tempA, index=tempS.index, name='AdjPval_within')\n",
    "    tempD[bbmi] = tempS\n",
    "tempS = pd.concat(list(tempD.values()), axis=0)\n",
    "tempDF = pd.merge(tempDF, tempS, left_index=True, right_index=True, how='left')\n",
    "\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",
    "\n",
    "display(tempDF)\n",
    "\n",
    "#Save\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220804_Multiomics-BMI-NatMed1stRevision_DeltaBMI-misclassification_ClinicalDefinition-ver2_'\n",
    "fileName = 'result-summary.tsv'\n",
    "tempDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "\n",
    "resDF = tempDF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2-2. Visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempD1 = {'MetBMI':'b', 'ProtBMI':'r', 'ChemBMI':'g', 'CombiBMI':'m'}\n",
    "tempD2 = {'Healthy':'0.8', 'Unhealthy':'crimson'}\n",
    "tempL1 = ['Normal', 'Obese']\n",
    "tempDF1 = bmiDF\n",
    "tempDF2 = metabDF\n",
    "tempDF3 = resDF\n",
    "\n",
    "#Prepare DF\n",
    "tempS = tempDF2['Metabolically']\n",
    "tempDF = pd.merge(tempDF1, tempS, left_index=True, right_index=True, how='left')\n",
    "##Select the target participants\n",
    "tempDF = tempDF.loc[tempDF['BMI_class'].isin(tempL1)]\n",
    "##Drop NaN in the metabolic health condition\n",
    "tempDF = tempDF.dropna()\n",
    "\n",
    "#Check sample size\n",
    "print('N (total):', len(tempDF))\n",
    "print(' - BMI class:', tempDF['BMI_class'].value_counts().sort_index(ascending=True).to_dict())\n",
    "for bmi_class in tempL1:\n",
    "    tempDF1 = tempDF.loc[tempDF['BMI_class']==bmi_class]\n",
    "    print('   - '+bmi_class+' BMI class - Metabolic condition:',\n",
    "          tempDF1['Metabolically'].value_counts().sort_index(ascending=True).to_dict())\n",
    "\n",
    "#Visualization\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "fig, axes = plt.subplots(nrows=1, ncols=len(tempD1),\n",
    "                         figsize=(7.1, 3), sharex=True, sharey=True,\n",
    "                         gridspec_kw={'width_ratios':[1, 1, 1, 1]})\n",
    "axis_ymin = -37.5\n",
    "axis_ymax = 50\n",
    "ymin = -30\n",
    "ymax = 45\n",
    "yinter = 15\n",
    "margin = 0.49\n",
    "#Set shared axis range\n",
    "plt.setp(axes, ylim=(axis_ymin, axis_ymax), yticks=np.arange(ymin, ymax+yinter/10, yinter))\n",
    "plt.setp(axes, xlim=(0-margin, len(tempD2)-1+margin))#To eliminate excess white space\n",
    "for ax_i, ax in enumerate(axes.flat):\n",
    "    bbmi = list(tempD1.keys())[ax_i]\n",
    "    sns.boxplot(data=tempDF, y='Delta'+bbmi, x='BMI_class', order=tempL1,\n",
    "                hue='Metabolically', hue_order=tempD2.keys(), dodge=True, palette=tempD2,\n",
    "                showfliers=False,#flierprops={'marker':'o', 'markerfacecolor':'gray', 'alpha':0.4},\n",
    "                showcaps=True, notch=True, ax=ax)\n",
    "    #Axis setting\n",
    "    if ax_i==0:\n",
    "        plt.setp(ax, xlabel='', ylabel=r'$\\Delta$'+'BMI [% BMI]')\n",
    "    else:\n",
    "        plt.setp(ax.get_yticklabels(), visible=False)\n",
    "        plt.setp(ax, xlabel='', ylabel='')\n",
    "    sns.despine()\n",
    "    plt.setp(ax.get_xticklabels(), rotation=70,\n",
    "             horizontalalignment='right', verticalalignment='center', rotation_mode='anchor')\n",
    "    #P-value annotation\n",
    "    lines = ax.get_lines()#Line2D: [[Q1, Q1-1.5IQR], [Q3, Q3+1.5IQR], [Q1, Q1], [Q3, Q3], [Med, Med], [flier]]\n",
    "    lines_unit = 5 + int(False)#showfliers=False\n",
    "    for class_i in range(len(tempL1)):\n",
    "        #Healthy\n",
    "        whisker_0 = lines[class_i*lines_unit*len(tempD2) + lines_unit*0 + 1]\n",
    "        xcoord_0 = whisker_0._x[1]#Q3+1.5IQR\n",
    "        ycoord_0 = whisker_0._y[1]#Q3+1.5IQR\n",
    "        #Unhealthy\n",
    "        whisker_1 = lines[class_i*lines_unit*len(tempD2) + lines_unit*1 + 1]\n",
    "        xcoord_1 = whisker_1._x[1]#Q3+1.5IQR\n",
    "        ycoord_1 = whisker_1._y[1]#Q3+1.5IQR\n",
    "        #Standard point for annotation\n",
    "        xcoord = (xcoord_0+xcoord_1)/2\n",
    "        ycoord = max(ycoord_0, ycoord_1)\n",
    "        #Add annotation lines\n",
    "        aline_offset = yinter/5\n",
    "        aline_length = yinter/5 + aline_offset/2\n",
    "        ax.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",
    "        bmi_class = tempL1[class_i]\n",
    "        pval = tempDF3.loc[(bmi_class, bbmi), 'AdjPval_all']\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/12\n",
    "            text_size = 'medium'\n",
    "        else:\n",
    "            text_offset = yinter/3\n",
    "            text_size = 'x-small'\n",
    "        ax.annotate(label, xy=(xcoord, ycoord+text_offset),\n",
    "                    horizontalalignment='center', verticalalignment='bottom',\n",
    "                    fontsize=text_size, color='k')\n",
    "    #Facet settings\n",
    "    ax.set_title(bbmi, {'fontsize':'medium'})\n",
    "    xoff = 0.025\n",
    "    yoff = 0.01\n",
    "    rect = plt.Rectangle((xoff, 1+yoff), 1-xoff, 0.15,#Manual adjustment\n",
    "                         transform=ax.transAxes, facecolor=tempD1[bbmi], alpha=0.3,\n",
    "                         clip_on=False, linewidth=0, zorder=0.5)\n",
    "    ax.add_patch(rect)\n",
    "    #Change the default boxplot settings\n",
    "    for line in lines:\n",
    "        line.set_color('k')\n",
    "    for box in ax.artists:\n",
    "        box.set_edgecolor('k')\n",
    "    #Legend\n",
    "    if ax_i==len(tempD1)-1:\n",
    "        ax.legend(title='Metabolic\\ncondition', title_fontsize='medium', fontsize='medium',\n",
    "                  bbox_to_anchor=(1, 0.5), loc='center left', borderaxespad=0.5,\n",
    "                  handlelength=1.5, handletextpad=0.5)\n",
    "    else:\n",
    "        ax.get_legend().remove()\n",
    "##Save\n",
    "fileDir = './ExportFigures/'\n",
    "ipynbName = '220804_Multiomics-BMI-NatMed1stRevision_DeltaBMI-misclassification_ClinicalDefinition-ver2_'\n",
    "fileName = 'DeltaBMI-all.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 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-autonumbering": false,
  "toc-showtags": false,
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}