Code10_DeltaBMI-misclassification_obesity.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multiomics BMI Paper — Relationships between ∆BMI-based Misclassification and Obesity-related Features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***by Kengo Watanabe***  \n",
    "\n",
    "This Jupyter Notebook (with Python 3 kernel) assessed the relationships between ∆BMI-based misclassification (i.e., BMI class vs. biological BMI class) and obesity-related features (BMI-associated physiological measures, obesity-related health markers) using regression analysis (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 numeric physiological features (preprocessed): 210106_Biological-BMI-paper_data-cleaning_OLS-regression_numDF.tsv  \n",
    "* Arivale baseline associations between BMI and numeric physiological features: 220803_Multiomics-BMI-NatMed1stRevision_BMI-regression_NumericFeatures_result-summary.tsv  \n",
    "* Arivale baseline clinical labs (preprocessed): 210106_Biological-BMI-paper_data-cleaning_OLS-regression_chemDF.tsv  \n",
    "* OLS linear regressions of the retained variables on BMI: 220802_Multiomics-BMI-NatMed1stRevision_BMI-LASSO-bcoef_ChemBMI-BothSex-OLS.tsv  \n",
    "\n",
    "Output figures and tables:  \n",
    "* Figure 3d, 3e  \n",
    "* Tables for Supplementary Data 6  \n",
    "\n",
    "Original notebook (memo for my future tracing):  \n",
    "* dalek:\\[JupyterLab HOME\\]/220621_Multiomics-BMI-NatMedRevision/220804_Multiomics-BMI-NatMed1stRevision_DeltaBMI-misclassification_Features-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-derived misclassification 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",
    "#Misclassification\n",
    "tempL = ['MetBMI', 'ProtBMI', 'ChemBMI', 'CombiBMI']\n",
    "for bbmi in tempL:\n",
    "    tempL1 = []\n",
    "    for row_i in range(len(tempDF)):\n",
    "        bmi_class = tempDF['BMI_class'].iloc[row_i]\n",
    "        bbmi_class = tempDF[bbmi+'_class'].iloc[row_i]\n",
    "        if bmi_class==bbmi_class:\n",
    "            tempL1.append('Matched')\n",
    "        else:\n",
    "            tempL1.append('Mismatched')\n",
    "    tempDF['vs_'+bbmi+'_class'] = tempL1\n",
    "\n",
    "display(tempDF)\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\n",
    "for bmi_class in ['Normal', 'Obese']:\n",
    "    print(bmi_class)\n",
    "    tempDF1 = tempDF.loc[tempDF['BMI_class']==bmi_class]\n",
    "    for bbmi in tempL:\n",
    "        display(tempDF1.groupby('vs_'+bbmi+'_class')['BMI'].describe())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> –> It would be safer to adjust the baseline BMI in the statistical tests.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF = bmiDF\n",
    "\n",
    "#Select the misclassification 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",
    "tempDF1 = tempDF[tempL]\n",
    "\n",
    "display(tempDF1)\n",
    "\n",
    "#Update\n",
    "bmiDF = tempDF1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1-2. BMI-associated features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Import the preprocessed DF for baseline numeric features\n",
    "fileDir = '../210104_Biological-BMI-paper/ExportData/'\n",
    "ipynbName = '210106_Biological-BMI-paper_data-cleaning_OLS-regression_'\n",
    "fileName = 'numDF.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 = tempDF.rename(columns={'BaseAge':'Age'})\n",
    "tempDF = tempDF.drop(columns=['Race'])#Race has NaN in this cohort\n",
    "\n",
    "#Eliminate the exceptional feature that should have been eliminated during data cleaning\n",
    "feature = 'activities.calories'\n",
    "tempDF = tempDF.drop(columns=[feature])\n",
    "\n",
    "#Drop the covariates, which were separately prepared in this notebook\n",
    "tempDF = tempDF.loc[:, ~tempDF.columns.isin(bmiDF.columns.tolist())]\n",
    "\n",
    "display(tempDF)\n",
    "print('NaN in DF:', tempDF.isnull().to_numpy().sum(axis=None))\n",
    "\n",
    "numDF = tempDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Import the regression result summary\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220803_Multiomics-BMI-NatMed1stRevision_BMI-regression_NumericFeatures_'\n",
    "fileName = 'result-summary.tsv'\n",
    "tempDF = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t')\n",
    "tempDF = tempDF.set_index('FeatureID')\n",
    "\n",
    "#Retrieve the features that were significantly associated with BMI\n",
    "tempDF = tempDF.loc[tempDF['BMI_AdjPval_all']<0.05]\n",
    "print('Significantly associated features with BMI (FDR < 0.05):', len(tempDF))\n",
    "\n",
    "#Eliminate unnecessary info (just for the display in Jupyter notebook)\n",
    "tempL = ['MetBMI', 'ProtBMI', 'ChemBMI', 'CombiBMI']\n",
    "for bbmi in tempL:\n",
    "    tempDF = tempDF.loc[:, ~tempDF.columns.str.contains(bbmi)]\n",
    "\n",
    "display(tempDF)\n",
    "\n",
    "numDF_meta = tempDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Update data DF\n",
    "tempDF1 = numDF\n",
    "tempDF2 = numDF_meta\n",
    "tempDF1 = tempDF1[tempDF2.index.tolist()]\n",
    "display(tempDF1)\n",
    "\n",
    "numDF = tempDF1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1-3. Obesity-related features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> The definition of obesity-related health marker would depend on researchers/doctors. So, in this study, the following 10 features in the available Arivale datasets are representatively used:  \n",
    "> * HDL-cholesterol  \n",
    "> * LDL-cholesterol  \n",
    "> * Triglycerides  \n",
    "> * Glucose  \n",
    "> * Insulin  \n",
    "> * HOMA-IR  \n",
    "> –> Note that the above features were also included for the model of standard clinical measures.  \n",
    "> * Glycohemoglobin (HbA1c)  \n",
    "> * High-sensitivity CRP  \n",
    "> * Adiponectin  \n",
    "> * Vitamin D (25(OH)D)  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Import the preprocessed DF for baseline clinical lab tests\n",
    "fileDir = '../210104_Biological-BMI-paper/ExportData/'\n",
    "ipynbName = '210106_Biological-BMI-paper_data-cleaning_OLS-regression_'\n",
    "fileName = 'chemDF.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 = tempDF.rename(columns={'BaseAge':'Age'})\n",
    "tempDF = tempDF.drop(columns=['Race'])#Race has NaN in this cohort\n",
    "\n",
    "#Drop the covariates, which were separately prepared in this notebook\n",
    "tempDF = tempDF.loc[:, ~tempDF.columns.isin(bmiDF.columns.tolist())]\n",
    "\n",
    "display(tempDF)\n",
    "print('NaN in DF:', tempDF.isnull().to_numpy().sum(axis=None))\n",
    "\n",
    "chemDF = tempDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Prepare display labels for the target features\n",
    "tempD = {'HDL CHOL DIRECT':'HDL-cholesterol',\n",
    "         'LDL-CHOL CALCULATION':'LDL-cholesterol',\n",
    "         'TRIGLYCERIDES':'Triglycerides',\n",
    "         'GLUCOSE':'Glucose',\n",
    "         'INSULIN':'Insulin',\n",
    "         'HOMA-IR':'HOMA-IR',\n",
    "         'GLYCOHEMOGLOBIN A1C':'HbA1c',\n",
    "         'CRP HIGH SENSITIVITY':'High-sensitivity CRP',\n",
    "         'ADIPONECTIN, SERUM':'Adiponectin',\n",
    "         'VITAMIN D, 25-OH TOT':'Vitamin D (25(OH)D)'}\n",
    "tempS = pd.Series(tempD, name='FeatureLabel')\n",
    "tempS.index = tempS.index.rename('FeatureID')\n",
    "\n",
    "#Prepare R2 calculated in the chemBMI model analysis\n",
    "##-> Just for checking reference; in fact, not all the features were assessed.\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220802_Multiomics-BMI-NatMed1stRevision_BMI-LASSO-bcoef_'\n",
    "fileName = 'ChemBMI-BothSex-OLS.tsv'\n",
    "tempDF = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t')\n",
    "tempDF = tempDF.rename(columns={'Variable':'FeatureID'})\n",
    "tempDF = tempDF.set_index('FeatureID')\n",
    "\n",
    "#Merge\n",
    "tempDF = pd.merge(tempS, tempDF, left_index=True, right_index=True, how='left')\n",
    "\n",
    "display(tempDF)\n",
    "\n",
    "chemDF_meta = tempDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Update data DF\n",
    "tempDF1 = chemDF\n",
    "tempDF2 = chemDF_meta\n",
    "tempDF1 = tempDF1[tempDF2.index.tolist()]\n",
    "display(tempDF1)\n",
    "\n",
    "chemDF = tempDF1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Regression analysis for the BMI-associated features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Because the datasets are almost ready for use, only the simple processing steps are required for OLS linear regression:  \n",
    "> * Missingness: To maximize the sample size for each regression, dropping NaN is performed after selecting feature.  \n",
    "> * Centering: standardization is applied to both dependent and independent variables (including covariates).  \n",
    ">\n",
    "> Hence, the remaining processing steps are implemented during the for-loop for each regression. Of note, because the feature values were already preprocessed to reduce skewness through the data cleaning steps, 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: Feature ~ b0 + b1\\*C(Misclassification) + b2\\*BMI + b3\\*C(Sex) + b4\\*Age + b5\\*AncestryPCs  \n",
    "> Main aim: Assess the difference in each feature between the matched and mismatched BMI class.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF1 = numDF\n",
    "tempDF2 = bmiDF\n",
    "tempL1 = ['Normal', 'Obese']\n",
    "tempL2 = ['MetBMI', 'ProtBMI', 'ChemBMI', 'CombiBMI']\n",
    "tempDF3 = numDF_meta\n",
    "\n",
    "t_start = time.time()\n",
    "tempD1 = {}\n",
    "for feature in tempDF1.columns.tolist():\n",
    "    tempD2 = {}\n",
    "    for bmi_class in tempL1:\n",
    "        #Processing for OLS linear regression\n",
    "        ##Gather all necessary variables into a single DF\n",
    "        tempS = tempDF1[feature]\n",
    "        tempDF = pd.merge(tempS, tempDF2, 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 feature values\n",
    "        tempDF = tempDF.dropna()\n",
    "        ##Z-score transformation\n",
    "        tempDF4 = tempDF.select_dtypes(include=[np.number])\n",
    "        scaler = StandardScaler(copy=True, with_mean=True, with_std=True)\n",
    "        tempA = scaler.fit_transform(tempDF4)#Column direction\n",
    "        tempDF4 = pd.DataFrame(data=tempA, index=tempDF4.index, columns=tempDF4.columns)\n",
    "        ###Recover the categorical variables\n",
    "        tempDF5 = tempDF.select_dtypes(exclude=[np.number])\n",
    "        tempDF = pd.merge(tempDF4, tempDF5, left_index=True, right_index=True, how='left')\n",
    "        ##Rename the dependent variable\n",
    "        tempDF = tempDF.rename(columns={feature:'Feature'})\n",
    "        ##Add a constant for the intercept\n",
    "        ###–> In statsmodels, a constant is automatically added as well as R!\n",
    "        \n",
    "        tempD3 = {}\n",
    "        for bbmi in tempL2:\n",
    "            #Processing specific to each bBMI\n",
    "            ##Sort to make bcoef = 0 and 1 for Matched and Mismatched\n",
    "            tempDF = tempDF.sort_values(by='vs_'+bbmi+'_class', ascending=True)\n",
    "            ##One-hot encoding for categorical covariates\n",
    "            ###–> In statsmodels, categorical variables are automatically recognized!\n",
    "            \n",
    "            #OLS linear regression\n",
    "            ##Fit univariate model\n",
    "            formula = 'Feature ~ C(vs_'+bbmi+'_class)'\n",
    "            fit_res1 = smf.ols(formula, data=tempDF).fit()\n",
    "            ##Fit full model\n",
    "            formula = 'Feature ~ C(vs_'+bbmi+'_class)'\\\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",
    "            tempS = pd.Series().astype('float64')\n",
    "            ##Save the sample size for each group\n",
    "            size1 = len(tempDF.loc[tempDF['vs_'+bbmi+'_class']=='Matched'])\n",
    "            size2 = len(tempDF.loc[tempDF['vs_'+bbmi+'_class']=='Mismatched'])\n",
    "            tempS.loc['vs'+bbmi+'class_nMatched'] = size1\n",
    "            tempS.loc['vs'+bbmi+'class_nMismatched'] = size2\n",
    "            ##Save R2 [%]\n",
    "            tempS.loc['vs'+bbmi+'class_UnivarR2'] = fit_res1.rsquared*100\n",
    "            tempS.loc['vs'+bbmi+'class_R2'] = fit_res2.rsquared*100\n",
    "            ##Save beta-coefficient of the target variable\n",
    "            tempS.loc['vs'+bbmi+'class_Bcoef'] = fit_res2.params['C(vs_'+bbmi+'_class)[T.Mismatched]']\n",
    "            tempS.loc['vs'+bbmi+'class_BcoefSE'] = fit_res2.bse['C(vs_'+bbmi+'_class)[T.Mismatched]']\n",
    "            ##Save t-statistic of the target variable\n",
    "            tempS.loc['vs'+bbmi+'class_tStat'] = fit_res2.tvalues['C(vs_'+bbmi+'_class)[T.Mismatched]']\n",
    "            ##Save P-value of the target variable\n",
    "            tempS.loc['vs'+bbmi+'class_Pval'] = fit_res2.pvalues['C(vs_'+bbmi+'_class)[T.Mismatched]']\n",
    "            ##Add dummy adjusted P-value rows for now\n",
    "            tempS.loc['vs'+bbmi+'class_AdjPval_within1'] = 1.0\n",
    "            tempS.loc['vs'+bbmi+'class_AdjPval_within2'] = 1.0\n",
    "            tempS.loc['vs'+bbmi+'class_AdjPval_all'] = 1.0\n",
    "            \n",
    "            tempD3[bbmi] = tempS\n",
    "        \n",
    "        #Clean the results (pd.Series) across bBMIs\n",
    "        ##Prepare common summary metrics: sample size, residual degrees of freedom\n",
    "        tempS1 = pd.Series().astype('float64')\n",
    "        tempS1.loc['N'] = len(tempDF)\n",
    "        tempS1.loc['DoF'] = int(fit_res2.df_resid)#Use the last result object but same b/w bBMIs\n",
    "        ##Combine each result\n",
    "        tempS2 = pd.concat(list(tempD3.values()), axis=0)\n",
    "        tempS = pd.concat([tempS1, tempS2], axis=0)\n",
    "        ##Convert to DF while transposing\n",
    "        tempDF = pd.DataFrame(tempS.to_dict(), index=[0])\n",
    "        ##Clean DF\n",
    "        tempDF['N'] = tempDF['N'].astype('int64')\n",
    "        tempDF['DoF'] = tempDF['DoF'].astype('int64')\n",
    "        tempL = tempDF.loc[:, tempDF.columns.str.contains('_nM.*ed')].columns.tolist()\n",
    "        for col_n in tempL:\n",
    "            tempDF[col_n] = tempDF[col_n].astype('int64')\n",
    "        tempDF['BMIclass'] = bmi_class\n",
    "        \n",
    "        tempD2[bmi_class] = tempDF\n",
    "    \n",
    "    #Clean the results (pd.DataFrame) across BMI classes\n",
    "    tempDF = pd.concat(list(tempD2.values()), axis=0)\n",
    "    tempDF['FeatureID'] = feature\n",
    "    \n",
    "    #P-value adjustment (across BMI classes and bBMIs within the feature) by using Benjamini–Hochberg method\n",
    "    tempL = tempDF.loc[:, tempDF.columns.str.contains('_Pval')].columns.tolist()\n",
    "    tempDF4 = tempDF.reset_index().melt(var_name='bBMI', value_name='Pval', value_vars=tempL,\n",
    "                                        id_vars=['BMIclass', 'FeatureID'])\n",
    "    tempDF4['AdjPval'] = multi.multipletests(tempDF4['Pval'], alpha=0.05, method='fdr_bh',\n",
    "                                             is_sorted=False, returnsorted=False)[1]\n",
    "    tempDF4 = tempDF4.pivot(index=['BMIclass', 'FeatureID'], columns='bBMI', values='AdjPval')\n",
    "    tempDF4.columns = tempDF4.columns.str.replace('_Pval', '_AdjPval_within1')\n",
    "    ##Replace the dummy values with the adjusted p-values\n",
    "    tempL = [(bmi_class, feature) for bmi_class in tempL1]\n",
    "    tempDF4 = tempDF4.loc[tempL]#Sort just in case\n",
    "    tempL = tempDF.loc[:, tempDF.columns.str.contains('_AdjPval_within1')].columns.tolist()\n",
    "    for col_n in tempL:\n",
    "        tempDF[col_n] = tempDF4[col_n].tolist()\n",
    "    \n",
    "    tempD1[feature] = tempDF\n",
    "t_elapsed = time.time() - t_start\n",
    "print('Elapsed time for',\n",
    "      len(tempDF1.columns)*len(tempL1)*len(tempL2), 'OLS linear regressions (',\n",
    "      len(tempDF1.columns), 'features x',\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 features\n",
    "tempDF = pd.concat(list(tempD1.values()), axis=0)\n",
    "##Add the display label for visualization etc.\n",
    "tempS = tempDF3['FeatureLabel']\n",
    "tempDF = pd.merge(tempS, tempDF, on='FeatureID', how='right')\n",
    "##Clean the column order by setting index\n",
    "tempDF = tempDF.set_index(['BMIclass', 'FeatureID'])\n",
    "\n",
    "#P-value adjustment (across BMI classes and features within each bBMI) by using Benjamini–Hochberg method\n",
    "for bbmi in tempL2:\n",
    "    tempDF['vs'+bbmi+'class_AdjPval_within2'] = multi.multipletests(tempDF['vs'+bbmi+'class_Pval'],\n",
    "                                                                    alpha=0.05, method='fdr_bh',\n",
    "                                                                    is_sorted=False, returnsorted=False)[1]\n",
    "\n",
    "#P-value adjustment (across all tests) by using Benjamini–Hochberg method\n",
    "tempL = tempDF.loc[:, tempDF.columns.str.contains('_Pval')].columns.tolist()\n",
    "tempDF4 = tempDF.reset_index().melt(var_name='bBMI', value_name='Pval', value_vars=tempL,\n",
    "                                    id_vars=['BMIclass', 'FeatureID'])\n",
    "tempDF4['AdjPval'] = multi.multipletests(tempDF4['Pval'], alpha=0.05, method='fdr_bh',\n",
    "                                         is_sorted=False, returnsorted=False)[1]\n",
    "tempDF4 = tempDF4.pivot(index=['BMIclass', 'FeatureID'], columns='bBMI', values='AdjPval')\n",
    "tempDF4.columns = tempDF4.columns.str.replace('_Pval', '_AdjPval_all')\n",
    "##Replace the dummy values with the adjusted p-values\n",
    "tempL = tempDF.loc[:, tempDF.columns.str.contains('_AdjPval_all')].columns.tolist()\n",
    "for col_n in tempL:\n",
    "    tempDF[col_n] = tempDF4[col_n]\n",
    "\n",
    "tempDF = tempDF.sort_index(axis=0, ascending=True, key=lambda x:x.str.lower())\n",
    "display(tempDF)\n",
    "\n",
    "#Save\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220804_Multiomics-BMI-NatMed1stRevision_DeltaBMI-misclassification_Features-ver2_'\n",
    "fileName = 'result-summary-NumericFeatures.tsv'\n",
    "tempDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "\n",
    "resDF = tempDF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2-2. Significantly different features by misclassification"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF = resDF\n",
    "tempL1 = ['Normal', 'Obese']\n",
    "tempL2 = ['MetBMI', 'ProtBMI', 'ChemBMI', 'CombiBMI']\n",
    "\n",
    "#Significantly different features\n",
    "print('Significantly different features by misclassification (FDR < 0.05)')\n",
    "tempD = {}\n",
    "for bmi_class in tempL1:\n",
    "    tempDF1 = tempDF.loc[bmi_class]#MultiIndex\n",
    "    print(' - '+bmi_class)\n",
    "    tempL = []\n",
    "    for bbmi in tempL2:\n",
    "        #Extact significant features\n",
    "        tempDF2 = tempDF1.loc[tempDF1['vs'+bbmi+'class_AdjPval_all']<0.05]\n",
    "        tempL.append(tempDF2.index.tolist())\n",
    "        \n",
    "        #Cf. Adjustment across BMI classes and bBMIs within the feature\n",
    "        tempDF3 = tempDF1.loc[tempDF1['vs'+bbmi+'class_AdjPval_within1']<0.05]\n",
    "        #Cf. Adjustment across BMI classes and features within each bBMI\n",
    "        tempDF4 = tempDF1.loc[tempDF1['vs'+bbmi+'class_AdjPval_within2']<0.05]\n",
    "        #Cf. Nominal P-value\n",
    "        tempDF5 = tempDF1.loc[tempDF1['vs'+bbmi+'class_Pval']<0.05]\n",
    "        #tempL.append(tempDF5.index.tolist())\n",
    "        \n",
    "        print('   - '+bbmi+':', len(tempDF2),\n",
    "              '(cf. within-adjustment 1:', len(tempDF3),\n",
    "              ', within-adjustment 2:', len(tempDF4),\n",
    "              ', nominal P-value:', len(tempDF5), ')')\n",
    "    #Flatten and drop multiplicates\n",
    "    tempL = list(set(row_n for sublist in tempL for row_n in sublist))\n",
    "    tempDF1 = tempDF1.loc[tempL]\n",
    "    \n",
    "    print('   - Union:', len(tempL))\n",
    "    tempD[bmi_class] = tempL\n",
    "    \n",
    "    #Clean (just for the display in Jupyter notebook)\n",
    "    tempL3 = tempDF1.loc[:, tempDF1.columns.str.contains('Bcoef$')]\n",
    "    tempL4 = tempDF1.loc[:, tempDF1.columns.str.contains('Pval')]\n",
    "    tempL = [col_n for sublist in [tempL3, tempL4] for col_n in sublist]\n",
    "    tempDF1 = tempDF1[tempL]\n",
    "    tempDF1 = tempDF1.sort_values(by='vsMetBMIclass_Pval', ascending=True)\n",
    "    display(tempDF1)\n",
    "\n",
    "#Flatten and drop multiplicates\n",
    "tempL = list(set(row_n for sublist in tempD.values() for row_n in sublist))\n",
    "print(' - Union:', len(tempL))\n",
    "tempL = [(bmi_class, feature) for bmi_class in tempL1 for feature in tempL]\n",
    "tempDF1 = tempDF.loc[tempL]\n",
    "#Clean (just for the display in Jupyter notebook)\n",
    "tempL3 = tempDF1.loc[:, tempDF1.columns.str.contains('Bcoef$')]\n",
    "tempL4 = tempDF1.loc[:, tempDF1.columns.str.contains('Pval')]\n",
    "tempL = [col_n for sublist in [tempL3, tempL4] for col_n in sublist]\n",
    "tempDF1 = tempDF1[tempL]\n",
    "tempDF1 = tempDF1.sort_values(by='vsMetBMIclass_Pval', ascending=True)\n",
    "display(tempDF1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> –> P-value adjustment seems too conservative, though...  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2-3. Visualization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Because the result summary will be provided as supplementary table, only the representative features are visualized.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempD1 = {'waist-to-height_ratio':'Waist-to-height ratio',\n",
    "          'heartrate.resting':'Heart rate (resting)',\n",
    "          'activities.minutesVeryActive':'Physical activity (very active minutes)',\n",
    "          'polygenic_bmi_raw':'BMI polygenic risk score'}\n",
    "tempD2 = {'Underweight':'blue', 'Normal':'green', 'Overweight':'orange', 'Obese':'red'}\n",
    "tempD3 = {'MetBMI':'b', 'ProtBMI':'r', 'ChemBMI':'g', 'CombiBMI':'m'}\n",
    "tempD4 = {'waist-to-height_ratio':'tab:red',\n",
    "          'heartrate.resting':'tab:red',\n",
    "          'activities.minutesVeryActive':'tab:blue',\n",
    "          'polygenic_bmi_raw':'0.8'}\n",
    "tempL1 = ['Normal', 'Obese']\n",
    "tempDF1 = numDF\n",
    "tempDF2 = bmiDF\n",
    "tempDF3 = resDF\n",
    "\n",
    "for feature_i, feature in enumerate(tempD1.keys()):\n",
    "    #Prepare DF\n",
    "    tempS = tempDF1[feature]\n",
    "    tempDF = pd.merge(tempS, tempDF2, left_index=True, right_index=True, how='left')\n",
    "    ##Drop NaN in the feature values\n",
    "    tempDF = tempDF.dropna()\n",
    "    ##Z-score transformation (based on the whole distribution)\n",
    "    tempDF4 = tempDF.select_dtypes(include=[np.number])\n",
    "    scaler = StandardScaler(copy=True, with_mean=True, with_std=True)\n",
    "    tempA = scaler.fit_transform(tempDF4)#Column direction\n",
    "    tempDF4 = pd.DataFrame(data=tempA, index=tempDF4.index, columns=tempDF4.columns)\n",
    "    ###Recover the categorical variables\n",
    "    tempDF5 = tempDF.select_dtypes(exclude=[np.number])\n",
    "    tempDF = pd.merge(tempDF4, tempDF5, left_index=True, right_index=True, how='left')\n",
    "    \n",
    "    #Check sample size\n",
    "    print(tempD1[feature])\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",
    "        print('   - '+bmi_class+' BMI class')\n",
    "        for bbmi in tempD3.keys():\n",
    "            tempS = tempDF['vs_'+bbmi+'_class'].loc[tempDF['BMI_class']==bmi_class]\n",
    "            print('     - vs. '+bbmi+' class:',\n",
    "                  tempS.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=1+len(tempD3),\n",
    "                             figsize=(9, 3), sharex=False, sharey=True,\n",
    "                             gridspec_kw={'width_ratios':[1, 1, 1, 1, 1]})\n",
    "    axis_ymin = -3.2\n",
    "    axis_ymax = 3.8\n",
    "    ymin = -3\n",
    "    ymax = 3\n",
    "    yinter = 1\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",
    "    tempL = []#For legend handles and labels\n",
    "    for ax_i, ax in enumerate(axes.flat):\n",
    "        if ax_i==0:#Overall\n",
    "            sns.boxplot(data=tempDF, y=feature, x='BMI_class', order=tempD2.keys(),\n",
    "                        hue='BMI_class', hue_order=tempD2.keys(), dodge=False, palette=tempD2,\n",
    "                        showfliers=False,#flierprops={'marker':'o', 'markerfacecolor':'gray', 'alpha':0.4},\n",
    "                        showcaps=True, notch=True, ax=ax)\n",
    "            tempL.append(ax.get_legend_handles_labels())\n",
    "        else:#Misclassification\n",
    "            bbmi = list(tempD3.keys())[ax_i-1]\n",
    "            tempDF4 = tempDF.loc[tempDF['BMI_class'].isin(tempL1)]\n",
    "            tempD = {'Matched':'0.8', 'Mismatched':tempD3[bbmi]}\n",
    "            sns.boxplot(data=tempDF4, y=feature, x='BMI_class', order=tempL1,\n",
    "                        hue='vs_'+bbmi+'_class', hue_order=tempD.keys(), dodge=True, palette=tempD,\n",
    "                        showfliers=False,#flierprops={'marker':'o', 'markerfacecolor':'gray', 'alpha':0.4},\n",
    "                        showcaps=True, notch=True, ax=ax)\n",
    "            tempL.append(ax.get_legend_handles_labels())\n",
    "        #Axis settings\n",
    "        if ax_i==0:\n",
    "            plt.setp(ax, xlim=(0-margin, len(tempD2)-1+margin))#To eliminate excess white space\n",
    "            plt.setp(ax, xlabel='', ylabel='Tarnsformed value [a.u.]\\n('+r'$Z$'+'-score)')\n",
    "        else:\n",
    "            plt.setp(ax, xlim=(0-margin, len(tempD)-1+margin))#To eliminate excess white space\n",
    "            plt.setp(ax.get_yticklabels(), visible=False)\n",
    "            plt.setp(ax, xlabel='', ylabel='')\n",
    "        #ax.grid(axis='y', linestyle='--', color='black')\n",
    "        sns.despine()\n",
    "        plt.setp(ax.get_xticklabels(), rotation=70,\n",
    "                 horizontalalignment='right', verticalalignment='center', rotation_mode='anchor')\n",
    "        #Annotation\n",
    "        lines = ax.get_lines()#Line2D: [[Q1, Q1-1.5IQR], [Q3, Q3+1.5IQR], [Q1, Q1], [Q3, Q3], [Med, Med], [flier]]\n",
    "        if ax_i!=0:\n",
    "            #P-value annotation\n",
    "            lines_unit = 5 + int(False)#showfliers=False\n",
    "            for class_i in range(len(tempL1)):\n",
    "                #Matched\n",
    "                whisker_0 = lines[class_i*lines_unit*len(tempD) + 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",
    "                #Mismatched\n",
    "                whisker_1 = lines[class_i*lines_unit*len(tempD) + 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, feature), 'vs'+bbmi+'class_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",
    "        if ax_i==0:\n",
    "            ax.set_title('Overall', {'fontsize':'medium'})\n",
    "        else:\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=tempD3[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",
    "        #Remove the default legend\n",
    "        ax.get_legend().remove()\n",
    "        #Save ax position for figure title\n",
    "        if ax_i==0:\n",
    "            ax_pos_l = ax.get_position().bounds\n",
    "        elif ax_i==len(tempD3):\n",
    "            ax_pos_r = ax.get_position().bounds\n",
    "    #Add legend manually\n",
    "    for legend_i in range(len(tempL)):\n",
    "        h, l = tempL[legend_i]\n",
    "        if legend_i==0:\n",
    "            title_text = 'BMI class (overall)'\n",
    "            position = (0.0, -0.25)\n",
    "        else:\n",
    "            bbmi =list(tempD3.keys())[legend_i-1]\n",
    "            title_text = 'vs. '+bbmi+' class'\n",
    "            position = (0.3+0.3*((legend_i-1)//2), -0.25-0.35*((legend_i-1)%2))\n",
    "        legend = fig.legend(handles=h, labels=l, fontsize='medium',\n",
    "                            title=title_text, title_fontsize='medium',\n",
    "                            bbox_to_anchor=position, loc='upper left',\n",
    "                            labelspacing=0.25, handletextpad=0.5,\n",
    "                            borderaxespad=1, frameon=False)\n",
    "        plt.gca().add_artist(legend)\n",
    "    #Add figure title\n",
    "    xcoord_0 = ax_pos_l[0]\n",
    "    xcoord_1 = ax_pos_r[0]+ax_pos_r[2]\n",
    "    ycoord = ax_pos_l[1]+ax_pos_l[3]\n",
    "    yoff = 0.15\n",
    "    if tempD4[feature] in ['tab:red', 'tab:blue']:\n",
    "        text_color = 'white'\n",
    "    else:\n",
    "        text_color = 'black'\n",
    "    fig.suptitle(tempD1[feature], x=(xcoord_0+xcoord_1)/2, y=ycoord+yoff,\n",
    "                 fontsize='large', fontweight='bold', color=text_color,\n",
    "                 horizontalalignment='center', verticalalignment='bottom')\n",
    "    yoff = yoff - 0.015#Manual adjustment\n",
    "    rect = plt.Rectangle((xcoord_0, ycoord+yoff), xcoord_1-xcoord_0, 0.125,#Manual adjustment\n",
    "                         transform=fig.transFigure, facecolor=tempD4[feature],\n",
    "                         clip_on=False, linewidth=0, zorder=0)\n",
    "    fig.patches.extend([rect])\n",
    "    ##Save\n",
    "    fileDir = './ExportFigures/'\n",
    "    ipynbName = '220804_Multiomics-BMI-NatMed1stRevision_DeltaBMI-misclassification_Features-ver2_'\n",
    "    fileName = 'NumericFeature'+str(feature_i+1).zfill(2)+'.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()\n",
    "    print('')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Regression analysis for the obesity-related features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Because the datasets are almost ready for use, only the simple processing steps are required for OLS linear regression:  \n",
    "> * Missingness: To maximize the sample size for each regression, dropping NaN is performed after selecting feature.  \n",
    "> * Centering: standardization is applied to both dependent and independent variables (including covariates).  \n",
    ">\n",
    "> Hence, the remaining processing steps are implemented during the for-loop for each regression. Of note, because the feature values were already preprocessed to reduce skewness through the data cleaning steps, OLS linear regression (i.e., GLM with Gaussian family) can be used simply.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3-1. Perform OLS linear regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Model: Feature ~ b0 + b1\\*C(Misclassification) + b2\\*BMI + b3\\*C(Sex) + b4\\*Age + b5\\*AncestryPCs  \n",
    "> Main aim: Assess the difference in each feature between the matched and mismatched BMI class.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF1 = chemDF\n",
    "tempDF2 = bmiDF\n",
    "tempL1 = ['Normal', 'Obese']\n",
    "tempL2 = ['MetBMI', 'ProtBMI']\n",
    "tempDF3 = chemDF_meta\n",
    "\n",
    "t_start = time.time()\n",
    "tempD1 = {}\n",
    "for feature in tempDF1.columns.tolist():\n",
    "    tempD2 = {}\n",
    "    for bmi_class in tempL1:\n",
    "        #Processing for OLS linear regression\n",
    "        ##Gather all necessary variables into a single DF\n",
    "        tempS = tempDF1[feature]\n",
    "        tempDF = pd.merge(tempS, tempDF2, 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 feature values\n",
    "        tempDF = tempDF.dropna()\n",
    "        ##Z-score transformation\n",
    "        tempDF4 = tempDF.select_dtypes(include=[np.number])\n",
    "        scaler = StandardScaler(copy=True, with_mean=True, with_std=True)\n",
    "        tempA = scaler.fit_transform(tempDF4)#Column direction\n",
    "        tempDF4 = pd.DataFrame(data=tempA, index=tempDF4.index, columns=tempDF4.columns)\n",
    "        ###Recover the categorical variables\n",
    "        tempDF5 = tempDF.select_dtypes(exclude=[np.number])\n",
    "        tempDF = pd.merge(tempDF4, tempDF5, left_index=True, right_index=True, how='left')\n",
    "        ##Rename the dependent variable\n",
    "        tempDF = tempDF.rename(columns={feature:'Feature'})\n",
    "        ##Add a constant for the intercept\n",
    "        ###–> In statsmodels, a constant is automatically added as well as R!\n",
    "        \n",
    "        tempD3 = {}\n",
    "        for bbmi in tempL2:\n",
    "            #Processing specific to each bBMI\n",
    "            ##Sort to make bcoef = 0 and 1 for Matched and Mismatched\n",
    "            tempDF = tempDF.sort_values(by='vs_'+bbmi+'_class', ascending=True)\n",
    "            ##One-hot encoding for categorical covariates\n",
    "            ###–> In statsmodels, categorical variables are automatically recognized!\n",
    "            \n",
    "            #OLS linear regression\n",
    "            ##Fit univariate model\n",
    "            formula = 'Feature ~ C(vs_'+bbmi+'_class)'\n",
    "            fit_res1 = smf.ols(formula, data=tempDF).fit()\n",
    "            ##Fit full model\n",
    "            formula = 'Feature ~ C(vs_'+bbmi+'_class)'\\\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",
    "            tempS = pd.Series().astype('float64')\n",
    "            ##Save the sample size for each group\n",
    "            size1 = len(tempDF.loc[tempDF['vs_'+bbmi+'_class']=='Matched'])\n",
    "            size2 = len(tempDF.loc[tempDF['vs_'+bbmi+'_class']=='Mismatched'])\n",
    "            tempS.loc['vs'+bbmi+'class_nMatched'] = size1\n",
    "            tempS.loc['vs'+bbmi+'class_nMismatched'] = size2\n",
    "            ##Save R2 [%]\n",
    "            tempS.loc['vs'+bbmi+'class_UnivarR2'] = fit_res1.rsquared*100\n",
    "            tempS.loc['vs'+bbmi+'class_R2'] = fit_res2.rsquared*100\n",
    "            ##Save beta-coefficient of the target variable\n",
    "            tempS.loc['vs'+bbmi+'class_Bcoef'] = fit_res2.params['C(vs_'+bbmi+'_class)[T.Mismatched]']\n",
    "            tempS.loc['vs'+bbmi+'class_BcoefSE'] = fit_res2.bse['C(vs_'+bbmi+'_class)[T.Mismatched]']\n",
    "            ##Save t-statistic of the target variable\n",
    "            tempS.loc['vs'+bbmi+'class_tStat'] = fit_res2.tvalues['C(vs_'+bbmi+'_class)[T.Mismatched]']\n",
    "            ##Save P-value of the target variable\n",
    "            tempS.loc['vs'+bbmi+'class_Pval'] = fit_res2.pvalues['C(vs_'+bbmi+'_class)[T.Mismatched]']\n",
    "            ##Add dummy adjusted P-value rows for now\n",
    "            tempS.loc['vs'+bbmi+'class_AdjPval_within1'] = 1.0\n",
    "            tempS.loc['vs'+bbmi+'class_AdjPval_within2'] = 1.0\n",
    "            tempS.loc['vs'+bbmi+'class_AdjPval_all'] = 1.0\n",
    "            \n",
    "            tempD3[bbmi] = tempS\n",
    "        \n",
    "        #Clean the results (pd.Series) across bBMIs\n",
    "        ##Prepare common summary metrics: sample size, residual degrees of freedom\n",
    "        tempS1 = pd.Series().astype('float64')\n",
    "        tempS1.loc['N'] = len(tempDF)\n",
    "        tempS1.loc['DoF'] = int(fit_res2.df_resid)#Use the last result object but same b/w bBMIs\n",
    "        ##Combine each result\n",
    "        tempS2 = pd.concat(list(tempD3.values()), axis=0)\n",
    "        tempS = pd.concat([tempS1, tempS2], axis=0)\n",
    "        ##Convert to DF while transposing\n",
    "        tempDF = pd.DataFrame(tempS.to_dict(), index=[0])\n",
    "        ##Clean DF\n",
    "        tempDF['N'] = tempDF['N'].astype('int64')\n",
    "        tempDF['DoF'] = tempDF['DoF'].astype('int64')\n",
    "        tempL = tempDF.loc[:, tempDF.columns.str.contains('_nM.*ed')].columns.tolist()\n",
    "        for col_n in tempL:\n",
    "            tempDF[col_n] = tempDF[col_n].astype('int64')\n",
    "        tempDF['BMIclass'] = bmi_class\n",
    "        \n",
    "        tempD2[bmi_class] = tempDF\n",
    "    \n",
    "    #Clean the results (pd.DataFrame) across BMI classes\n",
    "    tempDF = pd.concat(list(tempD2.values()), axis=0)\n",
    "    tempDF['FeatureID'] = feature\n",
    "    \n",
    "    #P-value adjustment (across BMI classes and bBMIs within the feature) by using Benjamini–Hochberg method\n",
    "    tempL = tempDF.loc[:, tempDF.columns.str.contains('_Pval')].columns.tolist()\n",
    "    tempDF4 = tempDF.reset_index().melt(var_name='bBMI', value_name='Pval', value_vars=tempL,\n",
    "                                        id_vars=['BMIclass', 'FeatureID'])\n",
    "    tempDF4['AdjPval'] = multi.multipletests(tempDF4['Pval'], alpha=0.05, method='fdr_bh',\n",
    "                                             is_sorted=False, returnsorted=False)[1]\n",
    "    tempDF4 = tempDF4.pivot(index=['BMIclass', 'FeatureID'], columns='bBMI', values='AdjPval')\n",
    "    tempDF4.columns = tempDF4.columns.str.replace('_Pval', '_AdjPval_within1')\n",
    "    ##Replace the dummy values with the adjusted p-values\n",
    "    tempL = [(bmi_class, feature) for bmi_class in tempL1]\n",
    "    tempDF4 = tempDF4.loc[tempL]#Sort just in case\n",
    "    tempL = tempDF.loc[:, tempDF.columns.str.contains('_AdjPval_within1')].columns.tolist()\n",
    "    for col_n in tempL:\n",
    "        tempDF[col_n] = tempDF4[col_n].tolist()\n",
    "    \n",
    "    tempD1[feature] = tempDF\n",
    "t_elapsed = time.time() - t_start\n",
    "print('Elapsed time for',\n",
    "      len(tempDF1.columns)*len(tempL1)*len(tempL2), 'OLS linear regressions (',\n",
    "      len(tempDF1.columns), 'features x',\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 features\n",
    "tempDF = pd.concat(list(tempD1.values()), axis=0)\n",
    "##Add the display label for visualization etc.\n",
    "tempS = tempDF3['FeatureLabel']\n",
    "tempDF = pd.merge(tempS, tempDF, on='FeatureID', how='right')\n",
    "##Clean the column order by setting index\n",
    "tempDF = tempDF.set_index(['BMIclass', 'FeatureID'])\n",
    "\n",
    "#P-value adjustment (across BMI classes and features within each bBMI) by using Benjamini–Hochberg method\n",
    "for bbmi in tempL2:\n",
    "    tempDF['vs'+bbmi+'class_AdjPval_within2'] = multi.multipletests(tempDF['vs'+bbmi+'class_Pval'],\n",
    "                                                                    alpha=0.05, method='fdr_bh',\n",
    "                                                                    is_sorted=False, returnsorted=False)[1]\n",
    "\n",
    "#P-value adjustment (across all tests) by using Benjamini–Hochberg method\n",
    "tempL = tempDF.loc[:, tempDF.columns.str.contains('_Pval')].columns.tolist()\n",
    "tempDF4 = tempDF.reset_index().melt(var_name='bBMI', value_name='Pval', value_vars=tempL,\n",
    "                                    id_vars=['BMIclass', 'FeatureID'])\n",
    "tempDF4['AdjPval'] = multi.multipletests(tempDF4['Pval'], alpha=0.05, method='fdr_bh',\n",
    "                                         is_sorted=False, returnsorted=False)[1]\n",
    "tempDF4 = tempDF4.pivot(index=['BMIclass', 'FeatureID'], columns='bBMI', values='AdjPval')\n",
    "tempDF4.columns = tempDF4.columns.str.replace('_Pval', '_AdjPval_all')\n",
    "##Replace the dummy values with the adjusted p-values\n",
    "tempL = tempDF.loc[:, tempDF.columns.str.contains('_AdjPval_all')].columns.tolist()\n",
    "for col_n in tempL:\n",
    "    tempDF[col_n] = tempDF4[col_n]\n",
    "\n",
    "tempDF = tempDF.sort_index(axis=0, ascending=True, key=lambda x:x.str.lower())\n",
    "display(tempDF)\n",
    "\n",
    "#Save\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220804_Multiomics-BMI-NatMed1stRevision_DeltaBMI-misclassification_Features-ver2_'\n",
    "fileName = 'result-summary-ClinicalFeatures.tsv'\n",
    "tempDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "\n",
    "resDF = tempDF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3-2. Significantly different features by misclassification"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF = resDF\n",
    "tempL1 = ['Normal', 'Obese']\n",
    "tempL2 = ['MetBMI', 'ProtBMI']\n",
    "\n",
    "#Significantly different features\n",
    "print('Significantly different features by misclassification (FDR < 0.05)')\n",
    "tempD = {}\n",
    "for bmi_class in tempL1:\n",
    "    tempDF1 = tempDF.loc[bmi_class]#MultiIndex\n",
    "    print(' - '+bmi_class)\n",
    "    tempL = []\n",
    "    for bbmi in tempL2:\n",
    "        #Extact significant features\n",
    "        tempDF2 = tempDF1.loc[tempDF1['vs'+bbmi+'class_AdjPval_all']<0.05]\n",
    "        tempL.append(tempDF2.index.tolist())\n",
    "        \n",
    "        #Cf. Adjustment across BMI classes and bBMIs within the feature\n",
    "        tempDF3 = tempDF1.loc[tempDF1['vs'+bbmi+'class_AdjPval_within1']<0.05]\n",
    "        #Cf. Adjustment across BMI classes and features within each bBMI\n",
    "        tempDF4 = tempDF1.loc[tempDF1['vs'+bbmi+'class_AdjPval_within2']<0.05]\n",
    "        #Cf. Nominal P-value\n",
    "        tempDF5 = tempDF1.loc[tempDF1['vs'+bbmi+'class_Pval']<0.05]\n",
    "        #tempL.append(tempDF5.index.tolist())\n",
    "        \n",
    "        print('   - '+bbmi+':', len(tempDF2),\n",
    "              '(cf. within-adjustment 1:', len(tempDF3),\n",
    "              ', within-adjustment 2:', len(tempDF4),\n",
    "              ', nominal P-value:', len(tempDF5), ')')\n",
    "    #Flatten and drop multiplicates\n",
    "    tempL = list(set(row_n for sublist in tempL for row_n in sublist))\n",
    "    tempDF1 = tempDF1.loc[tempL]\n",
    "    \n",
    "    print('   - Union:', len(tempL))\n",
    "    tempD[bmi_class] = tempL\n",
    "    \n",
    "    #Clean (just for the display in Jupyter notebook)\n",
    "    tempL3 = tempDF1.loc[:, tempDF1.columns.str.contains('Bcoef$')]\n",
    "    tempL4 = tempDF1.loc[:, tempDF1.columns.str.contains('Pval')]\n",
    "    tempL = [col_n for sublist in [tempL3, tempL4] for col_n in sublist]\n",
    "    tempDF1 = tempDF1[tempL]\n",
    "    tempDF1 = tempDF1.sort_values(by='vsMetBMIclass_Pval', ascending=True)\n",
    "    display(tempDF1)\n",
    "\n",
    "#Flatten and drop multiplicates\n",
    "tempL = list(set(row_n for sublist in tempD.values() for row_n in sublist))\n",
    "print(' - Union:', len(tempL))\n",
    "tempL = [(bmi_class, feature) for bmi_class in tempL1 for feature in tempL]\n",
    "tempDF1 = tempDF.loc[tempL]\n",
    "#Clean (just for the display in Jupyter notebook)\n",
    "tempL3 = tempDF1.loc[:, tempDF1.columns.str.contains('Bcoef$')]\n",
    "tempL4 = tempDF1.loc[:, tempDF1.columns.str.contains('Pval')]\n",
    "tempL = [col_n for sublist in [tempL3, tempL4] for col_n in sublist]\n",
    "tempDF1 = tempDF1[tempL]\n",
    "tempDF1 = tempDF1.sort_values(by='vsMetBMIclass_Pval', ascending=True)\n",
    "display(tempDF1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> –> P-value adjustment seems too conservative, though...  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3-3. Visualization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Because the result summary will be provided as supplementary table, only the representative features are visualized.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempD1 = chemDF_meta['FeatureLabel'].to_dict()\n",
    "tempD2 = {'Underweight':'blue', 'Normal':'green', 'Overweight':'orange', 'Obese':'red'}\n",
    "tempD3 = {'MetBMI':'b', 'ProtBMI':'r'}\n",
    "tempD4 = {'HDL CHOL DIRECT':'tab:blue',\n",
    "          'LDL-CHOL CALCULATION':'tab:red',\n",
    "          'TRIGLYCERIDES':'tab:red',\n",
    "          'GLUCOSE':'tab:red',\n",
    "          'INSULIN':'tab:red',\n",
    "          'HOMA-IR':'tab:red',\n",
    "          'GLYCOHEMOGLOBIN A1C':'tab:red',\n",
    "          'CRP HIGH SENSITIVITY':'tab:red',\n",
    "          'ADIPONECTIN, SERUM':'tab:blue',\n",
    "          'VITAMIN D, 25-OH TOT':'0.8'}\n",
    "tempL1 = ['Normal', 'Obese']\n",
    "tempDF1 = chemDF\n",
    "tempDF2 = bmiDF\n",
    "tempDF3 = resDF\n",
    "\n",
    "for feature_i, feature in enumerate(tempD1.keys()):\n",
    "    #Prepare DF\n",
    "    tempS = tempDF1[feature]\n",
    "    tempDF = pd.merge(tempS, tempDF2, left_index=True, right_index=True, how='left')\n",
    "    ##Drop NaN in the feature values\n",
    "    tempDF = tempDF.dropna()\n",
    "    ##Z-score transformation (based on the whole distribution)\n",
    "    tempDF4 = tempDF.select_dtypes(include=[np.number])\n",
    "    scaler = StandardScaler(copy=True, with_mean=True, with_std=True)\n",
    "    tempA = scaler.fit_transform(tempDF4)#Column direction\n",
    "    tempDF4 = pd.DataFrame(data=tempA, index=tempDF4.index, columns=tempDF4.columns)\n",
    "    ###Recover the categorical variables\n",
    "    tempDF5 = tempDF.select_dtypes(exclude=[np.number])\n",
    "    tempDF = pd.merge(tempDF4, tempDF5, left_index=True, right_index=True, how='left')\n",
    "    \n",
    "    #Check sample size\n",
    "    print(tempD1[feature])\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",
    "        print('   - '+bmi_class+' BMI class')\n",
    "        for bbmi in tempD3.keys():\n",
    "            tempS = tempDF['vs_'+bbmi+'_class'].loc[tempDF['BMI_class']==bmi_class]\n",
    "            print('     - vs. '+bbmi+' class:',\n",
    "                  tempS.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=1+len(tempD3),\n",
    "                             figsize=(5.25, 3), sharex=False, sharey=True,\n",
    "                             gridspec_kw={'width_ratios':[1, 1, 1]})\n",
    "    axis_ymin = -3.2\n",
    "    axis_ymax = 4.2\n",
    "    ymin = -3\n",
    "    ymax = 3\n",
    "    yinter = 1\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",
    "    tempL = []#For legend handles and labels\n",
    "    for ax_i, ax in enumerate(axes.flat):\n",
    "        if ax_i==0:#Overall\n",
    "            sns.boxplot(data=tempDF, y=feature, x='BMI_class', order=tempD2.keys(),\n",
    "                        hue='BMI_class', hue_order=tempD2.keys(), dodge=False, palette=tempD2,\n",
    "                        showfliers=False,#flierprops={'marker':'o', 'markerfacecolor':'gray', 'alpha':0.4},\n",
    "                        showcaps=True, notch=True, ax=ax)\n",
    "            tempL.append(ax.get_legend_handles_labels())\n",
    "        else:#Misclassification\n",
    "            bbmi = list(tempD3.keys())[ax_i-1]\n",
    "            tempDF4 = tempDF.loc[tempDF['BMI_class'].isin(tempL1)]\n",
    "            tempD = {'Matched':'0.8', 'Mismatched':tempD3[bbmi]}\n",
    "            sns.boxplot(data=tempDF4, y=feature, x='BMI_class', order=tempL1,\n",
    "                        hue='vs_'+bbmi+'_class', hue_order=tempD.keys(), dodge=True, palette=tempD,\n",
    "                        showfliers=False,#flierprops={'marker':'o', 'markerfacecolor':'gray', 'alpha':0.4},\n",
    "                        showcaps=True, notch=True, ax=ax)\n",
    "            tempL.append(ax.get_legend_handles_labels())\n",
    "        #Axis settings\n",
    "        if ax_i==0:\n",
    "            plt.setp(ax, xlim=(0-margin, len(tempD2)-1+margin))#To eliminate excess white space\n",
    "            plt.setp(ax, xlabel='', ylabel='Tarnsformed value [a.u.]\\n('+r'$Z$'+'-score)')\n",
    "        else:\n",
    "            plt.setp(ax, xlim=(0-margin, len(tempD)-1+margin))#To eliminate excess white space\n",
    "            plt.setp(ax.get_yticklabels(), visible=False)\n",
    "            plt.setp(ax, xlabel='', ylabel='')\n",
    "        #ax.grid(axis='y', linestyle='--', color='black')\n",
    "        sns.despine()\n",
    "        plt.setp(ax.get_xticklabels(), rotation=70,\n",
    "                 horizontalalignment='right', verticalalignment='center', rotation_mode='anchor')\n",
    "        #Annotation\n",
    "        lines = ax.get_lines()#Line2D: [[Q1, Q1-1.5IQR], [Q3, Q3+1.5IQR], [Q1, Q1], [Q3, Q3], [Med, Med], [flier]]\n",
    "        if ax_i!=0:\n",
    "            #P-value annotation\n",
    "            lines_unit = 5 + int(False)#showfliers=False\n",
    "            for class_i in range(len(tempL1)):\n",
    "                #Matched\n",
    "                whisker_0 = lines[class_i*lines_unit*len(tempD) + 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",
    "                #Mismatched\n",
    "                whisker_1 = lines[class_i*lines_unit*len(tempD) + 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, feature), 'vs'+bbmi+'class_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",
    "        if ax_i==0:\n",
    "            ax.set_title('Overall', {'fontsize':'medium'})\n",
    "        else:\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=tempD3[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",
    "        #Remove the default legend\n",
    "        ax.get_legend().remove()\n",
    "        #Save ax position for figure title\n",
    "        if ax_i==0:\n",
    "            ax_pos_l = ax.get_position().bounds\n",
    "        elif ax_i==len(tempD3):\n",
    "            ax_pos_r = ax.get_position().bounds\n",
    "    #Add legend manually\n",
    "    for legend_i in range(len(tempL)):\n",
    "        h, l = tempL[legend_i]\n",
    "        if legend_i==0:\n",
    "            title_text = 'BMI class (overall)'\n",
    "            position = (-0.075, -0.25)\n",
    "        else:\n",
    "            bbmi =list(tempD3.keys())[legend_i-1]\n",
    "            title_text = 'vs. '+bbmi+' class'\n",
    "            position = (0.4+0.3*((legend_i-1)//2), -0.25-0.35*((legend_i-1)%2))\n",
    "        legend = fig.legend(handles=h, labels=l, fontsize='medium',\n",
    "                            title=title_text, title_fontsize='medium',\n",
    "                            bbox_to_anchor=position, loc='upper left',\n",
    "                            labelspacing=0.25, handletextpad=0.5,\n",
    "                            borderaxespad=1, frameon=False)\n",
    "        plt.gca().add_artist(legend)\n",
    "    #Add figure title\n",
    "    xcoord_0 = ax_pos_l[0]\n",
    "    xcoord_1 = ax_pos_r[0]+ax_pos_r[2]\n",
    "    ycoord = ax_pos_l[1]+ax_pos_l[3]\n",
    "    yoff = 0.15\n",
    "    if tempD4[feature] in ['tab:red', 'tab:blue']:\n",
    "        text_color = 'white'\n",
    "    else:\n",
    "        text_color = 'black'\n",
    "    fig.suptitle(tempD1[feature], x=(xcoord_0+xcoord_1)/2, y=ycoord+yoff,\n",
    "                 fontsize='large', fontweight='bold', color=text_color,\n",
    "                 horizontalalignment='center', verticalalignment='bottom')\n",
    "    yoff = yoff - 0.015#Manual adjustment\n",
    "    rect = plt.Rectangle((xcoord_0, ycoord+yoff), xcoord_1-xcoord_0, 0.125,#Manual adjustment\n",
    "                         transform=fig.transFigure, facecolor=tempD4[feature],\n",
    "                         clip_on=False, linewidth=0, zorder=0)\n",
    "    fig.patches.extend([rect])\n",
    "    ##Save\n",
    "    fileDir = './ExportFigures/'\n",
    "    ipynbName = '220804_Multiomics-BMI-NatMed1stRevision_DeltaBMI-misclassification_Features-ver2_'\n",
    "    fileName = 'ClinicalFeature'+str(feature_i+1).zfill(2)+'.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()\n",
    "    print('')"
   ]
  },
  {
   "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
}