Code12_DeltaBMI-misclassification_TwinsUK.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multiomics BMI Paper — ∆BMI-based Misclassification in the TwinsUK Cohort"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***by Kengo Watanabe***  \n",
    "\n",
    "This Jupyter Notebook (with Python 3 kernel) validated the findings about ∆BMI-based misclassification (i.e., BMI class vs. biological BMI class) with the TwinsUK cohort.  \n",
    "\n",
    "Input files:  \n",
    "* TwinsUK covariates: 220916_Multiomics-BMI-NatMed1stRevision_TwinsUK-DataCleaning-ver3_general-data_final.tsv  \n",
    "* TwinsUK MetBMI predictions: 220918_Multiomics-BMI-NatMed1stRevision_TwinsUK-BMI-LASSO-ver3_MetBMI-BothSex-TwinsUK.tsv  \n",
    "* TwinsUK obesity-related features (preprocessed): 220919_Multiomics-BMI-NatMed1stRevision_TwinsUK-Preprocessing-ver2_preprocessed-phenotype-dataset.tsv  \n",
    "* Arivale baseline MetBMI predictions (full panel version): 220803_Multiomics-BMI-NatMed1stRevision_DeltaBMI-misclassification_biologicalBMI-baseline-summary-BothSex.tsv  \n",
    "* Arivale baseline MetBMI predictions (restricted panel version): 220918_Multiomics-BMI-NatMed1stRevision_TwinsUK-BMI-LASSO-ver3_MetBMI-BothSex-Arivale.tsv  \n",
    "* TwinsUK metabolic health condition: 220919_Multiomics-BMI-NatMed1stRevision_TwinsUK-Preprocessing-ver2_metabolic-health-summary.tsv  \n",
    "* Arivale baseline metabolic health condition: 220720_Multiomics-BMI-NatMedRevision_Misclassification_metabolic-health-summary.tsv  \n",
    "\n",
    "Output figures and tables:  \n",
    "* Supplementary Figure 6  \n",
    "* Tables for Supplementary Data 6, 10  \n",
    "\n",
    "Original notebook (memo for my future tracing):  \n",
    "* dalek:\\[JupyterLab HOME\\]/220621_Multiomics-BMI-NatMedRevision/220920_Multiomics-BMI-NatMed1stRevision_TwinsUK-Misclassification-DeltaBMI-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 the cleaned dataframe\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220916_Multiomics-BMI-NatMed1stRevision_TwinsUK-DataCleaning-ver3_'\n",
    "fileName = 'general-data_final.tsv'\n",
    "tempDF = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t', parse_dates=['VisitDate'])\n",
    "tempDF = tempDF.set_index('KeyIndex')\n",
    "\n",
    "display(tempDF)\n",
    "print(' - Unique KeyIndex:', len(tempDF.index.unique()))\n",
    "tempDF1 = tempDF.index.to_series().str.split(pat='_', expand=True)\n",
    "print(' - Unique participant:', len(tempDF1.iloc[:, 0].unique()))\n",
    "\n",
    "bmiDF = tempDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF = bmiDF\n",
    "bbmi = 'MetBMI'\n",
    "\n",
    "#Add bBMI prediction\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220918_Multiomics-BMI-NatMed1stRevision_TwinsUK-BMI-LASSO-ver3_'\n",
    "fileName = bbmi+'-BothSex-TwinsUK.tsv'\n",
    "tempDF1 = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t')\n",
    "tempDF1 = tempDF1.set_index('KeyIndex')\n",
    "tempDF = pd.merge(tempDF, tempDF1[bbmi], left_index=True, right_index=True, how='left')\n",
    "\n",
    "#Add bBMI-based class\n",
    "tempL = []\n",
    "for value in tempDF[bbmi].tolist():\n",
    "    if np.isnan(value):\n",
    "        tempL.append('NotCalculated')\n",
    "    elif value < 18.5:\n",
    "        tempL.append('Underweight')\n",
    "    elif value < 25:\n",
    "        tempL.append('Normal')\n",
    "    elif value < 30:\n",
    "        tempL.append('Overweight')\n",
    "    elif value >= 30:\n",
    "        tempL.append('Obese')\n",
    "    else:#Just in case\n",
    "        tempL.append('Error?')\n",
    "tempDF[bbmi+'_class'] = tempL\n",
    "##Check\n",
    "tempL = []\n",
    "for bmi in ['BMI', bbmi]:\n",
    "    tempS1 = tempDF[bmi+'_class'].value_counts()\n",
    "    tempS2 = tempS1 / len(tempDF) * 100\n",
    "    tempS2.name = tempS2.name+' [%]'\n",
    "    tempDF1 = pd.merge(tempS1, tempS2, left_index=True, right_index=True, how='outer')\n",
    "    tempL.append(tempDF1)\n",
    "tempDF1 = pd.concat(tempL, axis=1)\n",
    "display(tempDF1)\n",
    "\n",
    "#Misclassification\n",
    "tempL = [bbmi]\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",
    "tempS = tempDF.groupby(by='BMI_class')['vs_'+bbmi+'_class'].value_counts()\n",
    "display(tempS)\n",
    "\n",
    "#Save\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220920_Multiomics-BMI-NatMed1stRevision_TwinsUK-Misclassification-DeltaBMI-ver2_'\n",
    "fileName = 'general-data-with-'+bbmi+'-BothSex-TwinsUK.tsv'\n",
    "tempDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "\n",
    "#Update\n",
    "bmiDF = tempDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF = bmiDF\n",
    "bbmi = 'MetBMI'\n",
    "\n",
    "#Check BMI\n",
    "for bmi_class in ['Normal', 'Obese']:\n",
    "    print(bmi_class)\n",
    "    tempDF1 = tempDF.loc[tempDF['BMI_class']==bmi_class]\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": "markdown",
   "metadata": {},
   "source": [
    "### 1-2. Obesity-related features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Based on the data availability, the following 12 features are 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)~~  \n",
    "> –> Note that all the above features were used as the obesity-related health markers in the Arivale analysis.  \n",
    "> * Systolic blood pressure  \n",
    "> * Diastolic blood pressure  \n",
    "> –> Note that these two features were used as the BMI-associated physiological measures in the Arivale analysis.  \n",
    "> * DEXA total fat percentage  \n",
    "> * DEXA android-to-gynoid ratio  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Import the preprocessed dataframe\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220919_Multiomics-BMI-NatMed1stRevision_TwinsUK-Preprocessing-ver2_'\n",
    "fileName = 'preprocessed-phenotype-dataset.tsv'\n",
    "tempDF = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t')\n",
    "tempDF = tempDF.set_index('KeyIndex')\n",
    "\n",
    "display(tempDF)\n",
    "print('NaN in DF:', tempDF.isnull().to_numpy().sum(axis=None))\n",
    "\n",
    "featureDF = tempDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF = featureDF\n",
    "\n",
    "#Prepare display labels for the target features\n",
    "tempD = {'HDL-cholesterol':'HDL-cholesterol',\n",
    "         'LDL-cholesterol':'LDL-cholesterol',\n",
    "         'Triglycerides':'Triglycerides',\n",
    "         'Glucose':'Glucose',\n",
    "         'Insulin':'Insulin',\n",
    "         'HOMA-IR':'HOMA-IR',\n",
    "         'Hs-CRP':'Hs-CRP',\n",
    "         'Adiponectin':'Adiponectin',\n",
    "         'SystolicBP':'Systolic BP',\n",
    "         'DiastolicBP':'Diastolic BP',\n",
    "         'DEXA-TotalFatPercentage':'Percent total fat',\n",
    "         'DEXA-AndroidGynoidRatio':'Android-to-gynoid'}\n",
    "tempS = pd.Series(tempD, name='FeatureLabel')\n",
    "tempS.index = tempS.index.rename('FeatureID')\n",
    "\n",
    "#Add value range (Not necessary but just for making DF)\n",
    "tempDF1 = tempDF.describe().T\n",
    "tempDF1.index = tempDF1.index.rename('FeatureID')\n",
    "\n",
    "#Merge\n",
    "tempDF1 = pd.merge(tempS, tempDF1, left_index=True, right_index=True, how='left')\n",
    "\n",
    "display(tempDF1)\n",
    "\n",
    "featureDF_meta = tempDF1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. 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": [
    "### 2-1. Perform OLS linear regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Model: Feature ~ b0 + b1\\*C(Misclassification) + b2\\*BMI + b3\\*C(Sex) + b4\\*Age  \n",
    "> Main aim: Assess the difference in each feature between the matched and mismatched BMI class.  \n",
    "\n",
    "> In this TwinsUK analysis, ancestry PCs are NOT included as the covariates due to data availability. Also, because only the MetBMI is available, the reference P-value adjustment 2 (across BMI classes and features within each bBMI) is skipped.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF1 = featureDF\n",
    "tempDF2 = bmiDF\n",
    "tempL1 = ['Normal', 'Obese']\n",
    "tempL2 = ['MetBMI']\n",
    "tempDF3 = featureDF_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'\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 = '220920_Multiomics-BMI-NatMed1stRevision_TwinsUK-Misclassification-DeltaBMI-ver2_'\n",
    "fileName = 'regression-summary_feature.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']\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": [
    "### 2-3. Visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempD1 = featureDF_meta['FeatureLabel'].to_dict()\n",
    "tempD2 = {'Underweight':'blue', 'Normal':'green', 'Overweight':'orange', 'Obese':'red'}\n",
    "tempD3 = {'MetBMI':'b'}\n",
    "tempD4 = {'HDL-cholesterol':'tab:blue',\n",
    "          'LDL-cholesterol':'tab:red',\n",
    "          'Triglycerides':'tab:red',\n",
    "          'Glucose':'0.8',\n",
    "          'Insulin':'0.8',\n",
    "          'HOMA-IR':'tab:red',\n",
    "          'Hs-CRP':'tab:red',\n",
    "          'Adiponectin':'0.8',\n",
    "          'SystolicBP':'0.8',\n",
    "          'DiastolicBP':'0.8',\n",
    "          'DEXA-TotalFatPercentage':'tab:red',\n",
    "          'DEXA-AndroidGynoidRatio':'tab:red'}\n",
    "tempL1 = ['Normal', 'Obese']\n",
    "tempDF1 = featureDF\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=(3.5, 3), sharex=False, sharey=True,\n",
    "                             gridspec_kw={'width_ratios':[1, 1]})\n",
    "    axis_ymin = -4.4\n",
    "    axis_ymax = 4.4\n",
    "    ymin = -4\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.95, 0.9)\n",
    "        else:\n",
    "            bbmi =list(tempD3.keys())[legend_i-1]\n",
    "            title_text = 'vs. '+bbmi+' class'\n",
    "            position = (0.95, 0.3)\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=0.0, 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 = '220920_Multiomics-BMI-NatMed1stRevision_TwinsUK-Misclassification-DeltaBMI-ver2_'\n",
    "    fileName = 'Feature'+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. ∆BMI-derived misclassification rate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3-1. Prepare the misclassification in Arivale"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> To directly compare between Arivale and TwinsUK, MetBMI calculated from the restricted version is prioritized for Arivale. At the same time, MetBMI calculated from the full version is also prepared as reference; in fact, the final figure has no P-value and thus doesn't require P-value adjustment.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bbmi = 'MetBMI'\n",
    "\n",
    "print('Restricted version:')\n",
    "#Import the Arivale bBMI predicted by restricted version model\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220918_Multiomics-BMI-NatMed1stRevision_TwinsUK-BMI-LASSO-ver3_'\n",
    "fileName = bbmi+'-BothSex-Arivale.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",
    "#Add BMI and bBMI-based class\n",
    "for bmi in ['BMI', bbmi]:\n",
    "    tempL = []\n",
    "    for value in tempDF[bmi].tolist():\n",
    "        if np.isnan(value):\n",
    "            tempL.append('NotCalculated')\n",
    "        elif value < 18.5:\n",
    "            tempL.append('Underweight')\n",
    "        elif value < 25:\n",
    "            tempL.append('Normal')\n",
    "        elif value < 30:\n",
    "            tempL.append('Overweight')\n",
    "        elif value >= 30:\n",
    "            tempL.append('Obese')\n",
    "        else:#Just in case\n",
    "            tempL.append('Error?')\n",
    "    tempDF[bmi+'_class'] = tempL\n",
    "##Check\n",
    "tempL = []\n",
    "for bmi in ['BMI', bbmi]:\n",
    "    tempS1 = tempDF[bmi+'_class'].value_counts()\n",
    "    tempS2 = tempS1 / len(tempDF) * 100\n",
    "    tempS2.name = tempS2.name+' [%]'\n",
    "    tempDF1 = pd.merge(tempS1, tempS2, left_index=True, right_index=True, how='outer')\n",
    "    tempL.append(tempDF1)\n",
    "tempDF1 = pd.concat(tempL, axis=1)\n",
    "display(tempDF1)\n",
    "\n",
    "#Misclassification\n",
    "tempL = [bbmi]\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",
    "tempS = tempDF.groupby(by='BMI_class')['vs_'+bbmi+'_class'].value_counts()\n",
    "display(tempS)\n",
    "\n",
    "bmiDF_a_restricted = tempDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bbmi = 'MetBMI'\n",
    "\n",
    "print('Full version:')\n",
    "#Import the Arivale bBMI predicted by full version model\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",
    "#Change column names\n",
    "tempDF.columns = tempDF.columns.str.replace('Base', '')\n",
    "\n",
    "#Check\n",
    "tempL = []\n",
    "for bmi in ['BMI', bbmi]:\n",
    "    tempS1 = tempDF[bmi+'_class'].value_counts()\n",
    "    tempS2 = tempS1 / len(tempDF) * 100\n",
    "    tempS2.name = tempS2.name+' [%]'\n",
    "    tempDF1 = pd.merge(tempS1, tempS2, left_index=True, right_index=True, how='outer')\n",
    "    tempL.append(tempDF1)\n",
    "tempDF1 = pd.concat(tempL, axis=1)\n",
    "display(tempDF1)\n",
    "\n",
    "#Misclassification\n",
    "tempL = [bbmi]\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",
    "tempS = tempDF.groupby(by='BMI_class')['vs_'+bbmi+'_class'].value_counts()\n",
    "display(tempS)\n",
    "\n",
    "bmiDF_a_full = tempDF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3-2. Misclassification rate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempD1 = {'Arivale (full)':bmiDF_a_full,\n",
    "          'Arivale (restricted)':bmiDF_a_restricted,\n",
    "          'TwinsUK':bmiDF}\n",
    "tempL1 = ['MetBMI']\n",
    "tempL2 = ['Underweight', 'Normal', 'Overweight', 'Obese']\n",
    "tempD2 = {'Arivale (full): vs. MetBMI class':'Arivale (full)',\n",
    "          'Arivale (restricted): vs. MetBMI class':'Arivale (restricted)',\n",
    "          'TwinsUK: vs. MetBMI class':'TwinsUK'}\n",
    "tempD3 = {'Arivale (full)':plt.get_cmap('tab20')(1),\n",
    "          'Arivale (restricted)':'tab:blue',\n",
    "          'TwinsUK':'tab:orange'}\n",
    "legend_title = 'vs. MetBMI class'\n",
    "\n",
    "#Calculate misclassification rate based on each biological BMI class\n",
    "tempD = {}\n",
    "for cohort in tempD1.keys():\n",
    "    for bbmi in tempL1:\n",
    "        tempDF = tempD1[cohort]\n",
    "        \n",
    "        #Just in case\n",
    "        tempDF1 = tempDF.loc[tempDF[bbmi+'_class']=='NotCalculated']\n",
    "        if len(tempDF1)>0:\n",
    "            print('Check NotCalculated in '+cohort+' '+bbmi+' class.')\n",
    "            tempDF = tempDF.loc[tempDF['Base'+bbmi+'_class']!='NotCalculated']\n",
    "        \n",
    "        #Count misclassification\n",
    "        counter0 = 0\n",
    "        counter1 = 0\n",
    "        counter2 = 0\n",
    "        counter3 = 0\n",
    "        counter4 = 0\n",
    "        for row_n in tempDF.index.tolist():\n",
    "            bmi_class = tempDF.loc[row_n, 'BMI_class']\n",
    "            bbmi_class = tempDF.loc[row_n, bbmi+'_class']\n",
    "            if bmi_class!=bbmi_class:\n",
    "                counter0 += 1#Overall\n",
    "                if bmi_class=='Underweight':\n",
    "                    counter1 += 1\n",
    "                elif bmi_class=='Normal':\n",
    "                    counter2 += 1\n",
    "                elif bmi_class=='Overweight':\n",
    "                    counter3 += 1\n",
    "                elif bmi_class=='Obese':\n",
    "                    counter4 += 1\n",
    "                else:#Just in case\n",
    "                    print('Error?')\n",
    "        \n",
    "        #Clean the count result\n",
    "        tempS = pd.Series(name=cohort+': vs. '+bbmi+' class')\n",
    "        tempS.loc['Overall count'] = counter0\n",
    "        tempS.loc['Overall [%]'] = counter0/len(tempDF)*100\n",
    "        tempL = [counter1, counter2, counter3, counter4]\n",
    "        for counter_i in range(len(tempL)):\n",
    "            bmi_class = tempL2[counter_i]#Because of the same order\n",
    "            counter = tempL[counter_i]\n",
    "            tempS.loc[bmi_class+' count'] = counter\n",
    "            total = len(tempDF.loc[tempDF['BMI_class']==bmi_class])\n",
    "            tempS.loc[bmi_class+' [%]'] = counter / total * 100\n",
    "        tempD[cohort+': vs. '+bbmi+' class'] = tempS\n",
    "tempDF = pd.concat(list(tempD.values()), axis=1)\n",
    "display(tempDF)\n",
    "\n",
    "#Plot (without underweight)\n",
    "##Prepare DF\n",
    "tempDF = tempDF.loc[['Overall [%]', 'Normal [%]', 'Overweight [%]', 'Obese [%]']]\n",
    "tempDF.index = tempDF.index.str.replace(' \\[%\\]', '')\n",
    "tempDF.columns = tempDF.columns.map(tempD2)\n",
    "tempDF = tempDF.reset_index().melt(var_name='Category', value_name='Misclassification', id_vars='index')\n",
    "##Visualization\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "plt.figure(figsize=(3, 3))\n",
    "sns.barplot(data=tempDF, x='index', y='Misclassification',\n",
    "            hue='Category', dodge=True, palette=tempD3, edgecolor='black')\n",
    "sns.despine()\n",
    "plt.yticks(np.arange(0, 61, 10))\n",
    "plt.axvline(x=(0+1)/2, **{'linestyle':'--', 'color':'k'})\n",
    "##Add reference range\n",
    "plt.axhspan(ymin=28, ymax=48, facecolor='crimson', alpha=0.2, zorder=0)\n",
    "plt.axhline(y=28, **{'linestyle':'-', 'color':'crimson', 'zorder':0})\n",
    "plt.axhline(y=48, **{'linestyle':'-', 'color':'crimson', 'zorder':0})\n",
    "plt.ylabel('Misclassification [%]')\n",
    "#plt.xlabel('BMI class')\n",
    "plt.xlabel('')\n",
    "plt.xticks(rotation=70,\n",
    "           horizontalalignment='right', verticalalignment='center', rotation_mode='anchor')\n",
    "plt.legend(title=legend_title, bbox_to_anchor=(1, -0.6), loc='upper right',\n",
    "           handletextpad=0.5, borderaxespad=0)\n",
    "##Save\n",
    "fileDir = './ExportFigures/'\n",
    "ipynbName = '220920_Multiomics-BMI-NatMed1stRevision_TwinsUK-Misclassification-DeltaBMI-ver2_'\n",
    "fileName = 'misclassification-rate.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": [
    "## 4. ∆BMI and clinical definition-based misclassification"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4-1. Calculate ∆BMI and prepare covariates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempD1 = {'Arivale':bmiDF_a_restricted,#Restricted version of MetBMI\n",
    "          'TwinsUK':bmiDF}\n",
    "tempL1 = ['MetBMI']\n",
    "\n",
    "tempD2 = {}\n",
    "for cohort in tempD1.keys():\n",
    "    tempDF = tempD1[cohort]\n",
    "    #Calculate the rate of difference\n",
    "    for bbmi in tempL1:\n",
    "        tempDF['Delta'+bbmi] = (tempDF[bbmi] - tempDF['BMI']) / tempDF['BMI'] * 100\n",
    "    tempD2[cohort] = tempDF\n",
    "    \n",
    "    #Check skewness\n",
    "    tempDF1 = tempDF.select_dtypes(include=[np.number])\n",
    "    tempDF2 = tempDF1.describe()\n",
    "    tempDF2.loc['Skewness'] = stats.skew(tempDF1)\n",
    "    print(cohort)\n",
    "    display(tempDF2)\n",
    "#Update\n",
    "bmiDF_a_restricted = tempD2['Arivale']\n",
    "bmiDF = tempD2['TwinsUK']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> –> ∆BMI can be surly assumed as normal distribution.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF1 = bmiDF_a_restricted#Restricted version of MetBMI\n",
    "tempDF2 = bmiDF_a_full\n",
    "tempL = ['Sex', 'Age', 'PC1', 'PC2', 'PC3', 'PC4', 'PC5']\n",
    "\n",
    "#Add covariates\n",
    "tempDF = pd.merge(tempDF1, tempDF2[tempL], left_index=True, right_index=True, how='left')\n",
    "display(tempDF.describe(include='all'))\n",
    "print('NaN in DF:', tempDF.isnull().to_numpy().sum(axis=None))\n",
    "\n",
    "#Update\n",
    "bmiDF_a_restricted = tempDF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4-2. Metabolic health condiiton"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempD1 = {'Arivale':'220720_Multiomics-BMI-NatMedRevision_Misclassification_',\n",
    "          'TwinsUK':'220919_Multiomics-BMI-NatMed1stRevision_TwinsUK-Preprocessing-ver2_'}\n",
    "tempD2 = {'Arivale':'public_client_id',\n",
    "          'TwinsUK':'KeyIndex'}\n",
    "\n",
    "tempD3 = {}\n",
    "for cohort in tempD1.keys():\n",
    "    #Import cleaned table for metabolic health condition\n",
    "    fileDir = './ExportData/'\n",
    "    ipynbName = tempD1[cohort]\n",
    "    fileName = 'metabolic-health-summary.tsv'\n",
    "    tempDF = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t', dtype={tempD2[cohort]: str})\n",
    "    tempDF = tempDF.set_index(tempD2[cohort])\n",
    "    tempD3[cohort] = tempDF\n",
    "    \n",
    "    print(cohort)\n",
    "    display(tempDF)\n",
    "    display(tempDF.describe(include='all'))\n",
    "    print('NaN in DF:', tempDF.isnull().to_numpy().sum(axis=None))\n",
    "    print('')\n",
    "metabDF_a = tempD3['Arivale']\n",
    "metabDF_t = tempD3['TwinsUK']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> –> 6 participants in Arivale and 1,203 participants in TwinsUK had NaN for the metabolic health condition.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempD1 = {'Arivale':bmiDF_a_restricted,#Restricted version of MetBMI\n",
    "          'TwinsUK':bmiDF}\n",
    "tempD2 = {'Arivale':metabDF_a,\n",
    "          'TwinsUK':metabDF_t}\n",
    "\n",
    "#Check BMI between metabolically healthy and unhealthy groups\n",
    "for cohort in tempD1.keys():\n",
    "    print(cohort)\n",
    "    tempDF1 = tempD1[cohort]\n",
    "    tempDF2 = tempD2[cohort]\n",
    "    tempDF = pd.merge(tempDF1, tempDF2['Metabolically'],\n",
    "                      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())\n",
    "    print('')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> –> It would be safer to adjust the baseline BMI in the statistical tests.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4-3. Perform OLS linear regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Because ∆BMI values are assumed as normal distribution, OLS linear regression (i.e., GLM with Gaussian family) can be used simply.  \n",
    "\n",
    "> 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.  \n",
    ">\n",
    "> In this TwinsUK analysis, ancestry PCs are NOT included as the covariates due to data availability. Also, because only the MetBMI is assessed, the reference P-value adjustment (across bBMI classes within BMI class and within cohort) is skipped.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempD1 = {'Arivale':bmiDF_a_restricted,#Restricted version of MetBMI\n",
    "          'TwinsUK':bmiDF}\n",
    "tempD2 = {'Arivale':metabDF_a,\n",
    "          'TwinsUK':metabDF_t}\n",
    "tempD3 = {'Arivale':'+ BMI + C(Sex) + Age + PC1 + PC2 + PC3 + PC4 + PC5',\n",
    "          'TwinsUK':'+ BMI + C(Sex) + Age'}\n",
    "tempL1 = ['Normal', 'Obese']\n",
    "tempL2 = ['MetBMI']\n",
    "\n",
    "tempD4 = {}\n",
    "for cohort in tempD1.keys():\n",
    "    tempDF1 = tempD1[cohort]\n",
    "    tempDF2 = tempD2[cohort]\n",
    "    \n",
    "    t_start = time.time()\n",
    "    tempD5 = {}\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",
    "        tempD6 = {}\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)'+tempD3[cohort]\n",
    "            fit_res2 = smf.ols(formula, data=tempDF).fit()\n",
    "            \n",
    "            #Summarize the result\n",
    "            if cohort=='Arivale':\n",
    "                tempDF3 = pd.DataFrame({'DeltaBMI':[bbmi+' (restricted)']})\n",
    "            else:\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",
    "            tempD6[bbmi] = tempDF3\n",
    "        \n",
    "        #Clean the results (pd.DataFrame) across bBMIs\n",
    "        tempDF = pd.concat(list(tempD6.values()), axis=0)\n",
    "        tempDF['BMIclass'] = bmi_class\n",
    "        tempD5[bmi_class] = tempDF\n",
    "    t_elapsed = time.time() - t_start\n",
    "    print(cohort)\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(tempD5.values()), axis=0)\n",
    "    tempDF['Cohort'] = cohort\n",
    "    ##P-value adjustment (across BMI classes within cohort) by using Benjamini–Hochberg method\n",
    "    tempDF['AdjPval_within'] = multi.multipletests(tempDF['Pval'], alpha=0.05, method='fdr_bh',\n",
    "                                                   is_sorted=False, returnsorted=False)[1]\n",
    "    tempD4[cohort] = tempDF\n",
    "#Clean the results (pd.DataFrame) across cohorts\n",
    "tempDF = pd.concat(list(tempD4.values()), axis=0)\n",
    "##Clean the column order by setting index\n",
    "tempDF = tempDF.set_index(['Cohort', 'BMIclass', 'DeltaBMI'])\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 = '220920_Multiomics-BMI-NatMed1stRevision_TwinsUK-Misclassification-DeltaBMI-ver2_'\n",
    "fileName = 'regression-summary_DeltaBMI.tsv'\n",
    "tempDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "\n",
    "resDF = tempDF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4-4. Visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bbmi = 'MetBMI'\n",
    "tempD1 = {'Arivale':'tab:blue', 'TwinsUK':'tab:orange'}\n",
    "tempD2 = {'Healthy':'0.8', 'Unhealthy':'crimson'}\n",
    "tempL1 = ['Normal', 'Obese']\n",
    "tempD3 = {'Arivale':bmiDF_a_restricted,#Restricted version of MetBMI\n",
    "          'TwinsUK':bmiDF}\n",
    "tempD4 = {'Arivale':metabDF_a,\n",
    "          'TwinsUK':metabDF_t}\n",
    "tempDF3 = resDF\n",
    "\n",
    "#Visualization\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "fig, axes = plt.subplots(nrows=1, ncols=len(tempD1),\n",
    "                         figsize=(3.5, 3), sharex=True, sharey=True,\n",
    "                         gridspec_kw={'width_ratios':[1, 1]})\n",
    "axis_ymin = -35\n",
    "axis_ymax = 55\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",
    "    cohort = list(tempD1.keys())[ax_i]\n",
    "    tempDF1 = tempD3[cohort]\n",
    "    tempDF2 = tempD4[cohort]\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",
    "    #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",
    "    #Plot\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$'+bbmi+' [% 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[(cohort, bmi_class), 'AdjPval_all'].iloc[0]#3-level MultiIndex\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(cohort, {'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[cohort], 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 condition', title_fontsize='medium', fontsize='medium',\n",
    "                  bbox_to_anchor=(1, -0.45), loc='upper right', borderaxespad=0, ncol=2,\n",
    "                  handlelength=1.5, handletextpad=0.5, columnspacing=1.0)\n",
    "    else:\n",
    "        ax.get_legend().remove()\n",
    "##Save\n",
    "fileDir = './ExportFigures/'\n",
    "ipynbName = '220920_Multiomics-BMI-NatMed1stRevision_TwinsUK-Misclassification-DeltaBMI-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 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
}