Code07_BMI-LASSO-vs-others.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multiomics BMI Paper — Comparison of Biological BMI Models between LASSO and the Other Methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***by Kengo Watanabe***  \n",
    "\n",
    "This Jupyter Notebook (with Python 3 kernel) compared the blood omics-based BMI models (model performance, BMI predictions, and predictor variables) between LASSO and the other methods (elastic net, ridge, and random forest).  \n",
    "\n",
    "Input files:  \n",
    "* Arivale baseline BMI predictions with LASSO models: 220801_Multiomics-BMI-NatMed1stRevision_BMI-baseline-LASSO_\\[MetBMI/ProtBMI/ChemBMI/CombiBMI\\]-\\[Female/Male/BothSex\\].tsv  \n",
    "* Arivale baseline BMI predictions with elastic net models: 220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-ElasticNet_\\[MetBMI/ProtBMI/ChemBMI/CombiBMI\\]-\\[Female/Male/BothSex\\].tsv  \n",
    "* Arivale baseline BMI predictions with ridge models: 220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-Ridge_\\[MetBMI/ProtBMI/ChemBMI/CombiBMI\\]-\\[Female/Male/BothSex\\].tsv  \n",
    "* Arivale baseline BMI predictions with random forest models: 220828_Multiomics-BMI-NatMed1stRevision_BMI-baseline-RF-wenceslaus_\\[MetBMI/ProtBMI/ChemBMI/CombiBMI\\]-\\[Female/Male/BothSex\\].tsv  \n",
    "* LASSO models of biological BMI: 220801_Multiomics-BMI-NatMed1stRevision_BMI-baseline-LASSO_\\[MetBMI/ProtBMI/ChemBMI/CombiBMI\\]-BothSex-LASSObcoefs.tsv  \n",
    "* Elastic net models of biological BMI: 220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-ElasticNet_\\[MetBMI/ProtBMI/ChemBMI/CombiBMI\\]-BothSex-ElasticNetbcoefs.tsv  \n",
    "* Ridge models of biological BMI: 220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-Ridge_\\[MetBMI/ProtBMI/ChemBMI/CombiBMI\\]-BothSex-Ridgebcoefs.tsv  \n",
    "* Random forest models of biological BMI: 220828_Multiomics-BMI-NatMed1stRevision_BMI-baseline-RF-wenceslaus_\\[MetBMI/ProtBMI/ChemBMI/CombiBMI\\]-BothSex-feature-importance.tsv  \n",
    "\n",
    "Output figures and tables:  \n",
    "* Supplementary Figure 4  \n",
    "* Tables for Supplementary Data 10  \n",
    "\n",
    "Original notebook (memo for my future tracing):  \n",
    "* dalek:\\[JupyterLab HOME\\]/220621_Multiomics-BMI-NatMedRevision/220901_Multiomics-BMI-NatMed1stRevision_bBMI-LASSO-vs-others.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.metrics import r2_score\n",
    "from statsmodels.stats import weightstats\n",
    "from statsmodels.stats import multitest as multi\n",
    "from decimal import Decimal, ROUND_HALF_UP\n",
    "import sys\n",
    "#!pip install venn\n",
    "from venn import venn\n",
    "import matplotlib.patches as mpatches\n",
    "\n",
    "!conda list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Prepare predictions and reproduce out-of-sample R2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1-1. LASSO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#LASSO\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220801_Multiomics-BMI-NatMed1stRevision_BMI-baseline-LASSO_'\n",
    "yvar = 'BMI'\n",
    "tempL1 = ['Met'+yvar, 'Prot'+yvar, 'Chem'+yvar, 'Combi'+yvar]\n",
    "tempL2 = ['Testing', 'log_Base'+yvar, 'Base'+yvar]\n",
    "tempD1 = {}\n",
    "tempD2 = {}\n",
    "for sex in ['Female', 'Male', 'BothSex']:\n",
    "    #Clean all mesured and predicted values into a dataframe\n",
    "    for yvar_model in tempL1:\n",
    "        #Import the predictions dataframe\n",
    "        fileName = yvar_model+'-'+sex+'.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",
    "        #Put aside the common part\n",
    "        if yvar_model==tempL1[0]:\n",
    "            tempDF1 = tempDF[tempL2]#Initialize\n",
    "        tempDF = tempDF.drop(columns=tempL2)\n",
    "        #Merge\n",
    "        tempDF1 = pd.merge(tempDF1, tempDF, left_index=True, right_index=True, how='inner')\n",
    "    ##Change column names\n",
    "    tempDF1.columns = tempDF1.columns.str.replace('Base', '')\n",
    "    \n",
    "    #Reproduce out-of-sample R2\n",
    "    tempL = tempDF1['Testing'].unique().tolist()\n",
    "    tempDF2 = pd.DataFrame(index=pd.Index(tempL, name='Model'))\n",
    "    for yvar_model in tempL1:\n",
    "        tempL = []\n",
    "        for model_n in tempDF2.index.tolist():\n",
    "            tempS1 = tempDF1['log_'+yvar].loc[tempDF1['Testing']==model_n]\n",
    "            tempS2 = tempDF1['log_'+yvar_model].loc[tempDF1['Testing']==model_n]\n",
    "            tempL.append(r2_score(tempS1, tempS2, sample_weight=None))\n",
    "        tempDF2[yvar_model] = tempL\n",
    "    \n",
    "    tempD1[sex] = tempDF1\n",
    "    tempD2[sex] = tempDF2\n",
    "    \n",
    "    #Check\n",
    "    print(sex)\n",
    "    ##Reproduce the previous outputs (because R2 values are rounded in the following summary table)\n",
    "    for yvar_model in tempL1:\n",
    "        tempL = tempDF2[yvar_model].tolist()\n",
    "        print(' - '+yvar_model)\n",
    "        print('   - Out-of-sample R2 [Mean ± SD]:',\n",
    "              np.mean(tempL), '±', np.std(tempL, ddof=1))#Sample standard deviation\n",
    "        print('   - Out-of-sample R2 [Mean ± SEM]:',\n",
    "              np.mean(tempL), '±', np.std(tempL, ddof=1)/np.sqrt(len(tempL)))\n",
    "        print('   - Observed vs. predicted log_'+yvar+': (Pearson\\'s r, P) =',\n",
    "              stats.pearsonr(tempDF1['log_'+yvar], tempDF1['log_'+yvar_model]))\n",
    "        print('   - Observed vs. predicted '+yvar+': (Pearson\\'s r, P) =',\n",
    "              stats.pearsonr(tempDF1[yvar], tempDF1[yvar_model]))\n",
    "    print(' - Summary of out-of-sample R2 values')\n",
    "    display(tempDF2.describe(include='all'))\n",
    "    print(' - Summary of all measured and predicted values')\n",
    "    display(tempDF1.describe(include='all'))\n",
    "    print('')\n",
    "\n",
    "predictDF_F_LASSO = tempD1['Female']\n",
    "predictDF_M_LASSO = tempD1['Male']\n",
    "predictDF_B_LASSO = tempD1['BothSex']\n",
    "r2DF_F_LASSO = tempD2['Female']\n",
    "r2DF_M_LASSO = tempD2['Male']\n",
    "r2DF_B_LASSO = tempD2['BothSex']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> –> Surely reproduced, up to 15 decimal places. (The slight diffrence in the last decimal place was probably due to floating issue.)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1-2. Elastic net"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Elastic net\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-ElasticNet_'\n",
    "yvar = 'BMI'\n",
    "tempL1 = ['Met'+yvar, 'Prot'+yvar, 'Chem'+yvar, 'Combi'+yvar]\n",
    "tempL2 = ['Testing', 'log_Base'+yvar, 'Base'+yvar]\n",
    "tempD1 = {}\n",
    "tempD2 = {}\n",
    "for sex in ['Female', 'Male', 'BothSex']:\n",
    "    #Clean all mesured and predicted values into a dataframe\n",
    "    for yvar_model in tempL1:\n",
    "        #Import the predictions dataframe\n",
    "        fileName = yvar_model+'-'+sex+'.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",
    "        #Put aside the common part\n",
    "        if yvar_model==tempL1[0]:\n",
    "            tempDF1 = tempDF[tempL2]#Initialize\n",
    "        tempDF = tempDF.drop(columns=tempL2)\n",
    "        #Merge\n",
    "        tempDF1 = pd.merge(tempDF1, tempDF, left_index=True, right_index=True, how='inner')\n",
    "    ##Change column names\n",
    "    tempDF1.columns = tempDF1.columns.str.replace('Base', '')\n",
    "    \n",
    "    #Reproduce out-of-sample R2\n",
    "    tempL = tempDF1['Testing'].unique().tolist()\n",
    "    tempDF2 = pd.DataFrame(index=pd.Index(tempL, name='Model'))\n",
    "    for yvar_model in tempL1:\n",
    "        tempL = []\n",
    "        for model_n in tempDF2.index.tolist():\n",
    "            tempS1 = tempDF1['log_'+yvar].loc[tempDF1['Testing']==model_n]\n",
    "            tempS2 = tempDF1['log_'+yvar_model].loc[tempDF1['Testing']==model_n]\n",
    "            tempL.append(r2_score(tempS1, tempS2, sample_weight=None))\n",
    "        tempDF2[yvar_model] = tempL\n",
    "    \n",
    "    tempD1[sex] = tempDF1\n",
    "    tempD2[sex] = tempDF2\n",
    "    \n",
    "    #Check\n",
    "    print(sex)\n",
    "    ##Reproduce the previous outputs (because R2 values are rounded in the following summary table)\n",
    "    for yvar_model in tempL1:\n",
    "        tempL = tempDF2[yvar_model].tolist()\n",
    "        print(' - '+yvar_model)\n",
    "        print('   - Out-of-sample R2 [Mean ± SD]:',\n",
    "              np.mean(tempL), '±', np.std(tempL, ddof=1))#Sample standard deviation\n",
    "        print('   - Out-of-sample R2 [Mean ± SEM]:',\n",
    "              np.mean(tempL), '±', np.std(tempL, ddof=1)/np.sqrt(len(tempL)))\n",
    "        print('   - Observed vs. predicted log_'+yvar+': (Pearson\\'s r, P) =',\n",
    "              stats.pearsonr(tempDF1['log_'+yvar], tempDF1['log_'+yvar_model]))\n",
    "        print('   - Observed vs. predicted '+yvar+': (Pearson\\'s r, P) =',\n",
    "              stats.pearsonr(tempDF1[yvar], tempDF1[yvar_model]))\n",
    "    print(' - Summary of out-of-sample R2 values')\n",
    "    display(tempDF2.describe(include='all'))\n",
    "    print(' - Summary of all measured and predicted values')\n",
    "    display(tempDF1.describe(include='all'))\n",
    "    print('')\n",
    "\n",
    "predictDF_F_EN = tempD1['Female']\n",
    "predictDF_M_EN = tempD1['Male']\n",
    "predictDF_B_EN = tempD1['BothSex']\n",
    "r2DF_F_EN = tempD2['Female']\n",
    "r2DF_M_EN = tempD2['Male']\n",
    "r2DF_B_EN = tempD2['BothSex']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> –> Surely reproduced, up to 15 decimal places. (The slight diffrence in the last decimal place was probably due to floating issue.)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1-3. Ridge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Ridge\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-Ridge_'\n",
    "yvar = 'BMI'\n",
    "tempL1 = ['Met'+yvar, 'Prot'+yvar, 'Chem'+yvar, 'Combi'+yvar]\n",
    "tempL2 = ['Testing', 'log_Base'+yvar, 'Base'+yvar]\n",
    "tempD1 = {}\n",
    "tempD2 = {}\n",
    "for sex in ['Female', 'Male', 'BothSex']:\n",
    "    #Clean all mesured and predicted values into a dataframe\n",
    "    for yvar_model in tempL1:\n",
    "        #Import the predictions dataframe\n",
    "        fileName = yvar_model+'-'+sex+'.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",
    "        #Put aside the common part\n",
    "        if yvar_model==tempL1[0]:\n",
    "            tempDF1 = tempDF[tempL2]#Initialize\n",
    "        tempDF = tempDF.drop(columns=tempL2)\n",
    "        #Merge\n",
    "        tempDF1 = pd.merge(tempDF1, tempDF, left_index=True, right_index=True, how='inner')\n",
    "    ##Change column names\n",
    "    tempDF1.columns = tempDF1.columns.str.replace('Base', '')\n",
    "    \n",
    "    #Reproduce out-of-sample R2\n",
    "    tempL = tempDF1['Testing'].unique().tolist()\n",
    "    tempDF2 = pd.DataFrame(index=pd.Index(tempL, name='Model'))\n",
    "    for yvar_model in tempL1:\n",
    "        tempL = []\n",
    "        for model_n in tempDF2.index.tolist():\n",
    "            tempS1 = tempDF1['log_'+yvar].loc[tempDF1['Testing']==model_n]\n",
    "            tempS2 = tempDF1['log_'+yvar_model].loc[tempDF1['Testing']==model_n]\n",
    "            tempL.append(r2_score(tempS1, tempS2, sample_weight=None))\n",
    "        tempDF2[yvar_model] = tempL\n",
    "    \n",
    "    tempD1[sex] = tempDF1\n",
    "    tempD2[sex] = tempDF2\n",
    "    \n",
    "    #Check\n",
    "    print(sex)\n",
    "    ##Reproduce the previous outputs (because R2 values are rounded in the following summary table)\n",
    "    for yvar_model in tempL1:\n",
    "        tempL = tempDF2[yvar_model].tolist()\n",
    "        print(' - '+yvar_model)\n",
    "        print('   - Out-of-sample R2 [Mean ± SD]:',\n",
    "              np.mean(tempL), '±', np.std(tempL, ddof=1))#Sample standard deviation\n",
    "        print('   - Out-of-sample R2 [Mean ± SEM]:',\n",
    "              np.mean(tempL), '±', np.std(tempL, ddof=1)/np.sqrt(len(tempL)))\n",
    "        print('   - Observed vs. predicted log_'+yvar+': (Pearson\\'s r, P) =',\n",
    "              stats.pearsonr(tempDF1['log_'+yvar], tempDF1['log_'+yvar_model]))\n",
    "        print('   - Observed vs. predicted '+yvar+': (Pearson\\'s r, P) =',\n",
    "              stats.pearsonr(tempDF1[yvar], tempDF1[yvar_model]))\n",
    "    print(' - Summary of out-of-sample R2 values')\n",
    "    display(tempDF2.describe(include='all'))\n",
    "    print(' - Summary of all measured and predicted values')\n",
    "    display(tempDF1.describe(include='all'))\n",
    "    print('')\n",
    "\n",
    "predictDF_F_ridge = tempD1['Female']\n",
    "predictDF_M_ridge = tempD1['Male']\n",
    "predictDF_B_ridge = tempD1['BothSex']\n",
    "r2DF_F_ridge = tempD2['Female']\n",
    "r2DF_M_ridge = tempD2['Male']\n",
    "r2DF_B_ridge = tempD2['BothSex']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> –> Surely reproduced, up to 15 decimal places. (The slight diffrence in the last decimal place was probably due to floating issue.)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1-4. Random forest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Random forest\n",
    "fileDir = './ImportData/'\n",
    "ipynbName = '220828_Multiomics-BMI-NatMed1stRevision_BMI-baseline-RF-wenceslaus_'\n",
    "yvar = 'BMI'\n",
    "tempL1 = ['Met'+yvar, 'Prot'+yvar, 'Chem'+yvar, 'Combi'+yvar]\n",
    "tempL2 = ['Testing', 'log_Base'+yvar, 'Base'+yvar]\n",
    "tempD1 = {}\n",
    "tempD2 = {}\n",
    "for sex in ['Female', 'Male', 'BothSex']:\n",
    "    #Clean all mesured and predicted values into a dataframe\n",
    "    for yvar_model in tempL1:\n",
    "        #Import the predictions dataframe\n",
    "        fileName = yvar_model+'-'+sex+'.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",
    "        #Put aside the common part\n",
    "        if yvar_model==tempL1[0]:\n",
    "            tempDF1 = tempDF[tempL2]#Initialize\n",
    "        tempDF = tempDF.drop(columns=tempL2)\n",
    "        #Merge\n",
    "        tempDF1 = pd.merge(tempDF1, tempDF, left_index=True, right_index=True, how='inner')\n",
    "    ##Change column names\n",
    "    tempDF1.columns = tempDF1.columns.str.replace('Base', '')\n",
    "    \n",
    "    #Reproduce out-of-sample R2\n",
    "    tempL = tempDF1['Testing'].unique().tolist()\n",
    "    tempDF2 = pd.DataFrame(index=pd.Index(tempL, name='Model'))\n",
    "    for yvar_model in tempL1:\n",
    "        tempL = []\n",
    "        for model_n in tempDF2.index.tolist():\n",
    "            tempS1 = tempDF1['log_'+yvar].loc[tempDF1['Testing']==model_n]\n",
    "            tempS2 = tempDF1['log_'+yvar_model].loc[tempDF1['Testing']==model_n]\n",
    "            tempL.append(r2_score(tempS1, tempS2, sample_weight=None))\n",
    "        tempDF2[yvar_model] = tempL\n",
    "    \n",
    "    tempD1[sex] = tempDF1\n",
    "    tempD2[sex] = tempDF2\n",
    "    \n",
    "    #Check\n",
    "    print(sex)\n",
    "    ##Reproduce the previous outputs (because R2 values are rounded in the following summary table)\n",
    "    for yvar_model in tempL1:\n",
    "        tempL = tempDF2[yvar_model].tolist()\n",
    "        print(' - '+yvar_model)\n",
    "        print('   - Out-of-sample R2 [Mean ± SD]:',\n",
    "              np.mean(tempL), '±', np.std(tempL, ddof=1))#Sample standard deviation\n",
    "        print('   - Out-of-sample R2 [Mean ± SEM]:',\n",
    "              np.mean(tempL), '±', np.std(tempL, ddof=1)/np.sqrt(len(tempL)))\n",
    "        print('   - Observed vs. predicted log_'+yvar+': (Pearson\\'s r, P) =',\n",
    "              stats.pearsonr(tempDF1['log_'+yvar], tempDF1['log_'+yvar_model]))\n",
    "        print('   - Observed vs. predicted '+yvar+': (Pearson\\'s r, P) =',\n",
    "              stats.pearsonr(tempDF1[yvar], tempDF1[yvar_model]))\n",
    "    print(' - Summary of out-of-sample R2 values')\n",
    "    display(tempDF2.describe(include='all'))\n",
    "    print(' - Summary of all measured and predicted values')\n",
    "    display(tempDF1.describe(include='all'))\n",
    "    print('')\n",
    "\n",
    "predictDF_F_RF = tempD1['Female']\n",
    "predictDF_M_RF = tempD1['Male']\n",
    "predictDF_B_RF = tempD1['BothSex']\n",
    "r2DF_F_RF = tempD2['Female']\n",
    "r2DF_M_RF = tempD2['Male']\n",
    "r2DF_B_RF = tempD2['BothSex']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> –> Surely reproduced, up to 15 decimal places. (The slight diffrence in the last decimal place was probably due to floating issue.)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Out-of-sample R2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempD1 = {'LASSO':r2DF_B_LASSO, 'Elastic net':r2DF_B_EN,\n",
    "          'Ridge':r2DF_B_ridge, 'Random forest':r2DF_B_RF}\n",
    "tempD2 = {'LASSO':'LASSO', 'Elastic net':'EN',\n",
    "          'Ridge':'Ridge', 'Random forest':'RF'}\n",
    "tempD3 = {'MetBMI':'Metabolomics', 'ProtBMI':'Proteomics',\n",
    "          'ChemBMI':'Clinical labs', 'CombiBMI':'Combined omics'}\n",
    "tempD4 = {'Metabolomics':'b', 'Proteomics':'r',\n",
    "          'Clinical labs':'g', 'Combined omics':'m'}\n",
    "\n",
    "#Prepare DF\n",
    "tempD = {}\n",
    "for method in tempD1.keys():\n",
    "    tempDF = tempD1[method].copy()\n",
    "    method_label = tempD2[method]\n",
    "    tempDF.columns = method_label+':'+tempDF.columns\n",
    "    tempD[method] = tempDF\n",
    "tempDF1 = pd.concat(list(tempD.values()), axis=1)\n",
    "tempL = [method_label+':'+yvar_model for yvar_model in tempD3.keys() for method_label in tempD2.values()]\n",
    "tempDF1 = tempDF1[tempL]#Sort just for easily checking the summary\n",
    "display(tempDF1.describe())\n",
    "\n",
    "#Statistical tests\n",
    "tempD = {}\n",
    "for yvar_model in tempD3.keys():\n",
    "    control = list(tempD2.values())[0]\n",
    "    tempL = list(tempD2.values())[1:]\n",
    "    tempDF = pd.DataFrame(columns=['Control', 'Contrast', 'Control_N', 'contrast_N', 'DoF', 'tStat', 'Pval'])\n",
    "    for contrast in tempL:\n",
    "        tempS1 = tempDF1[control+':'+yvar_model]\n",
    "        tempS2 = tempDF1[contrast+':'+yvar_model]\n",
    "        #Two-sided Welch's t-test\n",
    "        tstat, pval, dof = weightstats.ttest_ind(tempS1, tempS2,\n",
    "                                                 alternative='two-sided', usevar='unequal')\n",
    "        size1 = len(tempS1)\n",
    "        size2 = len(tempS2)\n",
    "        tempDF.loc[contrast+'-vs-'+control] = [control, contrast, size1, size2, dof, tstat, pval]\n",
    "    ##P-value adjustment (within omics) by using Benjamini–Hochberg method\n",
    "    tempDF['AdjPval_omics'] = multi.multipletests(tempDF['Pval'], alpha=0.05, method='fdr_bh',\n",
    "                                                  is_sorted=False, returnsorted=False)[1]\n",
    "    tempDF['Category'] = tempD3[yvar_model]\n",
    "    tempD[yvar_model] = tempDF\n",
    "tempDF2 = pd.concat(list(tempD.values()), axis=0)\n",
    "##P-value adjustment (across all tests) by using Benjamini–Hochberg method\n",
    "tempDF2['AdjPval_all'] = multi.multipletests(tempDF2['Pval'], alpha=0.05, method='fdr_bh',\n",
    "                                             is_sorted=False, returnsorted=False)[1]\n",
    "tempDF2.index.rename('ComparisonLabel', inplace=True)\n",
    "tempDF2 = tempDF2.reset_index().set_index(['Category', 'ComparisonLabel'])\n",
    "display(tempDF2)\n",
    "##Save\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220901_Multiomics-BMI-NatMed1stRevision_bBMI-LASSO-vs-others_'\n",
    "fileName = 'R2-comparison-BothSex.tsv'\n",
    "tempDF2.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "\n",
    "#Prepare a dictionary for inverse mapping\n",
    "tempD = {}\n",
    "for key in tempD2.keys():\n",
    "    value = tempD2[key]\n",
    "    tempD[value] = key\n",
    "\n",
    "#Visualization\n",
    "axis_ymin = 0.0\n",
    "axis_ymax = 1.15\n",
    "ymin = 0.0\n",
    "ymax = 0.8\n",
    "yinter = 0.2\n",
    "aline_ymin = 0.825\n",
    "aline_yinter = 0.1\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "fig, axes = plt.subplots(nrows=1, ncols=len(tempD3),\n",
    "                         figsize=(1+2*len(tempD3), 3), sharex=True, sharey=True)\n",
    "for ax_i, ax in enumerate(axes.flat):\n",
    "    yvar_model = list(tempD3.keys())[ax_i]\n",
    "    category = tempD3[yvar_model]\n",
    "    #Prepare DF while cleaning label text\n",
    "    tempDF = tempDF1.loc[:, tempDF1.columns.str.contains(yvar_model)]\n",
    "    tempDF.columns = tempDF.columns.str.replace(':'+yvar_model, '')\n",
    "    tempDF = tempDF.melt(var_name='Method', value_name='R2', value_vars=tempDF.columns.tolist())\n",
    "    tempDF['Method'] = tempDF['Method'].map(tempD)\n",
    "    #Plot\n",
    "    sns.barplot(data=tempDF, y='R2', x='Method', order=tempD2.keys(),\n",
    "                color=tempD4[category], dodge=False,\n",
    "                ci=95, capsize=0.4, errwidth=1.5, errcolor='black', edgecolor='black', ax=ax)\n",
    "    sns.stripplot(data=tempDF, y='R2', x='Method', order=tempD2.keys(),\n",
    "                  dodge=False, jitter=0.3,\n",
    "                  size=5, edgecolor='black', color='gray', linewidth=1, alpha=0.4, ax=ax)\n",
    "    #P-value annotation\n",
    "    tempDF = tempDF2.loc[category]#MultiIndex\n",
    "    for row_i in range(len(tempDF)):\n",
    "        #Control\n",
    "        group_0 = tempDF['Control'].iloc[row_i]\n",
    "        xcoord_0 = list(tempD2.keys()).index(tempD[group_0])\n",
    "        #Contrast\n",
    "        group_1 = tempDF['Contrast'].iloc[row_i]\n",
    "        xcoord_1 = list(tempD2.keys()).index(tempD[group_1])\n",
    "        #Standard point of marker\n",
    "        xcoord = (xcoord_0+xcoord_1)/2\n",
    "        ycoord = aline_ymin + aline_yinter*row_i\n",
    "        #Add annotation lines\n",
    "        aline_offset = yinter/10\n",
    "        aline_length = yinter/10 + aline_offset\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",
    "        pval = tempDF['AdjPval_all'].iloc[row_i]\n",
    "        if pval<0.001:\n",
    "            label = '***'\n",
    "        elif pval<0.01:\n",
    "            label = '**'\n",
    "        elif pval<0.05:\n",
    "            label = '*'\n",
    "        else:\n",
    "            pval_text = str(Decimal(pval).quantize(Decimal('0.01'), rounding=ROUND_HALF_UP))\n",
    "            label = r'$P$'+' = '+pval_text\n",
    "        #Add annotation text\n",
    "        if label in ['***', '**', '*']:\n",
    "            text_offset = yinter/25\n",
    "            text_size = 'medium'\n",
    "        else:\n",
    "            text_offset = yinter/5\n",
    "            text_size = 'x-small'\n",
    "        ax.annotate(label, xy=(xcoord, ycoord+text_offset),\n",
    "                    horizontalalignment='center', verticalalignment='bottom',\n",
    "                    fontsize=text_size, color='k')\n",
    "    #Facet label\n",
    "    ax.set_title(category, {'fontsize':'medium'})\n",
    "    #Axis setting\n",
    "    plt.setp(ax.get_xticklabels(), rotation=70,\n",
    "             horizontalalignment='right', verticalalignment='center', rotation_mode='anchor')\n",
    "    if ax_i==0:\n",
    "        plt.setp(ax, xlabel='', ylabel='Out-of-sample '+r'$R^2$')\n",
    "    else:\n",
    "        plt.setp(ax, xlabel='', ylabel='')\n",
    "        plt.setp(ax.get_yticklabels(), visible=False)\n",
    "sns.despine()\n",
    "plt.setp(axes, ylim=(axis_ymin, axis_ymax), yticks=np.arange(ymin, ymax+yinter/10, yinter))\n",
    "##Save\n",
    "fileDir = './ExportFigures/'\n",
    "ipynbName = '220901_Multiomics-BMI-NatMed1stRevision_bBMI-LASSO-vs-others_'\n",
    "fileName = 'R2-comparison-BothSex.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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempD1 = {'LASSO':r2DF_F_LASSO, 'Elastic net':r2DF_F_EN,\n",
    "          'Ridge':r2DF_F_ridge, 'Random forest':r2DF_F_RF}\n",
    "tempD2 = {'LASSO':r2DF_M_LASSO, 'Elastic net':r2DF_M_EN,\n",
    "          'Ridge':r2DF_M_ridge, 'Random forest':r2DF_M_RF}\n",
    "tempD1 = {'Female':tempD1, 'Male':tempD2}\n",
    "tempD2 = {'LASSO':'LASSO', 'Elastic net':'EN',\n",
    "          'Ridge':'Ridge', 'Random forest':'RF'}\n",
    "tempD3 = {'MetBMI':'Metabolomics', 'ProtBMI':'Proteomics',\n",
    "          'ChemBMI':'Clinical labs', 'CombiBMI':'Combined omics'}\n",
    "tempD4 = {'Metabolomics':'b', 'Proteomics':'r',\n",
    "          'Clinical labs':'g', 'Combined omics':'m'}\n",
    "\n",
    "#Prepare DF\n",
    "tempD5 = {}\n",
    "for sex in tempD1.keys():\n",
    "    tempD6 = tempD1[sex]\n",
    "    tempD = {}\n",
    "    for method in tempD6.keys():\n",
    "        tempDF = tempD6[method].copy()\n",
    "        method_label = tempD2[method]\n",
    "        tempDF.columns = method_label+':'+tempDF.columns\n",
    "        tempD[method] = tempDF\n",
    "    tempDF1 = pd.concat(list(tempD.values()), axis=1)\n",
    "    tempDF1['Sex'] = sex\n",
    "    tempD5[sex] = tempDF1\n",
    "tempDF1 = pd.concat(list(tempD5.values()), axis=0)\n",
    "tempDF1 = tempDF1.reset_index().set_index(['Sex', 'Model'])\n",
    "tempL = [method_label+':'+yvar_model for yvar_model in tempD3.keys() for method_label in tempD2.values()]\n",
    "tempDF1 = tempDF1[tempL]#Sort just for easily checking the summary\n",
    "for sex in tempD1.keys():#To display all summary in the Jupyter Notebook\n",
    "    print(sex)\n",
    "    display(tempDF1.loc[sex].describe())\n",
    "\n",
    "#Statistical tests\n",
    "tempD5 = {}\n",
    "for sex in tempD1.keys():\n",
    "    tempDF2 = tempDF1.loc[sex]#MultiIndex\n",
    "    tempD = {}\n",
    "    for yvar_model in tempD3.keys():\n",
    "        control = list(tempD2.values())[0]\n",
    "        tempL = list(tempD2.values())[1:]\n",
    "        tempDF = pd.DataFrame(columns=['Control', 'Contrast', 'Control_N', 'contrast_N', 'DoF', 'tStat', 'Pval'])\n",
    "        for contrast in tempL:\n",
    "            tempS1 = tempDF2[control+':'+yvar_model]\n",
    "            tempS2 = tempDF2[contrast+':'+yvar_model]\n",
    "            #Two-sided Welch's t-test\n",
    "            tstat, pval, dof = weightstats.ttest_ind(tempS1, tempS2,\n",
    "                                                     alternative='two-sided', usevar='unequal')\n",
    "            size1 = len(tempS1)\n",
    "            size2 = len(tempS2)\n",
    "            tempDF.loc[contrast+'-vs-'+control] = [control, contrast, size1, size2, dof, tstat, pval]\n",
    "        ##P-value adjustment (within sex and omics) by using Benjamini–Hochberg method\n",
    "        tempDF['AdjPval_sex-omics'] = multi.multipletests(tempDF['Pval'], alpha=0.05, method='fdr_bh',\n",
    "                                                          is_sorted=False, returnsorted=False)[1]\n",
    "        tempDF['Category'] = tempD3[yvar_model]\n",
    "        tempD[yvar_model] = tempDF\n",
    "    tempDF = pd.concat(list(tempD.values()), axis=0)\n",
    "    ##P-value adjustment (within sex) by using Benjamini–Hochberg method\n",
    "    tempDF['AdjPval_sex'] = multi.multipletests(tempDF['Pval'], alpha=0.05, method='fdr_bh',\n",
    "                                                is_sorted=False, returnsorted=False)[1]\n",
    "    tempDF['Sex'] = sex\n",
    "    tempD5[sex] = tempDF\n",
    "tempDF2 = pd.concat(list(tempD5.values()), axis=0)\n",
    "##P-value adjustment (across all tests) by using Benjamini–Hochberg method\n",
    "tempDF2['AdjPval_all'] = multi.multipletests(tempDF2['Pval'], alpha=0.05, method='fdr_bh',\n",
    "                                             is_sorted=False, returnsorted=False)[1]\n",
    "tempDF2.index.rename('ComparisonLabel', inplace=True)\n",
    "tempDF2 = tempDF2.reset_index().set_index(['Sex', 'Category', 'ComparisonLabel'])\n",
    "display(tempDF2)\n",
    "##Save\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220901_Multiomics-BMI-NatMed1stRevision_bBMI-LASSO-vs-others_'\n",
    "fileName = 'R2-comparison-FemaleMale.tsv'\n",
    "tempDF2.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "\n",
    "#Prepare a dictionary for inverse mapping\n",
    "tempD = {}\n",
    "for key in tempD2.keys():\n",
    "    value = tempD2[key]\n",
    "    tempD[value] = key\n",
    "\n",
    "#Visualization (per sex, for now)\n",
    "for sex in tempD1.keys():\n",
    "    axis_ymin = 0.0\n",
    "    axis_ymax = 1.15\n",
    "    ymin = 0.0\n",
    "    ymax = 0.8\n",
    "    yinter = 0.2\n",
    "    aline_ymin = 0.825\n",
    "    aline_yinter = 0.1\n",
    "    sns.set(style='ticks', font='Arial', context='talk')\n",
    "    fig, axes = plt.subplots(nrows=1, ncols=len(tempD3),\n",
    "                             figsize=(1+2*len(tempD3), 3), sharex=True, sharey=True)\n",
    "    for ax_i, ax in enumerate(axes.flat):\n",
    "        yvar_model = list(tempD3.keys())[ax_i]\n",
    "        category = tempD3[yvar_model]\n",
    "        #Prepare DF while cleaning label text\n",
    "        tempDF = tempDF1.loc[sex]#MultiIndex\n",
    "        tempDF = tempDF.loc[:, tempDF.columns.str.contains(yvar_model)]\n",
    "        tempDF.columns = tempDF.columns.str.replace(':'+yvar_model, '')\n",
    "        tempDF = tempDF.melt(var_name='Method', value_name='R2', value_vars=tempDF.columns.tolist())\n",
    "        tempDF['Method'] = tempDF['Method'].map(tempD)\n",
    "        #Plot\n",
    "        sns.barplot(data=tempDF, y='R2', x='Method', order=tempD2.keys(),\n",
    "                    color=tempD4[category], dodge=False,\n",
    "                    ci=95, capsize=0.4, errwidth=1.5, errcolor='black', edgecolor='black', ax=ax)\n",
    "        sns.stripplot(data=tempDF, y='R2', x='Method', order=tempD2.keys(),\n",
    "                      dodge=False, jitter=0.3,\n",
    "                      size=5, edgecolor='black', color='gray', linewidth=1, alpha=0.4, ax=ax)\n",
    "        #P-value annotation\n",
    "        tempDF = tempDF2.loc[sex]#MultiIndex\n",
    "        tempDF = tempDF.loc[category]#MultiIndex\n",
    "        for row_i in range(len(tempDF)):\n",
    "            #Control\n",
    "            group_0 = tempDF['Control'].iloc[row_i]\n",
    "            xcoord_0 = list(tempD2.keys()).index(tempD[group_0])\n",
    "            #Contrast\n",
    "            group_1 = tempDF['Contrast'].iloc[row_i]\n",
    "            xcoord_1 = list(tempD2.keys()).index(tempD[group_1])\n",
    "            #Standard point of marker\n",
    "            xcoord = (xcoord_0+xcoord_1)/2\n",
    "            ycoord = aline_ymin + aline_yinter*row_i\n",
    "            #Add annotation lines\n",
    "            aline_offset = yinter/10\n",
    "            aline_length = yinter/10 + aline_offset\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",
    "            pval = tempDF['AdjPval_all'].iloc[row_i]\n",
    "            if pval<0.001:\n",
    "                label = '***'\n",
    "            elif pval<0.01:\n",
    "                label = '**'\n",
    "            elif pval<0.05:\n",
    "                label = '*'\n",
    "            else:\n",
    "                pval_text = str(Decimal(pval).quantize(Decimal('0.01'), rounding=ROUND_HALF_UP))\n",
    "                label = r'$P$'+' = '+pval_text\n",
    "            #Add annotation text\n",
    "            if label in ['***', '**', '*']:\n",
    "                text_offset = yinter/25\n",
    "                text_size = 'medium'\n",
    "            else:\n",
    "                text_offset = yinter/5\n",
    "                text_size = 'x-small'\n",
    "            ax.annotate(label, xy=(xcoord, ycoord+text_offset),\n",
    "                        horizontalalignment='center', verticalalignment='bottom',\n",
    "                        fontsize=text_size, color='k')\n",
    "        #Facet label\n",
    "        ax.set_title(category, {'fontsize':'medium'})\n",
    "        #Axis setting\n",
    "        plt.setp(ax.get_xticklabels(), rotation=70,\n",
    "                 horizontalalignment='right', verticalalignment='center', rotation_mode='anchor')\n",
    "        if ax_i==0:\n",
    "            plt.setp(ax, xlabel='', ylabel='Out-of-sample '+r'$R^2$')\n",
    "        else:\n",
    "            plt.setp(ax, xlabel='', ylabel='')\n",
    "            plt.setp(ax.get_yticklabels(), visible=False)\n",
    "    sns.despine()\n",
    "    plt.setp(axes, ylim=(axis_ymin, axis_ymax), yticks=np.arange(ymin, ymax+yinter/10, yinter))\n",
    "    fig.suptitle(sex, size='medium',\n",
    "                 verticalalignment='bottom', horizontalalignment='center', y=0.965)\n",
    "    ##Save\n",
    "    fileDir = './ExportFigures/'\n",
    "    ipynbName = '220901_Multiomics-BMI-NatMed1stRevision_bBMI-LASSO-vs-others_'\n",
    "    fileName = 'R2-comparison-'+sex+'.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": [
    "## 3. Predictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempD1 = {'LASSO':predictDF_B_LASSO, 'Elastic net':predictDF_B_EN,\n",
    "          'Ridge':predictDF_B_ridge, 'Random forest':predictDF_B_RF}\n",
    "tempD2 = {'Metabolomics':'MetBMI', 'Proteomics':'ProtBMI',\n",
    "          'Clinical labs':'ChemBMI', 'Combined omics':'CombiBMI'}\n",
    "tempD3 = {'Metabolomics':'b', 'Proteomics':'r',\n",
    "          'Clinical labs':'g', 'Combined omics':'m'}\n",
    "yvar = 'LASSO'\n",
    "unit_label = ' [kg m'+r'$^{-2}$'+']'\n",
    "\n",
    "#Prepare DFs per category\n",
    "tempD4 = {}\n",
    "for category in tempD2.keys():\n",
    "    bbmi = tempD2[category]\n",
    "    tempD = {}\n",
    "    for method in tempD1.keys():\n",
    "        tempDF = tempD1[method]\n",
    "        tempS = tempDF[bbmi].copy()\n",
    "        tempS.name = method\n",
    "        tempD[method] = tempS\n",
    "    tempDF = pd.concat(list(tempD.values()), axis=1)\n",
    "    print(category)\n",
    "    display(tempDF.describe())\n",
    "    tempD4[category] = tempDF\n",
    "\n",
    "#Statistical tests\n",
    "tempD = {}\n",
    "for category in tempD4.keys():\n",
    "    tempDF = tempD4[category]\n",
    "    tempDF1 = pd.DataFrame(columns=['Yvar', 'N', 'DoF', 'Pearson_r', 'Pval'])\n",
    "    for method in tempD1.keys():\n",
    "        if method!=yvar:\n",
    "            xvar = method\n",
    "            #Pearson's correlation\n",
    "            pearson_r, pval = stats.pearsonr(tempDF[xvar], tempDF[yvar])\n",
    "            size = len(tempDF)\n",
    "            dof = size - 2\n",
    "            tempDF1.loc[method] = [yvar, size, dof, pearson_r, pval]\n",
    "    ##P-value adjustment (within omics) by using Benjamini–Hochberg method\n",
    "    tempDF1['AdjPval_omics'] = multi.multipletests(tempDF1['Pval'], alpha=0.05, method='fdr_bh',\n",
    "                                                   is_sorted=False, returnsorted=False)[1]\n",
    "    tempDF1.index.rename('Xvar', inplace=True)\n",
    "    tempDF1['Category'] = category\n",
    "    tempD[category] = tempDF1\n",
    "tempDF1 = pd.concat(list(tempD.values()), axis=0)\n",
    "##P-value adjustment (across all tests) by using Benjamini–Hochberg method\n",
    "tempDF1['AdjPval_all'] = multi.multipletests(tempDF1['Pval'], alpha=0.05, method='fdr_bh',\n",
    "                                             is_sorted=False, returnsorted=False)[1]\n",
    "tempDF1['N'] = tempDF1['N'].astype('int64')#Otherwise, float64!\n",
    "tempDF1['DoF'] = tempDF1['DoF'].astype('int64')#Otherwise, float64!\n",
    "tempDF1 = tempDF1.reset_index().set_index(['Category', 'Xvar'])\n",
    "display(tempDF1)\n",
    "##Save\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220901_Multiomics-BMI-NatMed1stRevision_bBMI-LASSO-vs-others_'\n",
    "fileName = 'bBMI-comparison-BothSex.tsv'\n",
    "tempDF1.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "\n",
    "#Visualization (per category, for now)\n",
    "for category in tempD4.keys():\n",
    "    sns.set(style='ticks', font='Arial', context='talk')\n",
    "    fig, axes = plt.subplots(nrows=1, ncols=len(tempD1)-1,\n",
    "                             figsize=(2.8*(len(tempD1)-1), 2.8+0.8), sharex=True, sharey=True)\n",
    "    axis_xymin = 12.5\n",
    "    axis_xymax = 57.5\n",
    "    xymin = 20\n",
    "    xymax = 50\n",
    "    xyinter = 10\n",
    "    #Set axis range first; otherwise, regression line can be truncated differently\n",
    "    plt.setp(axes, xlim=(axis_xymin, axis_xymax), xticks=np.arange(xymin, xymax+xyinter/10, xyinter))\n",
    "    plt.setp(axes, ylim=(axis_xymin, axis_xymax), yticks=np.arange(xymin, xymax+xyinter/10, xyinter))\n",
    "    for ax_i, ax in enumerate(axes.flat):\n",
    "        xvar = list(tempD1.keys())[ax_i+1]\n",
    "        tempDF = tempD4[category]\n",
    "        #Scatterplot with regression line\n",
    "        sns.regplot(data=tempDF, x=xvar, y=yvar, color=tempD3[category],\n",
    "                    scatter=True, fit_reg=True, ci=95, truncate=False, marker='o',\n",
    "                    scatter_kws={'alpha':0.2, 'edgecolor':'k', 's':25}, ax=ax)\n",
    "        #Draw Y=X as reference\n",
    "        ax.plot([axis_xymin, axis_xymax], [axis_xymin, axis_xymax],\n",
    "                color='black', linestyle=(0, (1, 2)), zorder=0)\n",
    "        #Annotate Pearson's correlation\n",
    "        tempDF = tempDF1.loc[category]#MultiIndex\n",
    "        pearson_r = tempDF['Pearson_r'].loc[xvar]\n",
    "        r_text = str(Decimal(str(pearson_r)).quantize(Decimal('0.001'), rounding=ROUND_HALF_UP))\n",
    "        pval = tempDF['AdjPval_all'].loc[xvar]\n",
    "        below_limit = 0#Initialize\n",
    "        if pval==1.0:\n",
    "            pval_text = '1.0'\n",
    "        else:\n",
    "            if pval==0.0:#Due to smaller than the float minimum\n",
    "                pval = sys.float_info.min\n",
    "                print('P-value was smaller than the float minimum:', pval)\n",
    "                below_limit = 1\n",
    "            pval_text = f'{Decimal(str(pval)):.3E}'#Take more digits because rounding is bad here\n",
    "            significand, exponent = pval_text.split(sep='E-')\n",
    "            significand = str(Decimal(significand).quantize(Decimal('0.1'), rounding=ROUND_HALF_UP))\n",
    "            if significand=='10.0':\n",
    "                significand = '1.0'\n",
    "                exponent = str(int(exponent)-1)\n",
    "            if int(exponent)>2:\n",
    "                pval_text = significand+r'$\\times$'+'10'+r'$^{{-{0}}}$'.format(exponent)##Font is different in r'$ $'...\n",
    "            elif int(exponent)>0:\n",
    "                pval_text = '0.'+'0'*(int(exponent)-1)+significand.replace('.', '')\n",
    "            else:\n",
    "                pval_text = significand\n",
    "        if below_limit==1:\n",
    "            text = 'Pearson\\'s '+r'$r$'+' = '+r_text+'\\n'+r'$P$'+' < '+pval_text\n",
    "        else:\n",
    "            text = 'Pearson\\'s '+r'$r$'+' = '+r_text+'\\n'+r'$P$'+' = '+pval_text\n",
    "        ax.annotate(text, xy=(0.05, 0.95), xycoords='axes fraction',\n",
    "                    horizontalalignment='left', verticalalignment='top',\n",
    "                    multialignment='left', fontsize='small', color='k')\n",
    "        #Facet label (Fig title)\n",
    "        if ax_i==np.median(range(len(tempD1)-1)):\n",
    "            title = category+': '+tempD2[category]+unit_label\n",
    "            ax.set_title(title, {'fontsize':'large'})\n",
    "        #Axis setting\n",
    "        if ax_i%(len(tempD1)-1)==0:\n",
    "            plt.setp(ax, xlabel=xvar, ylabel=yvar)\n",
    "        else:\n",
    "            plt.setp(ax, xlabel=xvar, ylabel='')\n",
    "            plt.setp(ax.get_yticklabels(), visible=False)\n",
    "    sns.despine()\n",
    "    #Reset and generate common axis title\n",
    "    #plt.setp(axes, xlabel='', ylabel='')\n",
    "    fig.tight_layout(pad=0.75)\n",
    "    ##Save\n",
    "    fileDir = './ExportFigures/'\n",
    "    ipynbName = '220901_Multiomics-BMI-NatMed1stRevision_bBMI-LASSO-vs-others_'\n",
    "    fileName = tempD2[category]+'-comparison-BothSex.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. Variables between LASSO and EN models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4-1. Import beta-coefficients"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#LASSO\n",
    "print('LASSO models')\n",
    "tempD = {}\n",
    "for bbmi in ['MetBMI', 'ProtBMI', 'ChemBMI', 'CombiBMI']:\n",
    "    #Import the beta-coefficients\n",
    "    fileDir = './ExportData/'\n",
    "    ipynbName = '220801_Multiomics-BMI-NatMed1stRevision_BMI-baseline-LASSO_'\n",
    "    fileName = bbmi+'-BothSex-LASSObcoefs.tsv'\n",
    "    tempDF = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t').set_index('Variable')\n",
    "    tempDF = tempDF.drop(index=['Intercept'])\n",
    "    tempD[bbmi] = tempDF\n",
    "    \n",
    "    #Check\n",
    "    print(bbmi+':')\n",
    "    print(' - Variables:', len(tempDF))\n",
    "    #Variables with non-zero beta-coefficient\n",
    "    tempDF1 = tempDF.loc[tempDF['nZeros']!=10]\n",
    "    print(' - Variables with non-zero beta-coefficient:', len(tempDF1),\n",
    "          '(', len(tempDF1)/len(tempDF)*100, '%)')\n",
    "    #Extract robust beta-coefficient: no zeros in all 10 models\n",
    "    tempDF1 = tempDF.loc[tempDF['nZeros']==0]\n",
    "    tempDF1 = tempDF1.sort_values(by='Mean', ascending=False)\n",
    "    print(' - Variables with non-zero beta-coefficient in all 10 models:', len(tempDF1),\n",
    "          '(', len(tempDF1)/len(tempDF)*100, '%)')\n",
    "    print('')\n",
    "\n",
    "bcoefDF_met_LASSO = tempD['MetBMI']\n",
    "bcoefDF_prot_LASSO = tempD['ProtBMI']\n",
    "bcoefDF_chem_LASSO = tempD['ChemBMI']\n",
    "bcoefDF_combi_LASSO = tempD['CombiBMI']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Elastic net\n",
    "print('Elastic net models')\n",
    "tempD = {}\n",
    "for bbmi in ['MetBMI', 'ProtBMI', 'ChemBMI', 'CombiBMI']:\n",
    "    #Import the beta-coefficients\n",
    "    fileDir = './ExportData/'\n",
    "    ipynbName = '220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-ElasticNet_'\n",
    "    fileName = bbmi+'-BothSex-ElasticNetbcoefs.tsv'\n",
    "    tempDF = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t').set_index('Variable')\n",
    "    tempDF = tempDF.drop(index=['Intercept'])\n",
    "    tempD[bbmi] = tempDF\n",
    "    \n",
    "    #Check\n",
    "    print(bbmi+':')\n",
    "    print(' - Variables:', len(tempDF))\n",
    "    #Variables with non-zero beta-coefficient\n",
    "    tempDF1 = tempDF.loc[tempDF['nZeros']!=10]\n",
    "    print(' - Variables with non-zero beta-coefficient:', len(tempDF1),\n",
    "          '(', len(tempDF1)/len(tempDF)*100, '%)')\n",
    "    #Extract robust beta-coefficient: no zeros in all 10 models\n",
    "    tempDF1 = tempDF.loc[tempDF['nZeros']==0]\n",
    "    tempDF1 = tempDF1.sort_values(by='Mean', ascending=False)\n",
    "    print(' - Variables with non-zero beta-coefficient in all 10 models:', len(tempDF1),\n",
    "          '(', len(tempDF1)/len(tempDF)*100, '%)')\n",
    "    print('')\n",
    "\n",
    "bcoefDF_met_EN = tempD['MetBMI']\n",
    "bcoefDF_prot_EN = tempD['ProtBMI']\n",
    "bcoefDF_chem_EN = tempD['ChemBMI']\n",
    "bcoefDF_combi_EN = tempD['CombiBMI']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 4-2. Metabolites"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Prepare target analytes\n",
    "tempDF1 = bcoefDF_combi_LASSO.loc[bcoefDF_met_LASSO.index.tolist()]\n",
    "tempDF2 = bcoefDF_combi_EN.loc[bcoefDF_met_EN.index.tolist()]\n",
    "tempD1 = {'LASSO CombiBMI':tempDF1, 'LASSO MetBMI':bcoefDF_met_LASSO,\n",
    "          'Elastic net MetBMI':bcoefDF_met_EN, 'Elastic net CombiBMI':tempDF2}\n",
    "tempD2 = {'LASSO CombiBMI':'tab:blue', 'LASSO MetBMI':'tab:orange',\n",
    "          'Elastic net MetBMI':'tab:green', 'Elastic net CombiBMI':'tab:red'}\n",
    "analyte = 'metabolites'\n",
    "title = 'Robustly retained metabolites'\n",
    "\n",
    "#Prepare module sets\n",
    "print('Variables with non-zero beta-coefficient in all 10 models')\n",
    "tempD = {}\n",
    "for measure in tempD1.keys():\n",
    "    tempDF = tempD1[measure]\n",
    "    #Extract robust beta-coefficient: no zeros in all 10 models\n",
    "    tempDF1 = tempDF.loc[tempDF['nZeros']==0]\n",
    "    print(measure+':', len(tempDF1), 'per', len(tempDF), analyte+' (',\n",
    "          len(tempDF1)/len(tempDF)*100, '%)')\n",
    "    tempD[measure] = set(tempDF1.index.tolist())\n",
    "\n",
    "#Check common region (1,1,1,1)\n",
    "tempS = list(tempD.values())[0]#Initialize\n",
    "for measure in tempD.keys():\n",
    "    tempS = tempS & tempD[measure]\n",
    "print(' -> Common (1,1,1,1):', len(tempS))\n",
    "display(tempS)\n",
    "\n",
    "#Venn diagram\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "fig, ax = plt.subplots(figsize=(4, 4))\n",
    "venn(tempD, fmt='{size:,}', cmap=list(tempD2.values()), legend_loc=None, ax=ax)\n",
    "plt.setp(ax, ylim=(0.1, 0.875))#Otherwise, weird space...\n",
    "##Add legend annotation\n",
    "x_coord = [0.1, 0.1, 0.9, 0.9]\n",
    "y_coord = [0.25, 0.7, 0.7, 0.25]\n",
    "h_align = ['right', 'right', 'left', 'left']\n",
    "v_align = ['top', 'bottom', 'bottom', 'top']\n",
    "for i in range(len(tempD1)):\n",
    "    key = list(tempD2.keys())[i]\n",
    "    total = f'{len(tempD[key]):,}'\n",
    "    ax.text(x_coord[i], y_coord[i], key+'\\n('+total+' '+analyte+')',\n",
    "            fontsize='small', multialignment='center',\n",
    "            horizontalalignment=h_align[i], verticalalignment=v_align[i],\n",
    "            bbox={'boxstyle':'round', 'facecolor':tempD2[key], 'pad':0.2, 'alpha':0.5})\n",
    "ax.set_title(title, fontsize='medium')\n",
    "##Save\n",
    "fileDir = './ExportFigures/'\n",
    "ipynbName = '220901_Multiomics-BMI-NatMed1stRevision_bBMI-LASSO-vs-others_'\n",
    "fileName = 'LASSO-vs-ElasticNet-'+analyte+'.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": {
    "tags": []
   },
   "source": [
    "### 4-3. Proteins"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Prepare target analytes\n",
    "tempDF1 = bcoefDF_combi_LASSO.loc[bcoefDF_prot_LASSO.index.tolist()]\n",
    "tempDF2 = bcoefDF_combi_EN.loc[bcoefDF_prot_EN.index.tolist()]\n",
    "tempD1 = {'LASSO CombiBMI':tempDF1, 'LASSO ProtBMI':bcoefDF_prot_LASSO,\n",
    "          'Elastic net ProtBMI':bcoefDF_prot_EN, 'Elastic net CombiBMI':tempDF2}\n",
    "tempD2 = {'LASSO CombiBMI':'tab:blue', 'LASSO ProtBMI':'tab:orange',\n",
    "          'Elastic net ProtBMI':'tab:green', 'Elastic net CombiBMI':'tab:red'}\n",
    "analyte = 'proteins'\n",
    "title = 'Robustly retained proteins'\n",
    "\n",
    "#Prepare module sets\n",
    "print('Variables with non-zero beta-coefficient in all 10 models')\n",
    "tempD = {}\n",
    "for measure in tempD1.keys():\n",
    "    tempDF = tempD1[measure]\n",
    "    #Extract robust beta-coefficient: no zeros in all 10 models\n",
    "    tempDF1 = tempDF.loc[tempDF['nZeros']==0]\n",
    "    print(measure+':', len(tempDF1), 'per', len(tempDF), analyte+' (',\n",
    "          len(tempDF1)/len(tempDF)*100, '%)')\n",
    "    tempD[measure] = set(tempDF1.index.tolist())\n",
    "\n",
    "#Check common region (1,1,1,1)\n",
    "tempS = list(tempD.values())[0]#Initialize\n",
    "for measure in tempD.keys():\n",
    "    tempS = tempS & tempD[measure]\n",
    "print(' -> Common (1,1,1,1):', len(tempS))\n",
    "display(tempS)\n",
    "\n",
    "#Venn diagram\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "fig, ax = plt.subplots(figsize=(4, 4))\n",
    "venn(tempD, fmt='{size:,}', cmap=list(tempD2.values()), legend_loc=None, ax=ax)\n",
    "plt.setp(ax, ylim=(0.1, 0.875))#Otherwise, weird space...\n",
    "##Add legend annotation\n",
    "x_coord = [0.1, 0.1, 0.9, 0.9]\n",
    "y_coord = [0.25, 0.7, 0.7, 0.25]\n",
    "h_align = ['right', 'right', 'left', 'left']\n",
    "v_align = ['top', 'bottom', 'bottom', 'top']\n",
    "for i in range(len(tempD1)):\n",
    "    key = list(tempD2.keys())[i]\n",
    "    total = f'{len(tempD[key]):,}'\n",
    "    ax.text(x_coord[i], y_coord[i], key+'\\n('+total+' '+analyte+')',\n",
    "            fontsize='small', multialignment='center',\n",
    "            horizontalalignment=h_align[i], verticalalignment=v_align[i],\n",
    "            bbox={'boxstyle':'round', 'facecolor':tempD2[key], 'pad':0.2, 'alpha':0.5})\n",
    "ax.set_title(title, fontsize='medium')\n",
    "##Save\n",
    "fileDir = './ExportFigures/'\n",
    "ipynbName = '220901_Multiomics-BMI-NatMed1stRevision_bBMI-LASSO-vs-others_'\n",
    "fileName = 'LASSO-vs-ElasticNet-'+analyte+'.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": {
    "tags": []
   },
   "source": [
    "### 4-4. Clinical lab tests"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Prepare target analytes\n",
    "tempDF1 = bcoefDF_combi_LASSO.loc[bcoefDF_chem_LASSO.index.tolist()]\n",
    "tempDF2 = bcoefDF_combi_EN.loc[bcoefDF_chem_EN.index.tolist()]\n",
    "tempD1 = {'LASSO CombiBMI':tempDF1, 'LASSO ChemBMI':bcoefDF_chem_LASSO,\n",
    "          'Elastic net ChemBMI':bcoefDF_chem_EN, 'Elastic net CombiBMI':tempDF2}\n",
    "tempD2 = {'LASSO CombiBMI':'tab:blue', 'LASSO ChemBMI':'tab:orange',\n",
    "          'Elastic net ChemBMI':'tab:green', 'Elastic net CombiBMI':'tab:red'}\n",
    "analyte = 'tests'\n",
    "title = 'Robustly retained clinical lab tests'\n",
    "\n",
    "#Prepare module sets\n",
    "print('Variables with non-zero beta-coefficient in all 10 models')\n",
    "tempD = {}\n",
    "for measure in tempD1.keys():\n",
    "    tempDF = tempD1[measure]\n",
    "    #Extract robust beta-coefficient: no zeros in all 10 models\n",
    "    tempDF1 = tempDF.loc[tempDF['nZeros']==0]\n",
    "    print(measure+':', len(tempDF1), 'per', len(tempDF), analyte+' (',\n",
    "          len(tempDF1)/len(tempDF)*100, '%)')\n",
    "    tempD[measure] = set(tempDF1.index.tolist())\n",
    "\n",
    "#Check common region (1,1,1,1)\n",
    "tempS = list(tempD.values())[0]#Initialize\n",
    "for measure in tempD.keys():\n",
    "    tempS = tempS & tempD[measure]\n",
    "print(' -> Common (1,1,1,1):', len(tempS))\n",
    "display(tempS)\n",
    "\n",
    "#Venn diagram\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "fig, ax = plt.subplots(figsize=(4, 4))\n",
    "venn(tempD, fmt='{size:,}', cmap=list(tempD2.values()), legend_loc=None, ax=ax)\n",
    "plt.setp(ax, ylim=(0.1, 0.875))#Otherwise, weird space...\n",
    "##Add legend annotation\n",
    "x_coord = [0.1, 0.1, 0.9, 0.9]\n",
    "y_coord = [0.25, 0.7, 0.7, 0.25]\n",
    "h_align = ['right', 'right', 'left', 'left']\n",
    "v_align = ['top', 'bottom', 'bottom', 'top']\n",
    "for i in range(len(tempD1)):\n",
    "    key = list(tempD2.keys())[i]\n",
    "    total = f'{len(tempD[key]):,}'\n",
    "    ax.text(x_coord[i], y_coord[i], key+'\\n('+total+' '+analyte+')',\n",
    "            fontsize='small', multialignment='center',\n",
    "            horizontalalignment=h_align[i], verticalalignment=v_align[i],\n",
    "            bbox={'boxstyle':'round', 'facecolor':tempD2[key], 'pad':0.2, 'alpha':0.5})\n",
    "ax.set_title(title, fontsize='medium')\n",
    "##Save\n",
    "fileDir = './ExportFigures/'\n",
    "ipynbName = '220901_Multiomics-BMI-NatMed1stRevision_bBMI-LASSO-vs-others_'\n",
    "fileName = 'LASSO-vs-ElasticNet-'+analyte+'.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-5. Correlation of beta-coefficients"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#All variables\n",
    "tempD1 = {'LASSO':bcoefDF_met_LASSO, 'elastic net':bcoefDF_met_EN}\n",
    "tempD2 = {'LASSO':bcoefDF_prot_LASSO, 'elastic net':bcoefDF_prot_EN}\n",
    "tempD3 = {'LASSO':bcoefDF_chem_LASSO, 'elastic net':bcoefDF_chem_EN}\n",
    "tempD4 = {'LASSO':bcoefDF_combi_LASSO, 'elastic net':bcoefDF_combi_EN}\n",
    "tempD1 = {'MetBMI':tempD1, 'ProtBMI':tempD2, 'ChemBMI':tempD3, 'CombiBMI':tempD4}\n",
    "tempD2 = {'MetBMI':'Metabolomics', 'ProtBMI':'Proteomics',\n",
    "          'ChemBMI':'Clinical labs', 'CombiBMI':'Combined omics'}\n",
    "tempD3 = {'Metabolomics':'b', 'Proteomics':'r', 'Clinical labs':'g', 'Combined omics':'m'}\n",
    "yvar = 'LASSO'\n",
    "xvar = 'elastic net'\n",
    "xyvar_unit = '[log-scaled BMI (kg m'+r'$^{-2}$'+') per s.d.]'\n",
    "\n",
    "#Prepare DF\n",
    "tempD = {}\n",
    "for category in tempD1.keys():\n",
    "    tempD4 = tempD1[category]\n",
    "    tempDF1 = tempD4[xvar]\n",
    "    tempS1 = tempDF1['Mean'].copy()\n",
    "    tempS1.name = xvar+'_Bcoef'\n",
    "    tempDF2 = tempD4[yvar]\n",
    "    tempS2 = tempDF2['Mean'].copy()\n",
    "    tempS2.name = yvar+'_Bcoef'\n",
    "    tempDF = pd.merge(tempS1, tempS2, left_index=True, right_index=True, how='inner')\n",
    "    tempDF['Category'] = tempD2[category]\n",
    "    \n",
    "    tempD[category] = tempDF\n",
    "    \n",
    "    print(tempD2[category])\n",
    "    display(tempDF.describe(include='all'))\n",
    "tempDF = pd.concat(list(tempD.values()), axis=0)\n",
    "\n",
    "#Statistical tests\n",
    "tempDF1 = pd.DataFrame(columns=['N', 'DoF', 'Pearson_r', 'Pval'])\n",
    "for category in tempD3.keys():\n",
    "    tempDF2 = tempDF.loc[tempDF['Category']==category]\n",
    "    #Pearson's correlation\n",
    "    pearson_r, pval = stats.pearsonr(tempDF2[xvar+'_Bcoef'], tempDF2[yvar+'_Bcoef'])\n",
    "    size = len(tempDF2)\n",
    "    dof = size - 2\n",
    "    tempDF1.loc[category] = [size, dof, pearson_r, pval]\n",
    "##P-value adjustment by using Benjamini–Hochberg method\n",
    "tempDF1['AdjPval'] = multi.multipletests(tempDF1['Pval'], alpha=0.05, method='fdr_bh',\n",
    "                                         is_sorted=False, returnsorted=False)[1]\n",
    "tempDF1.index.rename('Category', inplace=True)\n",
    "tempDF1['N'] = tempDF1['N'].astype('int64')#Otherwise, float64!\n",
    "tempDF1['DoF'] = tempDF1['DoF'].astype('int64')#Otherwise, float64!\n",
    "display(tempDF1)\n",
    "##Save\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220901_Multiomics-BMI-NatMed1stRevision_bBMI-LASSO-vs-others_'\n",
    "fileName = 'LASSO-vs-ElasticNet-bcoefs_all.tsv'\n",
    "tempDF1.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "\n",
    "#Plot\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "fig, axes = plt.subplots(nrows=2, ncols=2,\n",
    "                         figsize=(3.5*2, 3.5*2-0.2), sharex=True, sharey=True)\n",
    "axis_xymin = -0.0425\n",
    "axis_xymax = 0.050\n",
    "xymin = -0.04\n",
    "xymax = 0.04\n",
    "xyinter = 0.02\n",
    "#Set axis range first; otherwise, regression line can be truncated differently\n",
    "plt.setp(axes, xlim=(axis_xymin, axis_xymax), xticks=np.arange(xymin, xymax+xyinter/10, xyinter))\n",
    "plt.setp(axes, ylim=(axis_xymin, axis_xymax), yticks=np.arange(xymin, xymax+xyinter/10, xyinter))\n",
    "for ax_i, ax in enumerate(axes.flat):\n",
    "    category = list(tempD3.keys())[ax_i]\n",
    "    #Prepare DF\n",
    "    tempDF2 = tempDF.loc[tempDF['Category']==category]\n",
    "    #Scatterplot with regression line\n",
    "    sns.regplot(data=tempDF2, x=xvar+'_Bcoef', y=yvar+'_Bcoef', color=tempD3[category],\n",
    "                scatter=True, fit_reg=True, ci=95, truncate=False, marker='o',\n",
    "                scatter_kws={'alpha':0.2, 'edgecolor':'k', 's':25}, ax=ax)\n",
    "    #Draw Y=X as reference\n",
    "    #ax.plot([axis_xymin, axis_xymax], [axis_xymin, axis_xymax],\n",
    "    #        color='black', linestyle=(0, (1, 2)), zorder=0)\n",
    "    #Draw Y=X=0 as reference\n",
    "    ax.axvline(x=0, color='black', linestyle=(0, (1, 2)), zorder=0)\n",
    "    ax.axhline(y=0, color='black', linestyle=(0, (1, 2)), zorder=0)\n",
    "    #Annotate Pearson's correlation\n",
    "    pearson_r = tempDF1['Pearson_r'].loc[category]\n",
    "    r_text = str(Decimal(str(pearson_r)).quantize(Decimal('0.001'), rounding=ROUND_HALF_UP))\n",
    "    pval = tempDF1['AdjPval'].loc[category]\n",
    "    below_limit = 0#Initialize\n",
    "    if pval==1.0:\n",
    "        pval_text = '1.0'\n",
    "    else:\n",
    "        if pval==0.0:#Due to smaller than the float minimum\n",
    "            pval = sys.float_info.min\n",
    "            print('P-value was smaller than the float minimum:', pval)\n",
    "            below_limit = 1\n",
    "        pval_text = f'{Decimal(str(pval)):.3E}'#Take more digits because rounding is bad here\n",
    "        significand, exponent = pval_text.split(sep='E-')\n",
    "        significand = str(Decimal(significand).quantize(Decimal('0.1'), rounding=ROUND_HALF_UP))\n",
    "        if significand=='10.0':\n",
    "            significand = '1.0'\n",
    "            exponent = str(int(exponent)-1)\n",
    "        if int(exponent)>2:\n",
    "            pval_text = significand+r'$\\times$'+'10'+r'$^{{-{0}}}$'.format(exponent)##Font is different in r'$ $'...\n",
    "        elif int(exponent)>0:\n",
    "            pval_text = '0.'+'0'*(int(exponent)-1)+significand.replace('.', '')\n",
    "        else:\n",
    "            pval_text = significand\n",
    "    if below_limit==1:\n",
    "        text = 'Pearson\\'s '+r'$r$'+' = '+r_text+'\\n'+r'$P$'+' < '+pval_text\n",
    "    else:\n",
    "        text = 'Pearson\\'s '+r'$r$'+' = '+r_text+'\\n'+r'$P$'+' = '+pval_text\n",
    "    ax.annotate(text, xy=(0.05, 0.95), xycoords='axes fraction',\n",
    "                horizontalalignment='left', verticalalignment='top',\n",
    "                multialignment='left', fontsize='small', color='k')\n",
    "    #Facet label\n",
    "    ax.set_title(category, {'fontsize':'large'})\n",
    "sns.despine()\n",
    "#Reset and generate common axis title\n",
    "plt.setp(axes, xlabel='', ylabel='')\n",
    "fig.tight_layout(pad=0.75)\n",
    "fig.text(x=0.54, y=0.02,#Manual adjustment\n",
    "         s='Mean of '+r'$\\beta$'+'-coefficients in '+xvar+' models\\n'+xyvar_unit,\n",
    "         fontsize='medium', verticalalignment='top', horizontalalignment='center')\n",
    "fig.text(x=0.0225, y=0.515,#Manual adjustment\n",
    "         s='Mean of '+r'$\\beta$'+'-coefficients in '+yvar+' models\\n'+xyvar_unit,\n",
    "         fontsize='medium', multialignment='center',\n",
    "         verticalalignment='center', horizontalalignment='right', rotation='vertical')\n",
    "##Save\n",
    "fileDir = './ExportFigures/'\n",
    "ipynbName = '220901_Multiomics-BMI-NatMed1stRevision_bBMI-LASSO-vs-others_'\n",
    "fileName = 'LASSO-vs-ElasticNet-bcoefs_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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Variables with non-zero beta-coefficients\n",
    "tempD1 = {'LASSO':bcoefDF_met_LASSO, 'elastic net':bcoefDF_met_EN}\n",
    "tempD2 = {'LASSO':bcoefDF_prot_LASSO, 'elastic net':bcoefDF_prot_EN}\n",
    "tempD3 = {'LASSO':bcoefDF_chem_LASSO, 'elastic net':bcoefDF_chem_EN}\n",
    "tempD4 = {'LASSO':bcoefDF_combi_LASSO, 'elastic net':bcoefDF_combi_EN}\n",
    "tempD1 = {'MetBMI':tempD1, 'ProtBMI':tempD2, 'ChemBMI':tempD3, 'CombiBMI':tempD4}\n",
    "tempD2 = {'MetBMI':'Metabolomics', 'ProtBMI':'Proteomics',\n",
    "          'ChemBMI':'Clinical labs', 'CombiBMI':'Combined omics'}\n",
    "tempD3 = {'Metabolomics':'b', 'Proteomics':'r', 'Clinical labs':'g', 'Combined omics':'m'}\n",
    "yvar = 'LASSO'\n",
    "xvar = 'elastic net'\n",
    "xyvar_unit = '[log-scaled BMI (kg m'+r'$^{-2}$'+') per s.d.]'\n",
    "\n",
    "#Prepare DF\n",
    "tempD = {}\n",
    "for category in tempD1.keys():\n",
    "    tempD4 = tempD1[category]\n",
    "    tempDF1 = tempD4[xvar]\n",
    "    tempS1 = tempDF1['Mean'].copy()\n",
    "    tempS1.name = xvar+'_Bcoef'\n",
    "    tempDF2 = tempD4[yvar]\n",
    "    tempS2 = tempDF2['Mean'].copy()\n",
    "    tempS2.name = yvar+'_Bcoef'\n",
    "    #Variables with non-zero beta-coefficients in either method\n",
    "    tempDF1 = tempDF1.loc[tempDF1['nZeros']!=10]\n",
    "    tempDF2 = tempDF2.loc[tempDF2['nZeros']!=10]\n",
    "    tempS = set(tempDF1.index.tolist()) | set(tempDF2.index.tolist())\n",
    "    tempDF = pd.merge(tempS1.loc[tempS], tempS2.loc[tempS],\n",
    "                      left_index=True, right_index=True, how='inner')\n",
    "    tempDF['Category'] = tempD2[category]\n",
    "    \n",
    "    tempD[category] = tempDF\n",
    "    \n",
    "    print(tempD2[category])\n",
    "    display(tempDF.describe(include='all'))\n",
    "tempDF = pd.concat(list(tempD.values()), axis=0)\n",
    "\n",
    "#Statistical tests\n",
    "tempDF1 = pd.DataFrame(columns=['N', 'DoF', 'Pearson_r', 'Pval'])\n",
    "for category in tempD3.keys():\n",
    "    tempDF2 = tempDF.loc[tempDF['Category']==category]\n",
    "    #Pearson's correlation\n",
    "    pearson_r, pval = stats.pearsonr(tempDF2[xvar+'_Bcoef'], tempDF2[yvar+'_Bcoef'])\n",
    "    size = len(tempDF2)\n",
    "    dof = size - 2\n",
    "    tempDF1.loc[category] = [size, dof, pearson_r, pval]\n",
    "##P-value adjustment by using Benjamini–Hochberg method\n",
    "tempDF1['AdjPval'] = multi.multipletests(tempDF1['Pval'], alpha=0.05, method='fdr_bh',\n",
    "                                         is_sorted=False, returnsorted=False)[1]\n",
    "tempDF1.index.rename('Category', inplace=True)\n",
    "tempDF1['N'] = tempDF1['N'].astype('int64')#Otherwise, float64!\n",
    "tempDF1['DoF'] = tempDF1['DoF'].astype('int64')#Otherwise, float64!\n",
    "display(tempDF1)\n",
    "##Save\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220901_Multiomics-BMI-NatMed1stRevision_bBMI-LASSO-vs-others_'\n",
    "fileName = 'LASSO-vs-ElasticNet-bcoefs_non-zero-in-any.tsv'\n",
    "tempDF1.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "\n",
    "#Plot\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "fig, axes = plt.subplots(nrows=2, ncols=2,\n",
    "                         figsize=(3.5*2, 3.5*2-0.2), sharex=True, sharey=True)\n",
    "axis_xymin = -0.0425\n",
    "axis_xymax = 0.050\n",
    "xymin = -0.04\n",
    "xymax = 0.04\n",
    "xyinter = 0.02\n",
    "#Set axis range first; otherwise, regression line can be truncated differently\n",
    "plt.setp(axes, xlim=(axis_xymin, axis_xymax), xticks=np.arange(xymin, xymax+xyinter/10, xyinter))\n",
    "plt.setp(axes, ylim=(axis_xymin, axis_xymax), yticks=np.arange(xymin, xymax+xyinter/10, xyinter))\n",
    "for ax_i, ax in enumerate(axes.flat):\n",
    "    category = list(tempD3.keys())[ax_i]\n",
    "    #Prepare DF\n",
    "    tempDF2 = tempDF.loc[tempDF['Category']==category]\n",
    "    #Scatterplot with regression line\n",
    "    sns.regplot(data=tempDF2, x=xvar+'_Bcoef', y=yvar+'_Bcoef', color=tempD3[category],\n",
    "                scatter=True, fit_reg=True, ci=95, truncate=False, marker='o',\n",
    "                scatter_kws={'alpha':0.2, 'edgecolor':'k', 's':25}, ax=ax)\n",
    "    #Draw Y=X as reference\n",
    "    #ax.plot([axis_xymin, axis_xymax], [axis_xymin, axis_xymax],\n",
    "    #        color='black', linestyle=(0, (1, 2)), zorder=0)\n",
    "    #Draw Y=X=0 as reference\n",
    "    ax.axvline(x=0, color='black', linestyle=(0, (1, 2)), zorder=0)\n",
    "    ax.axhline(y=0, color='black', linestyle=(0, (1, 2)), zorder=0)\n",
    "    #Annotate Pearson's correlation\n",
    "    pearson_r = tempDF1['Pearson_r'].loc[category]\n",
    "    r_text = str(Decimal(str(pearson_r)).quantize(Decimal('0.001'), rounding=ROUND_HALF_UP))\n",
    "    pval = tempDF1['AdjPval'].loc[category]\n",
    "    below_limit = 0#Initialize\n",
    "    if pval==1.0:\n",
    "        pval_text = '1.0'\n",
    "    else:\n",
    "        if pval==0.0:#Due to smaller than the float minimum\n",
    "            pval = sys.float_info.min\n",
    "            print('P-value was smaller than the float minimum:', pval)\n",
    "            below_limit = 1\n",
    "        pval_text = f'{Decimal(str(pval)):.3E}'#Take more digits because rounding is bad here\n",
    "        significand, exponent = pval_text.split(sep='E-')\n",
    "        significand = str(Decimal(significand).quantize(Decimal('0.1'), rounding=ROUND_HALF_UP))\n",
    "        if significand=='10.0':\n",
    "            significand = '1.0'\n",
    "            exponent = str(int(exponent)-1)\n",
    "        if int(exponent)>2:\n",
    "            pval_text = significand+r'$\\times$'+'10'+r'$^{{-{0}}}$'.format(exponent)##Font is different in r'$ $'...\n",
    "        elif int(exponent)>0:\n",
    "            pval_text = '0.'+'0'*(int(exponent)-1)+significand.replace('.', '')\n",
    "        else:\n",
    "            pval_text = significand\n",
    "    if below_limit==1:\n",
    "        text = 'Pearson\\'s '+r'$r$'+' = '+r_text+'\\n'+r'$P$'+' < '+pval_text\n",
    "    else:\n",
    "        text = 'Pearson\\'s '+r'$r$'+' = '+r_text+'\\n'+r'$P$'+' = '+pval_text\n",
    "    ax.annotate(text, xy=(0.05, 0.95), xycoords='axes fraction',\n",
    "                horizontalalignment='left', verticalalignment='top',\n",
    "                multialignment='left', fontsize='small', color='k')\n",
    "    #Facet label\n",
    "    ax.set_title(category, {'fontsize':'large'})\n",
    "sns.despine()\n",
    "#Reset and generate common axis title\n",
    "plt.setp(axes, xlabel='', ylabel='')\n",
    "fig.tight_layout(pad=0.75)\n",
    "fig.text(x=0.54, y=0.02,#Manual adjustment\n",
    "         s='Mean of '+r'$\\beta$'+'-coefficients in '+xvar+' models\\n'+xyvar_unit,\n",
    "         fontsize='medium', verticalalignment='top', horizontalalignment='center')\n",
    "fig.text(x=0.0225, y=0.515,#Manual adjustment\n",
    "         s='Mean of '+r'$\\beta$'+'-coefficients in '+yvar+' models\\n'+xyvar_unit,\n",
    "         fontsize='medium', multialignment='center',\n",
    "         verticalalignment='center', horizontalalignment='right', rotation='vertical')\n",
    "##Save\n",
    "fileDir = './ExportFigures/'\n",
    "ipynbName = '220901_Multiomics-BMI-NatMed1stRevision_bBMI-LASSO-vs-others_'\n",
    "fileName = 'LASSO-vs-ElasticNet-bcoefs_non-zero-in-any.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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Variables with non-zero beta-coefficient in all 10 models\n",
    "tempD1 = {'LASSO':bcoefDF_met_LASSO, 'elastic net':bcoefDF_met_EN}\n",
    "tempD2 = {'LASSO':bcoefDF_prot_LASSO, 'elastic net':bcoefDF_prot_EN}\n",
    "tempD3 = {'LASSO':bcoefDF_chem_LASSO, 'elastic net':bcoefDF_chem_EN}\n",
    "tempD4 = {'LASSO':bcoefDF_combi_LASSO, 'elastic net':bcoefDF_combi_EN}\n",
    "tempD1 = {'MetBMI':tempD1, 'ProtBMI':tempD2, 'ChemBMI':tempD3, 'CombiBMI':tempD4}\n",
    "tempD2 = {'MetBMI':'Metabolomics', 'ProtBMI':'Proteomics',\n",
    "          'ChemBMI':'Clinical labs', 'CombiBMI':'Combined omics'}\n",
    "tempD3 = {'Metabolomics':'b', 'Proteomics':'r', 'Clinical labs':'g', 'Combined omics':'m'}\n",
    "yvar = 'LASSO'\n",
    "xvar = 'elastic net'\n",
    "xyvar_unit = '[log-scaled BMI (kg m'+r'$^{-2}$'+') per s.d.]'\n",
    "\n",
    "#Prepare DF\n",
    "tempD = {}\n",
    "for category in tempD1.keys():\n",
    "    tempD4 = tempD1[category]\n",
    "    tempDF1 = tempD4[xvar]\n",
    "    tempS1 = tempDF1['Mean'].copy()\n",
    "    tempS1.name = xvar+'_Bcoef'\n",
    "    tempDF2 = tempD4[yvar]\n",
    "    tempS2 = tempDF2['Mean'].copy()\n",
    "    tempS2.name = yvar+'_Bcoef'\n",
    "    #Variables with non-zero beta-coefficients in bBMI or bWHtR\n",
    "    tempDF1 = tempDF1.loc[tempDF1['nZeros']==0]\n",
    "    tempDF2 = tempDF2.loc[tempDF2['nZeros']==0]\n",
    "    tempS = set(tempDF1.index.tolist()) | set(tempDF2.index.tolist())\n",
    "    tempDF = pd.merge(tempS1.loc[tempS], tempS2.loc[tempS],\n",
    "                      left_index=True, right_index=True, how='inner')\n",
    "    tempDF['Category'] = tempD2[category]\n",
    "    \n",
    "    tempD[category] = tempDF\n",
    "    \n",
    "    print(tempD2[category])\n",
    "    display(tempDF.describe(include='all'))\n",
    "tempDF = pd.concat(list(tempD.values()), axis=0)\n",
    "\n",
    "#Statistical tests\n",
    "tempDF1 = pd.DataFrame(columns=['N', 'DoF', 'Pearson_r', 'Pval'])\n",
    "for category in tempD3.keys():\n",
    "    tempDF2 = tempDF.loc[tempDF['Category']==category]\n",
    "    #Pearson's correlation\n",
    "    pearson_r, pval = stats.pearsonr(tempDF2[xvar+'_Bcoef'], tempDF2[yvar+'_Bcoef'])\n",
    "    size = len(tempDF2)\n",
    "    dof = size - 2\n",
    "    tempDF1.loc[category] = [size, dof, pearson_r, pval]\n",
    "##P-value adjustment by using Benjamini–Hochberg method\n",
    "tempDF1['AdjPval'] = multi.multipletests(tempDF1['Pval'], alpha=0.05, method='fdr_bh',\n",
    "                                         is_sorted=False, returnsorted=False)[1]\n",
    "tempDF1.index.rename('Category', inplace=True)\n",
    "tempDF1['N'] = tempDF1['N'].astype('int64')#Otherwise, float64!\n",
    "tempDF1['DoF'] = tempDF1['DoF'].astype('int64')#Otherwise, float64!\n",
    "display(tempDF1)\n",
    "##Save\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220901_Multiomics-BMI-NatMed1stRevision_bBMI-LASSO-vs-others_'\n",
    "fileName = 'LASSO-vs-ElasticNet-bcoefs_non-zero-in-all.tsv'\n",
    "tempDF1.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "\n",
    "#Plot\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "fig, axes = plt.subplots(nrows=2, ncols=2,\n",
    "                         figsize=(3.5*2, 3.5*2-0.2), sharex=True, sharey=True)\n",
    "axis_xymin = -0.0425\n",
    "axis_xymax = 0.050\n",
    "xymin = -0.04\n",
    "xymax = 0.04\n",
    "xyinter = 0.02\n",
    "#Set axis range first; otherwise, regression line can be truncated differently\n",
    "plt.setp(axes, xlim=(axis_xymin, axis_xymax), xticks=np.arange(xymin, xymax+xyinter/10, xyinter))\n",
    "plt.setp(axes, ylim=(axis_xymin, axis_xymax), yticks=np.arange(xymin, xymax+xyinter/10, xyinter))\n",
    "for ax_i, ax in enumerate(axes.flat):\n",
    "    category = list(tempD3.keys())[ax_i]\n",
    "    #Prepare DF\n",
    "    tempDF2 = tempDF.loc[tempDF['Category']==category]\n",
    "    #Scatterplot with regression line\n",
    "    sns.regplot(data=tempDF2, x=xvar+'_Bcoef', y=yvar+'_Bcoef', color=tempD3[category],\n",
    "                scatter=True, fit_reg=True, ci=95, truncate=False, marker='o',\n",
    "                scatter_kws={'alpha':0.2, 'edgecolor':'k', 's':25}, ax=ax)\n",
    "    #Draw Y=X as reference\n",
    "    #ax.plot([axis_xymin, axis_xymax], [axis_xymin, axis_xymax],\n",
    "    #        color='black', linestyle=(0, (1, 2)), zorder=0)\n",
    "    #Draw Y=X=0 as reference\n",
    "    ax.axvline(x=0, color='black', linestyle=(0, (1, 2)), zorder=0)\n",
    "    ax.axhline(y=0, color='black', linestyle=(0, (1, 2)), zorder=0)\n",
    "    #Annotate Pearson's correlation\n",
    "    pearson_r = tempDF1['Pearson_r'].loc[category]\n",
    "    r_text = str(Decimal(str(pearson_r)).quantize(Decimal('0.001'), rounding=ROUND_HALF_UP))\n",
    "    pval = tempDF1['AdjPval'].loc[category]\n",
    "    below_limit = 0#Initialize\n",
    "    if pval==1.0:\n",
    "        pval_text = '1.0'\n",
    "    else:\n",
    "        if pval==0.0:#Due to smaller than the float minimum\n",
    "            pval = sys.float_info.min\n",
    "            print('P-value was smaller than the float minimum:', pval)\n",
    "            below_limit = 1\n",
    "        pval_text = f'{Decimal(str(pval)):.3E}'#Take more digits because rounding is bad here\n",
    "        significand, exponent = pval_text.split(sep='E-')\n",
    "        significand = str(Decimal(significand).quantize(Decimal('0.1'), rounding=ROUND_HALF_UP))\n",
    "        if significand=='10.0':\n",
    "            significand = '1.0'\n",
    "            exponent = str(int(exponent)-1)\n",
    "        if int(exponent)>2:\n",
    "            pval_text = significand+r'$\\times$'+'10'+r'$^{{-{0}}}$'.format(exponent)##Font is different in r'$ $'...\n",
    "        elif int(exponent)>0:\n",
    "            pval_text = '0.'+'0'*(int(exponent)-1)+significand.replace('.', '')\n",
    "        else:\n",
    "            pval_text = significand\n",
    "    if below_limit==1:\n",
    "        text = 'Pearson\\'s '+r'$r$'+' = '+r_text+'\\n'+r'$P$'+' < '+pval_text\n",
    "    else:\n",
    "        text = 'Pearson\\'s '+r'$r$'+' = '+r_text+'\\n'+r'$P$'+' = '+pval_text\n",
    "    ax.annotate(text, xy=(0.05, 0.95), xycoords='axes fraction',\n",
    "                horizontalalignment='left', verticalalignment='top',\n",
    "                multialignment='left', fontsize='small', color='k')\n",
    "    #Facet label\n",
    "    ax.set_title(category, {'fontsize':'large'})\n",
    "sns.despine()\n",
    "#Reset and generate common axis title\n",
    "plt.setp(axes, xlabel='', ylabel='')\n",
    "fig.tight_layout(pad=0.75)\n",
    "fig.text(x=0.54, y=0.02,#Manual adjustment\n",
    "         s='Mean of '+r'$\\beta$'+'-coefficients in '+xvar+' models\\n'+xyvar_unit,\n",
    "         fontsize='medium', verticalalignment='top', horizontalalignment='center')\n",
    "fig.text(x=0.0225, y=0.515,#Manual adjustment\n",
    "         s='Mean of '+r'$\\beta$'+'-coefficients in '+yvar+' models\\n'+xyvar_unit,\n",
    "         fontsize='medium', multialignment='center',\n",
    "         verticalalignment='center', horizontalalignment='right', rotation='vertical')\n",
    "##Save\n",
    "fileDir = './ExportFigures/'\n",
    "ipynbName = '220901_Multiomics-BMI-NatMed1stRevision_bBMI-LASSO-vs-others_'\n",
    "fileName = 'LASSO-vs-ElasticNet-bcoefs_non-zero-in-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": [
    "## 5. Top 30 variables in ridge models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Because beta-coefficients can be adjusted among the high-collinear variables and thus almost all variables are retained across 10 models in Ridge, it is not so simple to compare the beta-coefficients between LASSO and ridge models.  \n",
    "> –> Hence, only the top 30 variables that had the highest absolute value of the mean beta-coefficient are checked.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5-1. Import beta-coefficients"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Ridge\n",
    "print('Ridge models')\n",
    "tempD = {}\n",
    "for bbmi in ['MetBMI', 'ProtBMI', 'ChemBMI', 'CombiBMI']:\n",
    "    #Import the beta-coefficients\n",
    "    fileDir = './ExportData/'\n",
    "    ipynbName = '220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-Ridge_'\n",
    "    fileName = bbmi+'-BothSex-Ridgebcoefs.tsv'\n",
    "    tempDF = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t').set_index('Variable')\n",
    "    tempDF = tempDF.drop(index=['Intercept'])\n",
    "    tempD[bbmi] = tempDF\n",
    "    \n",
    "    #Check\n",
    "    print(bbmi+':')\n",
    "    print(' - Variables:', len(tempDF))\n",
    "    #Variables with non-zero beta-coefficient\n",
    "    tempDF1 = tempDF.loc[tempDF['nZeros']!=10]\n",
    "    print(' - Variables with non-zero beta-coefficient:', len(tempDF1),\n",
    "          '(', len(tempDF1)/len(tempDF)*100, '%)')\n",
    "    #Extract robust beta-coefficient: no zeros in all 10 models\n",
    "    tempDF1 = tempDF.loc[tempDF['nZeros']==0]\n",
    "    tempDF1 = tempDF1.sort_values(by='Mean', ascending=False)\n",
    "    print(' - Variables with non-zero beta-coefficient in all 10 models:', len(tempDF1),\n",
    "          '(', len(tempDF1)/len(tempDF)*100, '%)')\n",
    "    print('')\n",
    "\n",
    "bcoefDF_met_ridge = tempD['MetBMI']\n",
    "bcoefDF_prot_ridge = tempD['ProtBMI']\n",
    "bcoefDF_chem_ridge = tempD['ChemBMI']\n",
    "bcoefDF_combi_ridge = tempD['CombiBMI']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5-2. Metabolomics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF = bcoefDF_met_ridge.copy()\n",
    "topX = 30\n",
    "method = 'Ridge'\n",
    "bbmi_color = 'b'\n",
    "bbmi = 'MetBMI'\n",
    "\n",
    "#Top X variables based on the absolute value of the mean beta-coefficient\n",
    "tempDF['AbsMean'] = np.abs(tempDF['Mean'])\n",
    "tempDF = tempDF.sort_values(by='AbsMean', ascending=False)\n",
    "tempDF = tempDF.iloc[:topX]\n",
    "print('Visualize top', len(tempDF), 'variables')\n",
    "\n",
    "#Prepare DF for plot\n",
    "tempDF = tempDF.sort_values(by='Mean', ascending=False)#Re-sort\n",
    "tempDF = tempDF.drop(columns=['Mean', 'SD', 'nZeros', 'AbsMean'])\n",
    "tempDF = tempDF.reset_index().melt(var_name='Model', value_name='bcoef', id_vars=['Variable'])\n",
    "\n",
    "#Plot\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "plt.figure(figsize=(5, 10))\n",
    "p = sns.boxplot(data=tempDF, y='Variable', x='bcoef', color=bbmi_color, dodge=False, saturation=1,\n",
    "                showfliers=True, flierprops={'marker':'o', 'markerfacecolor':'gray', 'alpha':0.4},\n",
    "                showcaps=True, notch=False)\n",
    "p.set(xlim=(-0.06, 0.09), xticks=np.arange(-0.04, 0.081, 0.04))#Fixed across omics\n",
    "p.grid(axis='x', linestyle='--', color='black')\n",
    "sns.despine()\n",
    "##Change default dull line color of sns.boxplot (saturation parameter is for patch)\n",
    "for line in p.get_lines():\n",
    "    line.set_color('k')\n",
    "for box in p.artists:\n",
    "    box.set_edgecolor('k')\n",
    "##Add background color\n",
    "for row_i in range(topX):\n",
    "    if row_i%2 == 0:\n",
    "        plt.axhspan(ymin=row_i-0.5, ymax=row_i+0.5, facecolor=bbmi_color, alpha=0.2, zorder=0)\n",
    "plt.ylabel('')\n",
    "plt.xlabel(r'$\\beta$'+'-coefficient in '+method.lower()+' model\\n [log-scaled BMI (kg m'+r'$^{-2}$'+') per s.d.]')\n",
    "##Save\n",
    "fileDir = './ExportFigures/'\n",
    "ipynbName = '220901_Multiomics-BMI-NatMed1stRevision_bBMI-LASSO-vs-others_'\n",
    "fileName = method+'-'+bbmi+'-bcoef_top'+str(topX)+'.tif'\n",
    "plt.gcf().savefig(fileDir+ipynbName+fileName, dpi=300, bbox_inches='tight', pad_inches=0.04,\n",
    "                  pil_kwargs={'compression':'tiff_lzw'})\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5-3. Proteomics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF = bcoefDF_prot_ridge.copy()\n",
    "topX = 30\n",
    "method = 'Ridge'\n",
    "bbmi_color = 'r'\n",
    "bbmi = 'ProtBMI'\n",
    "\n",
    "#Top X variables based on the absolute value of the mean beta-coefficient\n",
    "tempDF['AbsMean'] = np.abs(tempDF['Mean'])\n",
    "tempDF = tempDF.sort_values(by='AbsMean', ascending=False)\n",
    "tempDF = tempDF.iloc[:topX]\n",
    "print('Visualize top', len(tempDF), 'variables')\n",
    "\n",
    "#Prepare DF for plot\n",
    "tempDF = tempDF.sort_values(by='Mean', ascending=False)#Re-sort\n",
    "tempDF = tempDF.drop(columns=['Mean', 'SD', 'nZeros', 'AbsMean'])\n",
    "tempDF = tempDF.reset_index().melt(var_name='Model', value_name='bcoef', id_vars=['Variable'])\n",
    "\n",
    "#Plot\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "plt.figure(figsize=(5, 10))\n",
    "p = sns.boxplot(data=tempDF, y='Variable', x='bcoef', color=bbmi_color, dodge=False, saturation=1,\n",
    "                showfliers=True, flierprops={'marker':'o', 'markerfacecolor':'gray', 'alpha':0.4},\n",
    "                showcaps=True, notch=False)\n",
    "p.set(xlim=(-0.06, 0.09), xticks=np.arange(-0.04, 0.081, 0.04))#Fixed across omics\n",
    "p.grid(axis='x', linestyle='--', color='black')\n",
    "sns.despine()\n",
    "##Change default dull line color of sns.boxplot (saturation parameter is for patch)\n",
    "for line in p.get_lines():\n",
    "    line.set_color('k')\n",
    "for box in p.artists:\n",
    "    box.set_edgecolor('k')\n",
    "##Add background color\n",
    "for row_i in range(topX):\n",
    "    if row_i%2 == 0:\n",
    "        plt.axhspan(ymin=row_i-0.5, ymax=row_i+0.5, facecolor=bbmi_color, alpha=0.2, zorder=0)\n",
    "plt.ylabel('')\n",
    "plt.xlabel(r'$\\beta$'+'-coefficient in '+method.lower()+' model\\n [log-scaled BMI (kg m'+r'$^{-2}$'+') per s.d.]')\n",
    "##Save\n",
    "fileDir = './ExportFigures/'\n",
    "ipynbName = '220901_Multiomics-BMI-NatMed1stRevision_bBMI-LASSO-vs-others_'\n",
    "fileName = method+'-'+bbmi+'-bcoef_top'+str(topX)+'.tif'\n",
    "plt.gcf().savefig(fileDir+ipynbName+fileName, dpi=300, bbox_inches='tight', pad_inches=0.04,\n",
    "                  pil_kwargs={'compression':'tiff_lzw'})\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5-4. Clinical labs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF = bcoefDF_chem_ridge.copy()\n",
    "topX = 30\n",
    "method = 'Ridge'\n",
    "bbmi_color = 'g'\n",
    "bbmi = 'ChemBMI'\n",
    "\n",
    "#Top X variables based on the absolute value of the mean beta-coefficient\n",
    "tempDF['AbsMean'] = np.abs(tempDF['Mean'])\n",
    "tempDF = tempDF.sort_values(by='AbsMean', ascending=False)\n",
    "tempDF = tempDF.iloc[:topX]\n",
    "print('Visualize top', len(tempDF), 'variables')\n",
    "\n",
    "#Prepare DF for plot\n",
    "tempDF = tempDF.sort_values(by='Mean', ascending=False)#Re-sort\n",
    "tempDF = tempDF.drop(columns=['Mean', 'SD', 'nZeros', 'AbsMean'])\n",
    "tempDF = tempDF.reset_index().melt(var_name='Model', value_name='bcoef', id_vars=['Variable'])\n",
    "\n",
    "#Plot\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "plt.figure(figsize=(5, 10))\n",
    "p = sns.boxplot(data=tempDF, y='Variable', x='bcoef', color=bbmi_color, dodge=False, saturation=1,\n",
    "                showfliers=True, flierprops={'marker':'o', 'markerfacecolor':'gray', 'alpha':0.4},\n",
    "                showcaps=True, notch=False)\n",
    "p.set(xlim=(-0.06, 0.09), xticks=np.arange(-0.04, 0.081, 0.04))#Fixed across omics\n",
    "p.grid(axis='x', linestyle='--', color='black')\n",
    "sns.despine()\n",
    "##Change default dull line color of sns.boxplot (saturation parameter is for patch)\n",
    "for line in p.get_lines():\n",
    "    line.set_color('k')\n",
    "for box in p.artists:\n",
    "    box.set_edgecolor('k')\n",
    "##Add background color\n",
    "for row_i in range(topX):\n",
    "    if row_i%2 == 0:\n",
    "        plt.axhspan(ymin=row_i-0.5, ymax=row_i+0.5, facecolor=bbmi_color, alpha=0.2, zorder=0)\n",
    "plt.ylabel('')\n",
    "plt.xlabel(r'$\\beta$'+'-coefficient in '+method.lower()+' model\\n [log-scaled BMI (kg m'+r'$^{-2}$'+') per s.d.]')\n",
    "##Save\n",
    "fileDir = './ExportFigures/'\n",
    "ipynbName = '220901_Multiomics-BMI-NatMed1stRevision_bBMI-LASSO-vs-others_'\n",
    "fileName = method+'-'+bbmi+'-bcoef_top'+str(topX)+'.tif'\n",
    "plt.gcf().savefig(fileDir+ipynbName+fileName, dpi=300, bbox_inches='tight', pad_inches=0.04,\n",
    "                  pil_kwargs={'compression':'tiff_lzw'})\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5-5. Combined omics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF = bcoefDF_combi_ridge.copy()\n",
    "topX = 30\n",
    "method = 'Ridge'\n",
    "#bbmi_color = 'm'\n",
    "tempD1 = {'Metabolomics':bcoefDF_met_ridge.index.tolist(),\n",
    "          'Proteomics':bcoefDF_prot_ridge.index.tolist(),\n",
    "          'Clinical labs':bcoefDF_chem_ridge.index.tolist()}\n",
    "tempD2 = {'Metabolomics':'b', 'Proteomics':'r', 'Clinical labs':'g'}\n",
    "bbmi = 'CombiBMI'\n",
    "\n",
    "#Top X variables based on the absolute value of the mean beta-coefficient\n",
    "tempDF['AbsMean'] = np.abs(tempDF['Mean'])\n",
    "tempDF = tempDF.sort_values(by='AbsMean', ascending=False)\n",
    "tempDF = tempDF.iloc[:topX]\n",
    "print('Visualize top', len(tempDF), 'variables')\n",
    "\n",
    "#Prepare DF for plot\n",
    "tempDF = tempDF.sort_values(by='Mean', ascending=False)#Re-sort\n",
    "tempDF = tempDF.drop(columns=['Mean', 'SD', 'nZeros', 'AbsMean'])\n",
    "tempDF1 = tempDF.reset_index().melt(var_name='Model', value_name='bcoef', id_vars=['Variable'])\n",
    "\n",
    "#Prepare category color\n",
    "tempDF2 = pd.DataFrame(index=tempDF.index)\n",
    "tempL1 = []\n",
    "tempL2 = []\n",
    "for row_n in tempDF2.index.tolist():\n",
    "    count = 0#Initialize\n",
    "    for category in tempD1.keys():\n",
    "        if row_n in tempD1[category]:\n",
    "            tempL1.append(category)\n",
    "            tempL2.append(tempD2[category])\n",
    "            count += 1\n",
    "    if count==0:\n",
    "        print('Check error!')\n",
    "tempDF2['Category'] = tempL1\n",
    "tempDF2['Color'] = tempL2\n",
    "\n",
    "#Plot\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "plt.figure(figsize=(5, 10))\n",
    "p = sns.boxplot(data=tempDF1, y='Variable', order=tempDF2.index.tolist(), x='bcoef',\n",
    "                palette=tempDF2['Color'], dodge=False, saturation=1,\n",
    "                showfliers=True, flierprops={'marker':'o', 'markerfacecolor':'gray', 'alpha':0.4},\n",
    "                showcaps=True, notch=False)\n",
    "#p.set(xlim=(-0.06, 0.09), xticks=np.arange(-0.04, 0.081, 0.04))#Fixed across omics\n",
    "p.set(xlim=(-0.021, 0.021), xticks=np.arange(-0.02, 0.021, 0.01))#Based on the LASSO figure\n",
    "p.grid(axis='x', linestyle='--', color='black')\n",
    "sns.despine()\n",
    "##Change default dull line color of sns.boxplot (saturation parameter is for patch)\n",
    "for line in p.get_lines():\n",
    "    line.set_color('k')\n",
    "for box in p.artists:\n",
    "    box.set_edgecolor('k')\n",
    "##Add background color\n",
    "for row_i in range(len(tempDF2)):\n",
    "    cat_color = tempDF2['Color'].iloc[row_i]\n",
    "    if row_i%2 == 0:\n",
    "        plt.axhspan(ymin=row_i-0.5, ymax=row_i+0.5, facecolor=cat_color, alpha=0.4, zorder=0)\n",
    "    else:\n",
    "        plt.axhspan(ymin=row_i-0.5, ymax=row_i+0.5, facecolor=cat_color, alpha=0.4, zorder=0)\n",
    "plt.ylabel('')\n",
    "plt.xlabel(r'$\\beta$'+'-coefficient in '+method.lower()+' model\\n [log-scaled BMI (kg m'+r'$^{-2}$'+') per s.d.]')\n",
    "#Add legend\n",
    "tempL = []\n",
    "for category in tempD1.keys():\n",
    "    tempL.append(mpatches.Patch(facecolor=tempD2[category], edgecolor='k', label=category))\n",
    "plt.legend(handles=tempL, fontsize='large',\n",
    "           title='Omics category', title_fontsize='x-large',\n",
    "           bbox_to_anchor=(-0.35, 0), loc='center right', borderaxespad=0)#Manual adjustment\n",
    "##Save\n",
    "fileDir = './ExportFigures/'\n",
    "ipynbName = '220901_Multiomics-BMI-NatMed1stRevision_bBMI-LASSO-vs-others_'\n",
    "fileName = method+'-'+bbmi+'-bcoef_top'+str(topX)+'.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": [
    "## 6. Top 30 variables in RF models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Because RF model is not a linear model, it is not so simple to compare the predictor variables between LASSO and RF models.  \n",
    "> –> Hence, only the top 30 variables that had the highest feature importance are checked.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6-1. Import feature importance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> The feature importance that were obtained from sklearn.ensemble.RandomForestRegressor is based on impurity; i.e., the importance of a feature is computed as the (normalized) total reduction of the criterion (the mean squared error in this case) brought by that feature.  \n",
    "> –> Hence, the sum of feature importances in all variables within the model becomes 1.  \n",
    ">\n",
    "> Of note, it is known that this impurity-based feature importance can be misleading for high cardinality features, and alternative index such as permutation importance and SHAP (SHapley Additive exPlanations) feature importance would be better. In this study, however, the impurity-based feature importance is used because all variables are continuous (and because of a supplementary result).  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Random forest\n",
    "print('Random forest models')\n",
    "tempD = {}\n",
    "for bbmi in ['MetBMI', 'ProtBMI', 'ChemBMI', 'CombiBMI']:\n",
    "    #Import the feature importance\n",
    "    fileDir = './ImportData/'\n",
    "    ipynbName = '220828_Multiomics-BMI-NatMed1stRevision_BMI-baseline-RF-wenceslaus_'\n",
    "    fileName = bbmi+'-BothSex-feature-importance.tsv'\n",
    "    tempDF = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t').set_index('Variable')\n",
    "    tempD[bbmi] = tempDF\n",
    "    \n",
    "    #Check\n",
    "    print(bbmi+':')\n",
    "    print(' - Variables:', len(tempDF))\n",
    "    tempDF1 = tempDF.loc[tempDF['Mean']>0.01]\n",
    "    print(' - Variables with the mean of feature importances > 0.01:', len(tempDF1),\n",
    "          '(', len(tempDF1)/len(tempDF)*100, '%)')\n",
    "    tempDF1 = tempDF.loc[tempDF['Mean']>0.05]\n",
    "    print(' - Variables with the mean of feature importances > 0.05:', len(tempDF1),\n",
    "          '(', len(tempDF1)/len(tempDF)*100, '%)')\n",
    "    tempDF1 = tempDF.loc[:, tempDF.columns.str.contains('Model_')]\n",
    "    display(tempDF1.sum(axis=0))\n",
    "    print('')\n",
    "\n",
    "bcoefDF_met_RF = tempD['MetBMI']\n",
    "bcoefDF_prot_RF = tempD['ProtBMI']\n",
    "bcoefDF_chem_RF = tempD['ChemBMI']\n",
    "bcoefDF_combi_RF = tempD['CombiBMI']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6-2. Metabolomics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF = bcoefDF_met_RF.copy()\n",
    "topX = 30\n",
    "method = 'RF'\n",
    "bbmi_color = 'b'\n",
    "bbmi = 'MetBMI'\n",
    "\n",
    "#Top X variables based on feature importance\n",
    "tempDF = tempDF.sort_values(by='Mean', ascending=False)\n",
    "tempDF = tempDF.iloc[:topX]\n",
    "print('Visualize top', len(tempDF), 'variables')\n",
    "\n",
    "#Prepare DF for plot\n",
    "tempDF = tempDF.drop(columns=['Mean', 'SD'])\n",
    "tempDF = tempDF * 100#Convert to percentage\n",
    "tempDF = tempDF.reset_index().melt(var_name='Model', value_name='importance', id_vars=['Variable'])\n",
    "\n",
    "#Plot\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "plt.figure(figsize=(5, 10))\n",
    "sns.barplot(data=tempDF, y='Variable', x='importance',\n",
    "            color=bbmi_color, dodge=False, saturation=1,\n",
    "            ci=95, capsize=0.4, errwidth=1.5, errcolor='black', edgecolor='black')\n",
    "p = sns.stripplot(data=tempDF, y='Variable', x='importance',\n",
    "                  dodge=False, jitter=0.3,\n",
    "                  size=5, edgecolor='black', color='gray', linewidth=1, alpha=0.4)\n",
    "p.set(xlim=(0, 32.5), xticks=np.arange(0, 30.1, 10))#Fixed across omics\n",
    "p.grid(axis='x', linestyle='--', color='black')\n",
    "sns.despine()\n",
    "plt.ylabel('')\n",
    "plt.xlabel('Feature importance in '+method+' model [%]')\n",
    "##Save\n",
    "fileDir = './ExportFigures/'\n",
    "ipynbName = '220901_Multiomics-BMI-NatMed1stRevision_bBMI-LASSO-vs-others_'\n",
    "fileName = method+'-'+bbmi+'-importance_top'+str(topX)+'.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": [
    "### 6-3. Proteomics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF = bcoefDF_prot_RF.copy()\n",
    "topX = 30\n",
    "method = 'RF'\n",
    "bbmi_color = 'r'\n",
    "bbmi = 'ProtBMI'\n",
    "\n",
    "#Top X variables based on feature importance\n",
    "tempDF = tempDF.sort_values(by='Mean', ascending=False)\n",
    "tempDF = tempDF.iloc[:topX]\n",
    "print('Visualize top', len(tempDF), 'variables')\n",
    "\n",
    "#Prepare DF for plot\n",
    "tempDF = tempDF.drop(columns=['Mean', 'SD'])\n",
    "tempDF = tempDF * 100#Convert to percentage\n",
    "tempDF = tempDF.reset_index().melt(var_name='Model', value_name='importance', id_vars=['Variable'])\n",
    "\n",
    "#Plot\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "plt.figure(figsize=(5, 10))\n",
    "sns.barplot(data=tempDF, y='Variable', x='importance',\n",
    "            color=bbmi_color, dodge=False, saturation=1,\n",
    "            ci=95, capsize=0.4, errwidth=1.5, errcolor='black', edgecolor='black')\n",
    "p = sns.stripplot(data=tempDF, y='Variable', x='importance',\n",
    "                  dodge=False, jitter=0.3,\n",
    "                  size=5, edgecolor='black', color='gray', linewidth=1, alpha=0.4)\n",
    "p.set(xlim=(0, 32.5), xticks=np.arange(0, 30.1, 10))#Fixed across omics\n",
    "p.grid(axis='x', linestyle='--', color='black')\n",
    "sns.despine()\n",
    "plt.ylabel('')\n",
    "plt.xlabel('Feature importance in '+method+' model [%]')\n",
    "##Save\n",
    "fileDir = './ExportFigures/'\n",
    "ipynbName = '220901_Multiomics-BMI-NatMed1stRevision_bBMI-LASSO-vs-others_'\n",
    "fileName = method+'-'+bbmi+'-importance_top'+str(topX)+'.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": [
    "### 6-4. Clinical labs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF = bcoefDF_chem_RF.copy()\n",
    "topX = 30\n",
    "method = 'RF'\n",
    "bbmi_color = 'g'\n",
    "bbmi = 'ChemBMI'\n",
    "\n",
    "#Top X variables based on feature importance\n",
    "tempDF = tempDF.sort_values(by='Mean', ascending=False)\n",
    "tempDF = tempDF.iloc[:topX]\n",
    "print('Visualize top', len(tempDF), 'variables')\n",
    "\n",
    "#Prepare DF for plot\n",
    "tempDF = tempDF.drop(columns=['Mean', 'SD'])\n",
    "tempDF = tempDF * 100#Convert to percentage\n",
    "tempDF = tempDF.reset_index().melt(var_name='Model', value_name='importance', id_vars=['Variable'])\n",
    "\n",
    "#Plot\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "plt.figure(figsize=(5, 10))\n",
    "sns.barplot(data=tempDF, y='Variable', x='importance',\n",
    "            color=bbmi_color, dodge=False, saturation=1,\n",
    "            ci=95, capsize=0.4, errwidth=1.5, errcolor='black', edgecolor='black')\n",
    "p = sns.stripplot(data=tempDF, y='Variable', x='importance',\n",
    "                  dodge=False, jitter=0.3,\n",
    "                  size=5, edgecolor='black', color='gray', linewidth=1, alpha=0.4)\n",
    "p.set(xlim=(0, 32.5), xticks=np.arange(0, 30.1, 10))#Fixed across omics\n",
    "p.grid(axis='x', linestyle='--', color='black')\n",
    "sns.despine()\n",
    "plt.ylabel('')\n",
    "plt.xlabel('Feature importance in '+method+' model [%]')\n",
    "##Save\n",
    "fileDir = './ExportFigures/'\n",
    "ipynbName = '220901_Multiomics-BMI-NatMed1stRevision_bBMI-LASSO-vs-others_'\n",
    "fileName = method+'-'+bbmi+'-importance_top'+str(topX)+'.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": [
    "### 6-5. Combined omics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF = bcoefDF_combi_RF.copy()\n",
    "topX = 30\n",
    "method = 'RF'\n",
    "#bbmi_color = 'm'\n",
    "tempD1 = {'Metabolomics':bcoefDF_met_RF.index.tolist(),\n",
    "          'Proteomics':bcoefDF_prot_RF.index.tolist(),\n",
    "          'Clinical labs':bcoefDF_chem_RF.index.tolist()}\n",
    "tempD2 = {'Metabolomics':'b', 'Proteomics':'r', 'Clinical labs':'g'}\n",
    "bbmi = 'CombiBMI'\n",
    "\n",
    "#Top X variables based on feature importance\n",
    "tempDF = tempDF.sort_values(by='Mean', ascending=False)\n",
    "tempDF = tempDF.iloc[:topX]\n",
    "print('Visualize top', len(tempDF), 'variables')\n",
    "\n",
    "#Prepare DF for plot\n",
    "tempDF = tempDF.drop(columns=['Mean', 'SD'])\n",
    "tempDF = tempDF * 100#Convert to percentage\n",
    "tempDF1 = tempDF.reset_index().melt(var_name='Model', value_name='importance', id_vars=['Variable'])\n",
    "\n",
    "#Prepare category color\n",
    "tempDF2 = pd.DataFrame(index=tempDF.index)\n",
    "tempL1 = []\n",
    "tempL2 = []\n",
    "for row_n in tempDF2.index.tolist():\n",
    "    count = 0#Initialize\n",
    "    for category in tempD1.keys():\n",
    "        if row_n in tempD1[category]:\n",
    "            tempL1.append(category)\n",
    "            tempL2.append(tempD2[category])\n",
    "            count += 1\n",
    "    if count==0:\n",
    "        print('Check error!')\n",
    "tempDF2['Category'] = tempL1\n",
    "tempDF2['Color'] = tempL2\n",
    "\n",
    "#Plot\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "plt.figure(figsize=(5, 10))\n",
    "sns.barplot(data=tempDF1, y='Variable', order=tempDF2.index.tolist(), x='importance',\n",
    "            palette=tempDF2['Color'], dodge=False, saturation=1,\n",
    "            ci=95, capsize=0.4, errwidth=1.5, errcolor='black', edgecolor='black')\n",
    "p = sns.stripplot(data=tempDF1, y='Variable', order=tempDF2.index.tolist(), x='importance',\n",
    "                  dodge=False, jitter=0.3,\n",
    "                  size=5, edgecolor='black', color='gray', linewidth=1, alpha=0.4)\n",
    "#p.set(xlim=(0, 32.5), xticks=np.arange(0, 30.1, 10))#Fixed across omics\n",
    "p.grid(axis='x', linestyle='--', color='black')\n",
    "sns.despine()\n",
    "##Add background color\n",
    "for row_i in range(len(tempDF2)):\n",
    "    cat_color = tempDF2['Color'].iloc[row_i]\n",
    "    if row_i%2 == 0:\n",
    "        plt.axhspan(ymin=row_i-0.5, ymax=row_i+0.5, facecolor=cat_color, alpha=0.4, zorder=0)\n",
    "    else:\n",
    "        plt.axhspan(ymin=row_i-0.5, ymax=row_i+0.5, facecolor=cat_color, alpha=0.4, zorder=0)\n",
    "plt.margins(y=0)\n",
    "plt.ylabel('')\n",
    "plt.xlabel('Feature importance in '+method+' model [%]')\n",
    "#Add legend\n",
    "tempL = []\n",
    "for category in tempD1.keys():\n",
    "    tempL.append(mpatches.Patch(facecolor=tempD2[category], edgecolor='k', label=category))\n",
    "plt.legend(handles=tempL, fontsize='large',\n",
    "           title='Omics category', title_fontsize='x-large',\n",
    "           bbox_to_anchor=(1, 0), loc='lower right', borderaxespad=0.25)#Manual adjustment\n",
    "##Save\n",
    "fileDir = './ExportFigures/'\n",
    "ipynbName = '220901_Multiomics-BMI-NatMed1stRevision_bBMI-LASSO-vs-others_'\n",
    "fileName = method+'-'+bbmi+'-importance_top'+str(topX)+'.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"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}