Code05_BMI-Ridge_Arivale-baseline.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multiomics BMI Paper — Ridge-modeling for BMI from the Arivale Baseline Omics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***by Kengo Watanabe***  \n",
    "\n",
    "This Jupyter Notebook (with Python 3 kernel) generated the ridge linear regression models for predicting BMI (biological BMI) from each of the Arivale baseline blood omic datasets, and calculated the testing (hold-out) set-derived BMI predictions for the Arivale cohort.  \n",
    "\n",
    "Input files:  \n",
    "* Arivale baseline BMI and blood omics (preprocessed): 210104_Biological-BMI-paper_RF-imputation_baseline-\\[metDF/protDF/chemDF/combiDF\\]-with-RF-imputation.tsv  \n",
    "\n",
    "Output figures and tables:  \n",
    "* Intermediate tables for other notebooks (BMI predictions, beta-coefficients)  \n",
    "\n",
    "Original notebook (memo for my future tracing):  \n",
    "* dalek:\\[JupyterLab HOME\\]/220621_Multiomics-BMI-NatMedRevision/220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-Ridge.ipynb  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# packages in environment at /opt/conda/envs/arivale-py3:\n",
      "#\n",
      "# Name                    Version                   Build  Channel\n",
      "_libgcc_mutex             0.1                 conda_forge    conda-forge\n",
      "_openmp_mutex             4.5                       1_gnu    conda-forge\n",
      "analytics                 0.1                      pypi_0    pypi\n",
      "argon2-cffi               21.1.0           py39h3811e60_0    conda-forge\n",
      "arivale-data-interface    0.1.0                    pypi_0    pypi\n",
      "async_generator           1.10                       py_0    conda-forge\n",
      "atk-1.0                   2.36.0               h3371d22_4    conda-forge\n",
      "attrs                     21.2.0             pyhd8ed1ab_0    conda-forge\n",
      "backcall                  0.2.0              pyh9f0ad1d_0    conda-forge\n",
      "backports                 1.0                        py_2    conda-forge\n",
      "backports.functools_lru_cache 1.6.4              pyhd8ed1ab_0    conda-forge\n",
      "biopython                 1.79             py39h3811e60_0    conda-forge\n",
      "bleach                    4.1.0              pyhd8ed1ab_0    conda-forge\n",
      "bokeh                     2.4.1            py39hf3d152e_1    conda-forge\n",
      "boto                      2.49.0                     py_0    conda-forge\n",
      "boto3                     1.19.2                   pypi_0    pypi\n",
      "botocore                  1.22.2                   pypi_0    pypi\n",
      "brotlipy                  0.7.0           py39h3811e60_1001    conda-forge\n",
      "bzip2                     1.0.8                h7f98852_4    conda-forge\n",
      "c-ares                    1.17.2               h7f98852_0    conda-forge\n",
      "ca-certificates           2022.6.15            ha878542_0    conda-forge\n",
      "cachecontrol              0.12.6                     py_0    conda-forge\n",
      "cairo                     1.16.0            h6cf1ce9_1008    conda-forge\n",
      "certifi                   2022.6.15        py39hf3d152e_0    conda-forge\n",
      "cffi                      1.14.6           py39h4bc2ebd_1    conda-forge\n",
      "chardet                   4.0.0            py39hf3d152e_1    conda-forge\n",
      "charset-normalizer        2.0.0              pyhd8ed1ab_0    conda-forge\n",
      "click                     8.0.3            py39hf3d152e_0    conda-forge\n",
      "colorama                  0.4.4              pyh9f0ad1d_0    conda-forge\n",
      "cryptography              35.0.0           py39h95dcef6_1    conda-forge\n",
      "cycler                    0.10.0                     py_2    conda-forge\n",
      "cython                    0.29.24          py39he80948d_0    conda-forge\n",
      "debugpy                   1.4.1            py39he80948d_0    conda-forge\n",
      "decorator                 4.4.2                      py_0    conda-forge\n",
      "defusedxml                0.7.1              pyhd8ed1ab_0    conda-forge\n",
      "english                   2020.7.0           pyhd8ed1ab_0    conda-forge\n",
      "entrypoints               0.3             py39hde42818_1002    conda-forge\n",
      "et_xmlfile                1.0.1                   py_1001    conda-forge\n",
      "expat                     2.4.1                h9c3ff4c_0    conda-forge\n",
      "font-ttf-dejavu-sans-mono 2.37                 hab24e00_0    conda-forge\n",
      "font-ttf-inconsolata      3.000                h77eed37_0    conda-forge\n",
      "font-ttf-source-code-pro  2.038                h77eed37_0    conda-forge\n",
      "font-ttf-ubuntu           0.83                 hab24e00_0    conda-forge\n",
      "fontconfig                2.13.1            hba837de_1005    conda-forge\n",
      "fonts-conda-ecosystem     1                             0    conda-forge\n",
      "fonts-conda-forge         1                             0    conda-forge\n",
      "freetype                  2.10.4               h0708190_1    conda-forge\n",
      "fribidi                   1.0.10               h36c2ea0_0    conda-forge\n",
      "future                    0.18.2           py39hf3d152e_3    conda-forge\n",
      "gdk-pixbuf                2.42.6               h04a7f16_0    conda-forge\n",
      "gettext                   0.19.8.1          h73d1719_1008    conda-forge\n",
      "giflib                    5.2.1                h36c2ea0_2    conda-forge\n",
      "graphite2                 1.3.13            h58526e2_1001    conda-forge\n",
      "graphviz                  2.49.1               h85b4f2f_0    conda-forge\n",
      "gtk2                      2.24.33              h539f30e_1    conda-forge\n",
      "gts                       0.7.6                h64030ff_2    conda-forge\n",
      "harfbuzz                  3.0.0                h83ec7ef_1    conda-forge\n",
      "hdmedians                 0.14.2           py39hce5d2b2_0    conda-forge\n",
      "icu                       68.1                 h58526e2_0    conda-forge\n",
      "idna                      3.1                pyhd3deb0d_0    conda-forge\n",
      "importlib-metadata        4.8.1            py39hf3d152e_0    conda-forge\n",
      "iniconfig                 1.1.1              pyh9f0ad1d_0    conda-forge\n",
      "interchange               2021.0.4           pyhd8ed1ab_0    conda-forge\n",
      "ipykernel                 6.4.2            py39hef51801_0    conda-forge\n",
      "ipython                   7.28.0           py39hef51801_0    conda-forge\n",
      "ipython_genutils          0.2.0                      py_1    conda-forge\n",
      "ipywidgets                7.6.5              pyhd8ed1ab_0    conda-forge\n",
      "jbig                      2.1               h7f98852_2003    conda-forge\n",
      "jedi                      0.18.0           py39hf3d152e_2    conda-forge\n",
      "jinja2                    3.0.2              pyhd8ed1ab_0    conda-forge\n",
      "jmespath                  0.10.0                   pypi_0    pypi\n",
      "joblib                    1.1.0              pyhd8ed1ab_0    conda-forge\n",
      "jpeg                      9d                   h36c2ea0_0    conda-forge\n",
      "jsonschema                4.1.2              pyhd8ed1ab_0    conda-forge\n",
      "jupyter_client            7.0.6              pyhd8ed1ab_0    conda-forge\n",
      "jupyter_core              4.8.1            py39hf3d152e_0    conda-forge\n",
      "jupyterlab_pygments       0.1.2              pyh9f0ad1d_0    conda-forge\n",
      "jupyterlab_widgets        1.0.2              pyhd8ed1ab_0    conda-forge\n",
      "kiwisolver                1.3.2            py39h1a9c180_0    conda-forge\n",
      "krb5                      1.19.2               hcc1bbae_2    conda-forge\n",
      "lcms2                     2.12                 hddcbb42_0    conda-forge\n",
      "ld_impl_linux-64          2.36.1               hea4e1c9_2    conda-forge\n",
      "lerc                      2.2.1                h9c3ff4c_0    conda-forge\n",
      "libblas                   3.9.0           12_linux64_openblas    conda-forge\n",
      "libcblas                  3.9.0           12_linux64_openblas    conda-forge\n",
      "libcurl                   7.79.1               h2574ce0_1    conda-forge\n",
      "libdeflate                1.7                  h7f98852_5    conda-forge\n",
      "libedit                   3.1.20191231         he28a2e2_2    conda-forge\n",
      "libev                     4.33                 h516909a_1    conda-forge\n",
      "libffi                    3.4.2                h9c3ff4c_4    conda-forge\n",
      "libgcc-ng                 11.2.0              h1d223b6_11    conda-forge\n",
      "libgd                     2.3.3                h6ad9fb6_0    conda-forge\n",
      "libgfortran-ng            11.2.0              h69a702a_11    conda-forge\n",
      "libgfortran5              11.2.0              h5c6108e_11    conda-forge\n",
      "libglib                   2.70.0               h174f98d_1    conda-forge\n",
      "libgomp                   11.2.0              h1d223b6_11    conda-forge\n",
      "libiconv                  1.16                 h516909a_0    conda-forge\n",
      "liblapack                 3.9.0           12_linux64_openblas    conda-forge\n",
      "libnghttp2                1.43.0               h812cca2_1    conda-forge\n",
      "libopenblas               0.3.18          pthreads_h8fe5266_0    conda-forge\n",
      "libpng                    1.6.37               h21135ba_2    conda-forge\n",
      "libpq                     13.3                 hd57d9b9_1    conda-forge\n",
      "librsvg                   2.52.2               hc3c00ef_0    conda-forge\n",
      "libsodium                 1.0.18               h36c2ea0_1    conda-forge\n",
      "libssh2                   1.10.0               ha56f1ee_2    conda-forge\n",
      "libstdcxx-ng              11.2.0              he4da1e4_11    conda-forge\n",
      "libtiff                   4.3.0                hf544144_1    conda-forge\n",
      "libtool                   2.4.6             h9c3ff4c_1008    conda-forge\n",
      "libuuid                   2.32.1            h7f98852_1000    conda-forge\n",
      "libwebp                   1.2.1                h3452ae3_0    conda-forge\n",
      "libwebp-base              1.2.1                h7f98852_0    conda-forge\n",
      "libxcb                    1.13              h7f98852_1003    conda-forge\n",
      "libxml2                   2.9.12               h72842e0_0    conda-forge\n",
      "libxslt                   1.1.33               h15afd5d_2    conda-forge\n",
      "libzlib                   1.2.11            h36c2ea0_1013    conda-forge\n",
      "lockfile                  0.12.2                     py_1    conda-forge\n",
      "lxml                      4.6.3            py39h107f48f_0    conda-forge\n",
      "lz4-c                     1.9.3                h9c3ff4c_1    conda-forge\n",
      "markupsafe                2.0.1            py39h3811e60_0    conda-forge\n",
      "matplotlib-base           3.4.3            py39h2fa2bec_1    conda-forge\n",
      "matplotlib-inline         0.1.3              pyhd8ed1ab_0    conda-forge\n",
      "matplotlib-venn           0.11.7             pyhd8ed1ab_0    conda-forge\n",
      "missingpy                 0.2.0                    pypi_0    pypi\n",
      "mistune                   0.8.4           py39h3811e60_1004    conda-forge\n",
      "more-itertools            8.10.0             pyhd8ed1ab_0    conda-forge\n",
      "mscorefonts               0.0.1                         3    conda-forge\n",
      "msgpack-python            1.0.2            py39h1a9c180_1    conda-forge\n",
      "natsort                   7.1.1              pyhd8ed1ab_0    conda-forge\n",
      "nbclient                  0.5.4              pyhd8ed1ab_0    conda-forge\n",
      "nbconvert                 6.2.0            py39hf3d152e_0    conda-forge\n",
      "nbformat                  5.1.3              pyhd8ed1ab_0    conda-forge\n",
      "ncurses                   6.2                  h58526e2_4    conda-forge\n",
      "neo4j-python-driver       4.3.7              pyhd8ed1ab_0    conda-forge\n",
      "neobolt                   1.7.17           py39h3811e60_2    conda-forge\n",
      "neotime                   1.7.4                      py_0    conda-forge\n",
      "nest-asyncio              1.5.1              pyhd8ed1ab_0    conda-forge\n",
      "networkx                  2.6.3              pyhd8ed1ab_0    conda-forge\n",
      "notebook                  6.4.5              pyha770c72_0    conda-forge\n",
      "numpy                     1.21.3           py39hdbf815f_0    conda-forge\n",
      "olefile                   0.46               pyh9f0ad1d_1    conda-forge\n",
      "openjpeg                  2.4.0                hb52868f_1    conda-forge\n",
      "openpyxl                  3.0.9              pyhd8ed1ab_0    conda-forge\n",
      "openssl                   1.1.1o               h166bdaf_0    conda-forge\n",
      "packaging                 21.0               pyhd8ed1ab_0    conda-forge\n",
      "pandas                    1.3.4            py39hde0f152_0    conda-forge\n",
      "pandoc                    2.15                 h7f98852_0    conda-forge\n",
      "pandocfilters             1.5.0              pyhd8ed1ab_0    conda-forge\n",
      "pango                     1.48.10              h54213e6_2    conda-forge\n",
      "parso                     0.8.2              pyhd8ed1ab_0    conda-forge\n",
      "patsy                     0.5.2              pyhd8ed1ab_0    conda-forge\n",
      "pcre                      8.45                 h9c3ff4c_0    conda-forge\n",
      "pexpect                   4.8.0              pyh9f0ad1d_2    conda-forge\n",
      "pickleshare               0.7.5           py39hde42818_1002    conda-forge\n",
      "pillow                    8.3.2            py39ha612740_0    conda-forge\n",
      "pip                       21.3.1             pyhd8ed1ab_0    conda-forge\n",
      "pixman                    0.40.0               h36c2ea0_0    conda-forge\n",
      "plotly                    5.3.1              pyhd8ed1ab_0    conda-forge\n",
      "pluggy                    1.0.0            py39hf3d152e_1    conda-forge\n",
      "prometheus_client         0.11.0             pyhd8ed1ab_0    conda-forge\n",
      "prompt-toolkit            3.0.21             pyha770c72_0    conda-forge\n",
      "prompt_toolkit            3.0.21               hd8ed1ab_0    conda-forge\n",
      "psycopg2                  2.9.1            py39h3811e60_0    conda-forge\n",
      "pthread-stubs             0.4               h36c2ea0_1001    conda-forge\n",
      "ptyprocess                0.7.0              pyhd3deb0d_0    conda-forge\n",
      "py                        1.10.0             pyhd3deb0d_0    conda-forge\n",
      "py2neo                    2021.2.3           pyhd8ed1ab_0    conda-forge\n",
      "pycparser                 2.20               pyh9f0ad1d_2    conda-forge\n",
      "pygments                  2.10.0             pyhd8ed1ab_0    conda-forge\n",
      "pygraphviz                1.7              py39h78163bd_0    conda-forge\n",
      "pyopenssl                 21.0.0             pyhd8ed1ab_0    conda-forge\n",
      "pyparsing                 2.4.7              pyh9f0ad1d_0    conda-forge\n",
      "pyrsistent                0.17.3           py39h3811e60_2    conda-forge\n",
      "pysam                     0.17.0           py39h051187c_0    bioconda\n",
      "pysocks                   1.7.1            py39hf3d152e_3    conda-forge\n",
      "pytabix                   0.1              py39h98c8e45_1    bioconda\n",
      "pytest                    6.2.5            py39hf3d152e_0    conda-forge\n",
      "python                    3.9.7           hb7a2778_3_cpython    conda-forge\n",
      "python-dateutil           2.8.2              pyhd8ed1ab_0    conda-forge\n",
      "python-graphviz           0.17               pyhaef67bd_0    conda-forge\n",
      "python_abi                3.9                      2_cp39    conda-forge\n",
      "pytz                      2021.3             pyhd8ed1ab_0    conda-forge\n",
      "pyvcf                     0.6.8           py39hde42818_1002    conda-forge\n",
      "pyyaml                    6.0              py39h3811e60_0    conda-forge\n",
      "pyzmq                     22.3.0           py39h37b5a0c_0    conda-forge\n",
      "readline                  8.1                  h46c0cb4_0    conda-forge\n",
      "requests                  2.26.0             pyhd8ed1ab_0    conda-forge\n",
      "s3transfer                0.5.0                    pypi_0    pypi\n",
      "scikit-bio                0.5.6            py39h16ac069_4    conda-forge\n",
      "scikit-learn              1.0.1            py39h7c5d8c9_0    conda-forge\n",
      "scipy                     1.7.1            py39hee8e79c_0    conda-forge\n",
      "seaborn                   0.11.2               hd8ed1ab_0    conda-forge\n",
      "seaborn-base              0.11.2             pyhd8ed1ab_0    conda-forge\n",
      "send2trash                1.8.0              pyhd8ed1ab_0    conda-forge\n",
      "setuptools                58.2.0           py39hf3d152e_0    conda-forge\n",
      "six                       1.16.0             pyh6c4a22f_0    conda-forge\n",
      "sqlite                    3.36.0               h9cd32fc_2    conda-forge\n",
      "statsmodels               0.13.0           py39hce5d2b2_0    conda-forge\n",
      "tenacity                  8.0.1              pyhd8ed1ab_0    conda-forge\n",
      "terminado                 0.12.1           py39hf3d152e_0    conda-forge\n",
      "testpath                  0.5.0              pyhd8ed1ab_0    conda-forge\n",
      "threadpoolctl             3.0.0              pyh8a188c0_0    conda-forge\n",
      "tk                        8.6.11               h27826a3_1    conda-forge\n",
      "toml                      0.10.2             pyhd8ed1ab_0    conda-forge\n",
      "tornado                   6.1              py39h3811e60_1    conda-forge\n",
      "traitlets                 5.1.0              pyhd8ed1ab_0    conda-forge\n",
      "typing_extensions         3.10.0.2           pyha770c72_0    conda-forge\n",
      "tzdata                    2021e                he74cb21_0    conda-forge\n",
      "urllib3                   1.26.7             pyhd8ed1ab_0    conda-forge\n",
      "venn                      0.1.3                    pypi_0    pypi\n",
      "wcwidth                   0.2.5              pyh9f0ad1d_2    conda-forge\n",
      "webencodings              0.5.1                      py_1    conda-forge\n",
      "wheel                     0.37.0             pyhd8ed1ab_1    conda-forge\n",
      "widgetsnbextension        3.5.1            py39hf3d152e_4    conda-forge\n",
      "wordcloud                 1.8.1                    pypi_0    pypi\n",
      "xlrd                      2.0.1              pyhd8ed1ab_3    conda-forge\n",
      "xorg-kbproto              1.0.7             h7f98852_1002    conda-forge\n",
      "xorg-libice               1.0.10               h7f98852_0    conda-forge\n",
      "xorg-libsm                1.2.3             hd9c2040_1000    conda-forge\n",
      "xorg-libx11               1.7.2                h7f98852_0    conda-forge\n",
      "xorg-libxau               1.0.9                h7f98852_0    conda-forge\n",
      "xorg-libxdmcp             1.1.3                h7f98852_0    conda-forge\n",
      "xorg-libxext              1.3.4                h7f98852_1    conda-forge\n",
      "xorg-libxrender           0.9.10            h7f98852_1003    conda-forge\n",
      "xorg-renderproto          0.11.1            h7f98852_1002    conda-forge\n",
      "xorg-xextproto            7.3.0             h7f98852_1002    conda-forge\n",
      "xorg-xproto               7.0.31            h7f98852_1007    conda-forge\n",
      "xz                        5.2.5                h516909a_1    conda-forge\n",
      "yaml                      0.2.5                h516909a_0    conda-forge\n",
      "zeromq                    4.3.4                h9c3ff4c_1    conda-forge\n",
      "zipp                      3.6.0              pyhd8ed1ab_0    conda-forge\n",
      "zlib                      1.2.11            h36c2ea0_1013    conda-forge\n",
      "zstd                      1.5.0                ha95c52a_0    conda-forge\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scipy.stats as stats\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import seaborn as sns\n",
    "#For Arial font\n",
    "#!conda install -c conda-forge -y mscorefonts\n",
    "##-> The below was also needed in matplotlib 3.4.2\n",
    "#import shutil\n",
    "#import matplotlib\n",
    "#shutil.rmtree(matplotlib.get_cachedir())\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "from IPython.display import display\n",
    "import time\n",
    "\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.linear_model import RidgeCV\n",
    "from sklearn.linear_model import LinearRegression\n",
    "from statsmodels.stats import weightstats\n",
    "from statsmodels.stats import multitest as multi\n",
    "from decimal import Decimal, ROUND_HALF_UP\n",
    "\n",
    "!conda list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Data preparation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> The following code is completely same with the one used in the baseline LASSO modeling. Hence, the correspondence between participant and testing (hold-out) set for each model is maintained.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1-1. Import the cleaned dataframes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Import the baseline BMI dataframe\n",
    "fileDir = '../210104_Biological-BMI-paper/ExportData/'\n",
    "ipynbName = '210104_Biological-BMI-paper_RF-imputation_'\n",
    "fileName = 'baseline-combiDF-with-RF-imputation.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",
    "##Take BMI and general covariates (without Race in this study)\n",
    "tempL = ['log_BaseBMI', 'Sex', 'BaseAge', 'PC1', 'PC2', 'PC3', 'PC4', 'PC5']\n",
    "tempDF = tempDF[tempL]\n",
    "\n",
    "display(tempDF)\n",
    "\n",
    "bmiDF = tempDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Import the cleaned baseline omics dataframes\n",
    "fileDir = '../210104_Biological-BMI-paper/ExportData/'\n",
    "ipynbName = '210104_Biological-BMI-paper_RF-imputation_'\n",
    "tempL = ['log_BaseBMI', 'Sex', 'BaseAge', 'PC1', 'PC2', 'PC3', 'PC4', 'PC5', 'Race']\n",
    "tempD = {}\n",
    "for df_n in ['metDF', 'chemDF', 'protDF', 'combiDF']:\n",
    "    fileName = 'baseline-'+df_n+'-with-RF-imputation.tsv'\n",
    "    tempDF = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t', dtype={'public_client_id': str})\n",
    "    tempDF = tempDF.set_index('public_client_id')\n",
    "    print(df_n+' original shape:', tempDF.shape)\n",
    "    #Drop BMI and covariates\n",
    "    tempDF = tempDF.drop(columns=tempL)\n",
    "    display(tempDF)\n",
    "    tempD[df_n] = tempDF\n",
    "\n",
    "metDF = tempD['metDF']\n",
    "chemDF = tempD['chemDF']\n",
    "protDF = tempD['protDF']\n",
    "combiDF = tempD['combiDF']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1-2. Stratification with sex"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Stratify the cohort with sex\n",
    "bmiDF_F = bmiDF.loc[bmiDF['Sex']=='F']\n",
    "bmiDF_M = bmiDF.loc[bmiDF['Sex']=='M']\n",
    "bmiDF_B = bmiDF#Not copy just rename\n",
    "print('Female, Male, Both sex = ', len(bmiDF_F), ', ', len(bmiDF_M), ', ', len(bmiDF_B))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1-3. Split the cohort into 10 sets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempD1 = {'Female':bmiDF_F, 'Male':bmiDF_M, 'BothSex':bmiDF_B}\n",
    "nmodels = 10\n",
    "tempD2 = {}\n",
    "for sex in tempD1.keys():\n",
    "    tempDF = tempD1[sex]\n",
    "    #Split cohort to define the training and testing (hold-out) sets\n",
    "    tempL = np.array_split(tempDF, nmodels)#List of DFs\n",
    "    tempD = {}\n",
    "    for model_k in range(nmodels):\n",
    "        tempDF1 = tempL[model_k]\n",
    "        model_n = 'Model_'+str(model_k+1).zfill(2)\n",
    "        tempS = pd.Series(np.repeat(model_n, len(tempDF1)),\n",
    "                          index=tempDF1.index, name='Testing')\n",
    "        tempD[model_k] = tempS\n",
    "    tempS = pd.concat(list(tempD.values()), axis=0)\n",
    "    #Add the info to bmiDF\n",
    "    tempDF = pd.merge(tempDF, tempS, left_index=True, right_index=True, how='left')\n",
    "    tempD2[sex] = tempDF\n",
    "    print(sex)\n",
    "    display(tempDF)\n",
    "    display(tempDF['Testing'].value_counts())\n",
    "    print('')\n",
    "#Update\n",
    "bmiDF_F = tempD2['Female']\n",
    "bmiDF_M = tempD2['Male']\n",
    "bmiDF_B = tempD2['BothSex']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1-4. Check independency among the sets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> See the notebook for LASSO.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Metabolomics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2-1. Standardization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Note that all preprocessing, including outlier elimination, missingness filtering, imputation, and standardization, should be performed within the cross-validation fold, not across whole dataset, to minimize the potential data leakage. In this study, however, an external dataset is used for “validation\" of the findings (i.e., a testing set for the fitted models). Hence, the robustness of preprocessing is prioritized.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF = metDF\n",
    "tempD1 = {'Female':bmiDF_F, 'Male':bmiDF_M, 'BothSex':bmiDF_B}\n",
    "tempD2 = {}\n",
    "for sex in tempD1.keys():\n",
    "    tempDF1 = tempD1[sex]\n",
    "    tempDF1 = tempDF.loc[tempDF1.index.tolist()]\n",
    "    #Z-score transformation\n",
    "    scaler = StandardScaler(copy=True, with_mean=True, with_std=True)\n",
    "    tempA = scaler.fit_transform(tempDF1)#Column direction\n",
    "    tempDF1 = pd.DataFrame(data=tempA, index=tempDF1.index, columns=tempDF1.columns)\n",
    "    tempD2[sex] = tempDF1\n",
    "    \n",
    "    #Confirmation\n",
    "    print(sex+':', tempDF1.shape)\n",
    "    display(tempDF1.describe())\n",
    "    sns.set(style='ticks', font='Arial', context='notebook')\n",
    "    plt.figure(figsize=(4, 3))\n",
    "    for col_i in range(3):\n",
    "        sns.distplot(tempDF1.iloc[:, col_i], label=tempDF1.columns[col_i])\n",
    "    sns.despine()\n",
    "    plt.xlabel(r'$Z$'+'-score')\n",
    "    plt.ylabel('Density')\n",
    "    plt.legend(bbox_to_anchor=(1, 0.5), loc='center left', borderaxespad=1)\n",
    "    plt.show()\n",
    "    print('')\n",
    "\n",
    "metDF_F = tempD2['Female']\n",
    "metDF_M = tempD2['Male']\n",
    "metDF_B = tempD2['BothSex']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2-2. Ridge with cross-validation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Note that, in contrast to LassoCV and similarly to LinearRegression, the current RidgeCV seems to produce 2D ndarray although target y is 1D.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Female model\n",
    "tempDF1 = metDF_F#Standardized independent variables\n",
    "tempDF2 = bmiDF_F#Not-standardized dependent variable and info about testing set\n",
    "yvar = 'BaseBMI'\n",
    "nmodels = len(tempDF2['Testing'].unique())#Reset just in case\n",
    "ncvs = 10\n",
    "model = RidgeCV(alphas=np.geomspace(1, 1000, num=200), fit_intercept=True,\n",
    "                normalize=False, gcv_mode='auto', cv=ncvs)\n",
    "yvar_model = 'BaseMetBMI'\n",
    "\n",
    "#Perform Ridge\n",
    "bcoefDF = pd.DataFrame(index=tempDF1.columns).astype('float64')#For beta-coefficients\n",
    "bcoefDF.index.rename('Variable', inplace=True)\n",
    "interceptDF = pd.DataFrame(index=['Intercept']).astype('float64')#For intercept\n",
    "interceptDF.index.rename('Variable', inplace=True)\n",
    "scoreL = []#For the coefficient of determination R2\n",
    "predictS = pd.Series(name='log_'+yvar_model).astype('float64')#For predictions\n",
    "t_start = time.time()\n",
    "for model_k in range(nmodels):\n",
    "    #Prepare training and testing (hold-out) datasets in model k\n",
    "    model_n = 'Model_'+str(model_k+1).zfill(2)\n",
    "    yDF_train = tempDF2.loc[tempDF2['Testing']!=model_n]\n",
    "    yDF_test = tempDF2.loc[tempDF2['Testing']==model_n]\n",
    "    xDF_train = tempDF1.loc[yDF_train.index.tolist()]\n",
    "    xDF_test = tempDF1.loc[yDF_test.index.tolist()]\n",
    "    #Retrieve the dependent variable\n",
    "    yDF_train = pd.DataFrame(yDF_train['log_'+yvar])#Not Series but DF\n",
    "    yDF_test = pd.DataFrame(yDF_test['log_'+yvar])#Not Series but DF\n",
    "    #Fitting model with cross-validation using training dataset (i.e., internal training/validation datasets)\n",
    "    model.fit(xDF_train, yDF_train)\n",
    "    #Check the penalization amount decided by cross validation\n",
    "    print(model_n+': alpha =', model.alpha_)\n",
    "    #Save parameters\n",
    "    bcoefDF[model_n] = model.coef_.flatten()#weighted vector\n",
    "    interceptDF[model_n] = model.intercept_#Independent term in decision function\n",
    "    #Evaluation with testing (hold-out) dataset\n",
    "    scoreL.append(model.score(xDF_test, yDF_test))\n",
    "    #Prediction for testing dataset using the fitted model k\n",
    "    tempS = pd.Series(model.predict(xDF_test).flatten(),\n",
    "                      index=xDF_test.index, name='log_'+yvar_model)\n",
    "    predictS = pd.concat([predictS, tempS], axis=0)\n",
    "t_elapsed = time.time() - t_start\n",
    "print('Elapsed time for '+str(nmodels)+' models of '+str(ncvs)+'-fold CV Ridge:',\n",
    "      round(t_elapsed//60), 'min', round(t_elapsed%60, 1), 'sec')\n",
    "\n",
    "#Clean prediction DF\n",
    "tempDF = pd.merge(tempDF2[['Testing', 'log_'+yvar]], predictS,\n",
    "                  left_index=True, right_index=True, how='left')\n",
    "##Convert to original scale\n",
    "tempDF[yvar] = np.e**tempDF['log_'+yvar]\n",
    "tempDF[yvar_model] = np.e**tempDF['log_'+yvar_model]\n",
    "display(tempDF)\n",
    "\n",
    "#Save the cleaned prediction DF\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-Ridge_'\n",
    "fileName = yvar_model.replace('Base', '')+'-Female.tsv'\n",
    "tempDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "metBMI_F_bcoefs = bcoefDF\n",
    "metBMI_F_intercept = interceptDF\n",
    "metBMI_F_R2 = scoreL\n",
    "metBMI_F = tempDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Male model\n",
    "tempDF1 = metDF_M#Standardized independent variables\n",
    "tempDF2 = bmiDF_M#Not-standardized dependent variable and info about testing set\n",
    "yvar = 'BaseBMI'\n",
    "nmodels = len(tempDF2['Testing'].unique())#Reset just in case\n",
    "ncvs = 10\n",
    "model = RidgeCV(alphas=np.geomspace(1, 1000, num=200), fit_intercept=True,\n",
    "                normalize=False, gcv_mode='auto', cv=ncvs)\n",
    "yvar_model = 'BaseMetBMI'\n",
    "\n",
    "#Perform Ridge\n",
    "bcoefDF = pd.DataFrame(index=tempDF1.columns).astype('float64')#For beta-coefficients\n",
    "bcoefDF.index.rename('Variable', inplace=True)\n",
    "interceptDF = pd.DataFrame(index=['Intercept']).astype('float64')#For intercept\n",
    "interceptDF.index.rename('Variable', inplace=True)\n",
    "scoreL = []#For the coefficient of determination R2\n",
    "predictS = pd.Series(name='log_'+yvar_model).astype('float64')#For predictions\n",
    "t_start = time.time()\n",
    "for model_k in range(nmodels):\n",
    "    #Prepare training and testing (hold-out) datasets in model k\n",
    "    model_n = 'Model_'+str(model_k+1).zfill(2)\n",
    "    yDF_train = tempDF2.loc[tempDF2['Testing']!=model_n]\n",
    "    yDF_test = tempDF2.loc[tempDF2['Testing']==model_n]\n",
    "    xDF_train = tempDF1.loc[yDF_train.index.tolist()]\n",
    "    xDF_test = tempDF1.loc[yDF_test.index.tolist()]\n",
    "    #Retrieve the dependent variable\n",
    "    yDF_train = pd.DataFrame(yDF_train['log_'+yvar])#Not Series but DF\n",
    "    yDF_test = pd.DataFrame(yDF_test['log_'+yvar])#Not Series but DF\n",
    "    #Fitting model with cross-validation using training dataset (i.e., internal training/validation datasets)\n",
    "    model.fit(xDF_train, yDF_train)\n",
    "    #Check the penalization amount decided by cross validation\n",
    "    print(model_n+': alpha =', model.alpha_)\n",
    "    #Save parameters\n",
    "    bcoefDF[model_n] = model.coef_.flatten()#weighted vector\n",
    "    interceptDF[model_n] = model.intercept_#Independent term in decision function\n",
    "    #Evaluation with testing (hold-out) dataset\n",
    "    scoreL.append(model.score(xDF_test, yDF_test))\n",
    "    #Prediction for testing dataset using the fitted model k\n",
    "    tempS = pd.Series(model.predict(xDF_test).flatten(),\n",
    "                      index=xDF_test.index, name='log_'+yvar_model)\n",
    "    predictS = pd.concat([predictS, tempS], axis=0)\n",
    "t_elapsed = time.time() - t_start\n",
    "print('Elapsed time for '+str(nmodels)+' models of '+str(ncvs)+'-fold CV Ridge:',\n",
    "      round(t_elapsed//60), 'min', round(t_elapsed%60, 1), 'sec')\n",
    "\n",
    "#Clean prediction DF\n",
    "tempDF = pd.merge(tempDF2[['Testing', 'log_'+yvar]], predictS,\n",
    "                  left_index=True, right_index=True, how='left')\n",
    "##Convert to original scale\n",
    "tempDF[yvar] = np.e**tempDF['log_'+yvar]\n",
    "tempDF[yvar_model] = np.e**tempDF['log_'+yvar_model]\n",
    "display(tempDF)\n",
    "\n",
    "#Save the cleaned prediction DF\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-Ridge_'\n",
    "fileName = yvar_model.replace('Base', '')+'-Male.tsv'\n",
    "tempDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "metBMI_M_bcoefs = bcoefDF\n",
    "metBMI_M_intercept = interceptDF\n",
    "metBMI_M_R2 = scoreL\n",
    "metBMI_M = tempDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Both sex model\n",
    "tempDF1 = metDF_B#Standardized independent variables\n",
    "tempDF2 = bmiDF_B#Not-standardized dependent variable and info about testing set\n",
    "yvar = 'BaseBMI'\n",
    "nmodels = len(tempDF2['Testing'].unique())#Reset just in case\n",
    "ncvs = 10\n",
    "model = RidgeCV(alphas=np.geomspace(1, 1000, num=200), fit_intercept=True,\n",
    "                normalize=False, gcv_mode='auto', cv=ncvs)\n",
    "yvar_model = 'BaseMetBMI'\n",
    "\n",
    "#Perform Ridge\n",
    "bcoefDF = pd.DataFrame(index=tempDF1.columns).astype('float64')#For beta-coefficients\n",
    "bcoefDF.index.rename('Variable', inplace=True)\n",
    "interceptDF = pd.DataFrame(index=['Intercept']).astype('float64')#For intercept\n",
    "interceptDF.index.rename('Variable', inplace=True)\n",
    "scoreL = []#For the coefficient of determination R2\n",
    "predictS = pd.Series(name='log_'+yvar_model).astype('float64')#For predictions\n",
    "t_start = time.time()\n",
    "for model_k in range(nmodels):\n",
    "    #Prepare training and testing (hold-out) datasets in model k\n",
    "    model_n = 'Model_'+str(model_k+1).zfill(2)\n",
    "    yDF_train = tempDF2.loc[tempDF2['Testing']!=model_n]\n",
    "    yDF_test = tempDF2.loc[tempDF2['Testing']==model_n]\n",
    "    xDF_train = tempDF1.loc[yDF_train.index.tolist()]\n",
    "    xDF_test = tempDF1.loc[yDF_test.index.tolist()]\n",
    "    #Retrieve the dependent variable\n",
    "    yDF_train = pd.DataFrame(yDF_train['log_'+yvar])#Not Series but DF\n",
    "    yDF_test = pd.DataFrame(yDF_test['log_'+yvar])#Not Series but DF\n",
    "    #Fitting model with cross-validation using training dataset (i.e., internal training/validation datasets)\n",
    "    model.fit(xDF_train, yDF_train)\n",
    "    #Check the penalization amount decided by cross validation\n",
    "    print(model_n+': alpha =', model.alpha_)\n",
    "    #Save parameters\n",
    "    bcoefDF[model_n] = model.coef_.flatten()#weighted vector\n",
    "    interceptDF[model_n] = model.intercept_#Independent term in decision function\n",
    "    #Evaluation with testing (hold-out) dataset\n",
    "    scoreL.append(model.score(xDF_test, yDF_test))\n",
    "    #Prediction for testing dataset using the fitted model k\n",
    "    tempS = pd.Series(model.predict(xDF_test).flatten(),\n",
    "                      index=xDF_test.index, name='log_'+yvar_model)\n",
    "    predictS = pd.concat([predictS, tempS], axis=0)\n",
    "t_elapsed = time.time() - t_start\n",
    "print('Elapsed time for '+str(nmodels)+' models of '+str(ncvs)+'-fold CV Ridge:',\n",
    "      round(t_elapsed//60), 'min', round(t_elapsed%60, 1), 'sec')\n",
    "\n",
    "#Clean prediction DF\n",
    "tempDF = pd.merge(tempDF2[['Testing', 'log_'+yvar]], predictS,\n",
    "                  left_index=True, right_index=True, how='left')\n",
    "##Convert to original scale\n",
    "tempDF[yvar] = np.e**tempDF['log_'+yvar]\n",
    "tempDF[yvar_model] = np.e**tempDF['log_'+yvar_model]\n",
    "display(tempDF)\n",
    "\n",
    "#Save the cleaned prediction DF\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-Ridge_'\n",
    "fileName = yvar_model.replace('Base', '')+'-BothSex.tsv'\n",
    "tempDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "metBMI_B_bcoefs = bcoefDF\n",
    "metBMI_B_intercept = interceptDF\n",
    "metBMI_B_R2 = scoreL\n",
    "metBMI_B = tempDF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2-3. Prediction accuracy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Summary\n",
    "tempD1 = {'Female':metBMI_F_R2, 'Male':metBMI_M_R2, 'Both sex':metBMI_B_R2}\n",
    "tempD2 = {'Female':metBMI_F, 'Male':metBMI_M, 'Both sex':metBMI_B}\n",
    "yvar = 'BaseBMI'\n",
    "yvar_model = 'BaseMetBMI'\n",
    "\n",
    "for sex in tempD1.keys():\n",
    "    tempL = tempD1[sex]\n",
    "    print(sex+' 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",
    "    tempDF = tempD2[sex]\n",
    "    print(' - Observed vs. predicted log_'+yvar+': (Pearson\\'s r, P) =',\n",
    "          stats.pearsonr(tempDF['log_'+yvar], tempDF['log_'+yvar_model]))\n",
    "    print(' - Observed vs. predicted '+yvar+': (Pearson\\'s r, P) =',\n",
    "          stats.pearsonr(tempDF[yvar], tempDF[yvar_model]))\n",
    "    display(tempDF.describe())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Check difference between sex-specific and sex-mixed models for now. Note that this is a rough comparison because the sample size is different.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Prepare DF\n",
    "tempDF = pd.DataFrame({'Female':metBMI_F_R2, 'Male':metBMI_M_R2, 'Both sex':metBMI_B_R2})\n",
    "tempDF = tempDF.melt(var_name='Cohort', value_name='R2', value_vars=tempDF.columns.tolist())\n",
    "\n",
    "#Plot\n",
    "sns.set(style='ticks', font='Arial', context='notebook')\n",
    "plt.figure(figsize=(3, 1))\n",
    "sns.barplot(data=tempDF, y='Cohort', x='R2', hue='Cohort', palette='Set1', dodge=False, edgecolor='black',\n",
    "            ci=95, capsize=0.4, errwidth=1.5, errcolor='black')\n",
    "sns.stripplot(data=tempDF, y='Cohort', x='R2', hue='Cohort', dodge=False, size=8, edgecolor='black',\n",
    "              linewidth=1, alpha=0.4, palette={'Female':'gray', 'Male':'gray', 'Both sex':'gray'})\n",
    "sns.despine()\n",
    "plt.xlabel('Out-of-sample '+r'$R^2$'+' in 10 Ridge models\\n[mean with 95% CI]')\n",
    "plt.ylabel('')\n",
    "plt.legend('', frameon=False)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempD1 = {'Female':metBMI_F, 'Male':metBMI_M, 'Both sex':metBMI_B}\n",
    "yvar = 'BaseBMI'\n",
    "yvar_model = 'BaseMetBMI'\n",
    "range_min = np.min([df[var].min() for df in tempD1.values() for var in [yvar, yvar_model]])\n",
    "range_max = np.max([df[var].max() for df in tempD1.values() for var in [yvar, yvar_model]])\n",
    "axis_label = 'BMI [kg m'+r'$^{-2}$'+']'\n",
    "\n",
    "#Plot measured vs. predicted per model\n",
    "tempD2 = {'Female':'tab:red', 'Male':'tab:blue', 'Both sex':'tab:green'}\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "fig, axes = plt.subplots(nrows=1, ncols=len(tempD2), figsize=(10, 3), sharex=True, sharey=True,\n",
    "                         gridspec_kw={'hspace':0.1, 'wspace':0.1})\n",
    "for ax_i, ax in enumerate(axes.flat):\n",
    "    sex = list(tempD2.keys())[ax_i]\n",
    "    tempDF = tempD1[sex]\n",
    "    #Y=X as reference\n",
    "    ax.plot([range_min, range_max], [range_min, range_max], color='black', linestyle=(0, (1, 2)))\n",
    "    #Regplot\n",
    "    sns.regplot(data=tempDF, x=yvar_model, y=yvar, color=tempD2[sex],\n",
    "                scatter=True, fit_reg=True, ci=95, truncate=False, marker='o',\n",
    "                scatter_kws={'alpha':0.2, 'edgecolor':'k', 's':30}, ax=ax)\n",
    "    if ax_i%len(tempD2)==0:\n",
    "        plt.setp(ax, xlabel='', ylabel='Measured '+axis_label)\n",
    "    else:\n",
    "        plt.setp(ax, xlabel='', ylabel='')\n",
    "        plt.setp(ax.get_yticklabels(), visible=False)\n",
    "    ##Annotate Pearson's correlation\n",
    "    pearson_r, pval = stats.pearsonr(tempDF[yvar_model], tempDF[yvar])\n",
    "    r_text = str(Decimal(str(pearson_r)).quantize(Decimal('0.001'), rounding=ROUND_HALF_UP))\n",
    "    if pval==1.0:\n",
    "        pval_text = '1.0'\n",
    "    elif pval==0.0:\n",
    "        pval_text = '0.0'\n",
    "    else:\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",
    "    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='x-small', color='k')\n",
    "    ##Facet label\n",
    "    ax.set_title(sex, {'fontsize':'medium'})\n",
    "    #Save position to generate facet and legend later\n",
    "    if ax_i ==0:\n",
    "        ax0_pos = ax.get_position().bounds\n",
    "    elif ax_i==2:\n",
    "        ax2_pos = ax.get_position().bounds\n",
    "sns.despine()\n",
    "fig.text(x=(ax0_pos[0]+(ax2_pos[0]+ax2_pos[2]))/2, y=ax0_pos[1]-ax0_pos[3]*0.2,#Minor manual adjustment\n",
    "         s='Predicted '+axis_label, fontsize='medium',\n",
    "         verticalalignment='top', horizontalalignment='center', rotation='horizontal')\n",
    "plt.show()\n",
    "\n",
    "#Plot difference b/w sex-specific and sex-mixed models\n",
    "tempD2 = {'Female':'tab:red', 'Male':'tab:blue'}\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "fig, axes = plt.subplots(nrows=1, ncols=len(tempD2), figsize=(6.5, 3), sharex=True, sharey=True,\n",
    "                         gridspec_kw={'hspace':0.1, 'wspace':0.1})\n",
    "for ax_i, ax in enumerate(axes.flat):\n",
    "    sex = list(tempD2.keys())[ax_i]\n",
    "    #Prepare DF\n",
    "    tempS1 = tempD1[sex][yvar_model]\n",
    "    tempS1.name = 'Sex-specific'\n",
    "    tempS2 = tempD1['Both sex'][yvar_model]\n",
    "    tempS2.name = 'Sex-mixed'\n",
    "    tempDF = pd.merge(tempS1, tempS2, left_index=True, right_index=True, how='inner')\n",
    "    #Y=X as reference\n",
    "    ax.plot([range_min, range_max], [range_min, range_max], color='black', linestyle=(0, (1, 2)))\n",
    "    #Regplot\n",
    "    sns.regplot(data=tempDF, x='Sex-specific', y='Sex-mixed', color=tempD2[sex],\n",
    "                scatter=True, fit_reg=True, ci=95, truncate=False, marker='o',\n",
    "                scatter_kws={'alpha':0.2, 'edgecolor':'k', 's':30}, ax=ax)\n",
    "    if ax_i%len(tempD2)==0:\n",
    "        plt.setp(ax, xlabel='', ylabel='Sex-mixed b'+axis_label)\n",
    "    else:\n",
    "        plt.setp(ax, xlabel='', ylabel='')\n",
    "        plt.setp(ax.get_yticklabels(), visible=False)\n",
    "    ##Annotate Pearson's correlation\n",
    "    pearson_r, pval = stats.pearsonr(tempDF['Sex-specific'], tempDF['Sex-mixed'])\n",
    "    r_text = str(Decimal(str(pearson_r)).quantize(Decimal('0.001'), rounding=ROUND_HALF_UP))\n",
    "    if pval==1.0:\n",
    "        pval_text = '1.0'\n",
    "    elif pval==0.0:\n",
    "        pval_text = '0.0'\n",
    "    else:\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",
    "    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='x-small', color='k')\n",
    "    ##Facet label\n",
    "    ax.set_title(sex, {'fontsize':'medium'})\n",
    "    #Save position to generate facet and legend later\n",
    "    if ax_i ==0:\n",
    "        ax0_pos = ax.get_position().bounds\n",
    "    elif ax_i==1:\n",
    "        ax1_pos = ax.get_position().bounds\n",
    "sns.despine()\n",
    "fig.text(x=(ax0_pos[0]+(ax1_pos[0]+ax1_pos[2]))/2, y=ax0_pos[1]-ax0_pos[3]*0.2,#Minor manual adjustment\n",
    "         s='Sex-specific b'+axis_label, fontsize='medium',\n",
    "         verticalalignment='top', horizontalalignment='center', rotation='horizontal')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2-4. Clean beta-coefficient dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempD1 = {'Female':metBMI_F_bcoefs, 'Male':metBMI_M_bcoefs, 'BothSex':metBMI_B_bcoefs}\n",
    "tempD2 = {'Female':metBMI_F_intercept, 'Male':metBMI_M_intercept, 'BothSex':metBMI_B_intercept}\n",
    "yvar_model = 'MetBMI'\n",
    "model_method = 'Ridge'\n",
    "tempD3 = {}\n",
    "for sex in tempD1.keys():\n",
    "    #Combine variables and intercept\n",
    "    tempDF1 = tempD1[sex]\n",
    "    tempDF2 = tempD2[sex]\n",
    "    tempDF = pd.concat([tempDF1, tempDF2], axis=0)\n",
    "    #Summarize\n",
    "    tempL1 = []\n",
    "    tempL2 = []\n",
    "    tempL3 = []\n",
    "    for row_n in tempDF.index.tolist():\n",
    "        tempL1.append(tempDF.loc[row_n].mean())\n",
    "        tempL2.append(tempDF.loc[row_n].std(ddof=1))#Sample standard deviation\n",
    "        tempL3.append((tempDF.loc[row_n]==0.0).astype('int64').sum())\n",
    "    tempDF['Mean'] = tempL1\n",
    "    tempDF['SD'] = tempL2\n",
    "    tempDF['nZeros'] = tempL3\n",
    "    #Save the cleaned DF\n",
    "    fileDir = './ExportData/'\n",
    "    ipynbName = '220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-Ridge_'\n",
    "    fileName = yvar_model+'-'+sex+'-'+model_method+'bcoefs.tsv'\n",
    "    tempDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "    \n",
    "    #Check\n",
    "    tempDF2 = tempDF.loc[tempDF.index.isin(['Intercept'])]#Retrieve as pd.DataFrame\n",
    "    tempDF = tempDF.drop(index=['Intercept'])\n",
    "    print(sex+':')\n",
    "    print(' - Variables:', len(tempDF))\n",
    "    tempDF1 = tempDF.loc[tempDF['nZeros']!=10]\n",
    "    print(' - Variables with non-zero beta-coefficient:', len(tempDF1),\n",
    "          '(', len(tempDF1)/len(tempDF)*100, '%)')\n",
    "    tempDF1 = tempDF.loc[tempDF['nZeros']==0]\n",
    "    print(' - Variables with non-zero beta-coefficient in all 10 models:', len(tempDF1),\n",
    "          '(', len(tempDF1)/len(tempDF)*100, '%)')\n",
    "    tempDF1 = tempDF1.sort_values(by='Mean', ascending=False)\n",
    "    display(tempDF1)\n",
    "    display(tempDF2)#Intercept\n",
    "    print('')\n",
    "    \n",
    "    tempD3[sex] = tempDF\n",
    "#Update for using R2 transition analysis\n",
    "metBMI_F_bcoefs = tempD3['Female']\n",
    "metBMI_M_bcoefs = tempD3['Male']\n",
    "metBMI_B_bcoefs = tempD3['BothSex']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Proteomics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 3-1. Standardization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Note that all preprocessing, including outlier elimination, missingness filtering, imputation, and standardization, should be performed within the cross-validation fold, not across whole dataset, to minimize the potential data leakage. In this study, however, an external dataset is used for “validation\" of the findings (i.e., a testing set for the fitted models). Hence, the robustness of preprocessing is prioritized.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF = protDF\n",
    "tempD1 = {'Female':bmiDF_F, 'Male':bmiDF_M, 'BothSex':bmiDF_B}\n",
    "tempD2 = {}\n",
    "for sex in tempD1.keys():\n",
    "    tempDF1 = tempD1[sex]\n",
    "    tempDF1 = tempDF.loc[tempDF1.index.tolist()]\n",
    "    #Z-score transformation\n",
    "    scaler = StandardScaler(copy=True, with_mean=True, with_std=True)\n",
    "    tempA = scaler.fit_transform(tempDF1)#Column direction\n",
    "    tempDF1 = pd.DataFrame(data=tempA, index=tempDF1.index, columns=tempDF1.columns)\n",
    "    tempD2[sex] = tempDF1\n",
    "    \n",
    "    #Confirmation\n",
    "    print(sex+':', tempDF1.shape)\n",
    "    display(tempDF1.describe())\n",
    "    sns.set(style='ticks', font='Arial', context='notebook')\n",
    "    plt.figure(figsize=(4, 3))\n",
    "    for col_i in range(3):\n",
    "        sns.distplot(tempDF1.iloc[:, col_i], label=tempDF1.columns[col_i])\n",
    "    sns.despine()\n",
    "    plt.xlabel(r'$Z$'+'-score')\n",
    "    plt.ylabel('Density')\n",
    "    plt.legend(bbox_to_anchor=(1, 0.5), loc='center left', borderaxespad=1)\n",
    "    plt.show()\n",
    "    print('')\n",
    "\n",
    "protDF_F = tempD2['Female']\n",
    "protDF_M = tempD2['Male']\n",
    "protDF_B = tempD2['BothSex']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3-2. Ridge with cross-validation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Note that, in contrast to LassoCV and similarly to LinearRegression, the current RidgeCV seems to produce 2D ndarray although target y is 1D.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "#Female model\n",
    "tempDF1 = protDF_F#Standardized independent variables\n",
    "tempDF2 = bmiDF_F#Not-standardized dependent variable and info about testing set\n",
    "yvar = 'BaseBMI'\n",
    "nmodels = len(tempDF2['Testing'].unique())#Reset just in case\n",
    "ncvs = 10\n",
    "model = RidgeCV(alphas=np.geomspace(1, 1000, num=200), fit_intercept=True,\n",
    "                normalize=False, gcv_mode='auto', cv=ncvs)\n",
    "yvar_model = 'BaseProtBMI'\n",
    "\n",
    "#Perform Ridge\n",
    "bcoefDF = pd.DataFrame(index=tempDF1.columns).astype('float64')#For beta-coefficients\n",
    "bcoefDF.index.rename('Variable', inplace=True)\n",
    "interceptDF = pd.DataFrame(index=['Intercept']).astype('float64')#For intercept\n",
    "interceptDF.index.rename('Variable', inplace=True)\n",
    "scoreL = []#For the coefficient of determination R2\n",
    "predictS = pd.Series(name='log_'+yvar_model).astype('float64')#For predictions\n",
    "t_start = time.time()\n",
    "for model_k in range(nmodels):\n",
    "    #Prepare training and testing (hold-out) datasets in model k\n",
    "    model_n = 'Model_'+str(model_k+1).zfill(2)\n",
    "    yDF_train = tempDF2.loc[tempDF2['Testing']!=model_n]\n",
    "    yDF_test = tempDF2.loc[tempDF2['Testing']==model_n]\n",
    "    xDF_train = tempDF1.loc[yDF_train.index.tolist()]\n",
    "    xDF_test = tempDF1.loc[yDF_test.index.tolist()]\n",
    "    #Retrieve the dependent variable\n",
    "    yDF_train = pd.DataFrame(yDF_train['log_'+yvar])#Not Series but DF\n",
    "    yDF_test = pd.DataFrame(yDF_test['log_'+yvar])#Not Series but DF\n",
    "    #Fitting model with cross-validation using training dataset (i.e., internal training/validation datasets)\n",
    "    model.fit(xDF_train, yDF_train)\n",
    "    #Check the penalization amount decided by cross validation\n",
    "    print(model_n+': alpha =', model.alpha_)\n",
    "    #Save parameters\n",
    "    bcoefDF[model_n] = model.coef_.flatten()#weighted vector\n",
    "    interceptDF[model_n] = model.intercept_#Independent term in decision function\n",
    "    #Evaluation with testing (hold-out) dataset\n",
    "    scoreL.append(model.score(xDF_test, yDF_test))\n",
    "    #Prediction for testing dataset using the fitted model k\n",
    "    tempS = pd.Series(model.predict(xDF_test).flatten(),\n",
    "                      index=xDF_test.index, name='log_'+yvar_model)\n",
    "    predictS = pd.concat([predictS, tempS], axis=0)\n",
    "t_elapsed = time.time() - t_start\n",
    "print('Elapsed time for '+str(nmodels)+' models of '+str(ncvs)+'-fold CV Ridge:',\n",
    "      round(t_elapsed//60), 'min', round(t_elapsed%60, 1), 'sec')\n",
    "\n",
    "#Clean prediction DF\n",
    "tempDF = pd.merge(tempDF2[['Testing', 'log_'+yvar]], predictS,\n",
    "                  left_index=True, right_index=True, how='left')\n",
    "##Convert to original scale\n",
    "tempDF[yvar] = np.e**tempDF['log_'+yvar]\n",
    "tempDF[yvar_model] = np.e**tempDF['log_'+yvar_model]\n",
    "display(tempDF)\n",
    "\n",
    "#Save the cleaned prediction DF\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-Ridge_'\n",
    "fileName = yvar_model.replace('Base', '')+'-Female.tsv'\n",
    "tempDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "protBMI_F_bcoefs = bcoefDF\n",
    "protBMI_F_intercept = interceptDF\n",
    "protBMI_F_R2 = scoreL\n",
    "protBMI_F = tempDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "#Male model\n",
    "tempDF1 = protDF_M#Standardized independent variables\n",
    "tempDF2 = bmiDF_M#Not-standardized dependent variable and info about testing set\n",
    "yvar = 'BaseBMI'\n",
    "nmodels = len(tempDF2['Testing'].unique())#Reset just in case\n",
    "ncvs = 10\n",
    "model = RidgeCV(alphas=np.geomspace(1, 1000, num=200), fit_intercept=True,\n",
    "                normalize=False, gcv_mode='auto', cv=ncvs)\n",
    "yvar_model = 'BaseProtBMI'\n",
    "\n",
    "#Perform Ridge\n",
    "bcoefDF = pd.DataFrame(index=tempDF1.columns).astype('float64')#For beta-coefficients\n",
    "bcoefDF.index.rename('Variable', inplace=True)\n",
    "interceptDF = pd.DataFrame(index=['Intercept']).astype('float64')#For intercept\n",
    "interceptDF.index.rename('Variable', inplace=True)\n",
    "scoreL = []#For the coefficient of determination R2\n",
    "predictS = pd.Series(name='log_'+yvar_model).astype('float64')#For predictions\n",
    "t_start = time.time()\n",
    "for model_k in range(nmodels):\n",
    "    #Prepare training and testing (hold-out) datasets in model k\n",
    "    model_n = 'Model_'+str(model_k+1).zfill(2)\n",
    "    yDF_train = tempDF2.loc[tempDF2['Testing']!=model_n]\n",
    "    yDF_test = tempDF2.loc[tempDF2['Testing']==model_n]\n",
    "    xDF_train = tempDF1.loc[yDF_train.index.tolist()]\n",
    "    xDF_test = tempDF1.loc[yDF_test.index.tolist()]\n",
    "    #Retrieve the dependent variable\n",
    "    yDF_train = pd.DataFrame(yDF_train['log_'+yvar])#Not Series but DF\n",
    "    yDF_test = pd.DataFrame(yDF_test['log_'+yvar])#Not Series but DF\n",
    "    #Fitting model with cross-validation using training dataset (i.e., internal training/validation datasets)\n",
    "    model.fit(xDF_train, yDF_train)\n",
    "    #Check the penalization amount decided by cross validation\n",
    "    print(model_n+': alpha =', model.alpha_)\n",
    "    #Save parameters\n",
    "    bcoefDF[model_n] = model.coef_.flatten()#weighted vector\n",
    "    interceptDF[model_n] = model.intercept_#Independent term in decision function\n",
    "    #Evaluation with testing (hold-out) dataset\n",
    "    scoreL.append(model.score(xDF_test, yDF_test))\n",
    "    #Prediction for testing dataset using the fitted model k\n",
    "    tempS = pd.Series(model.predict(xDF_test).flatten(),\n",
    "                      index=xDF_test.index, name='log_'+yvar_model)\n",
    "    predictS = pd.concat([predictS, tempS], axis=0)\n",
    "t_elapsed = time.time() - t_start\n",
    "print('Elapsed time for '+str(nmodels)+' models of '+str(ncvs)+'-fold CV Ridge:',\n",
    "      round(t_elapsed//60), 'min', round(t_elapsed%60, 1), 'sec')\n",
    "\n",
    "#Clean prediction DF\n",
    "tempDF = pd.merge(tempDF2[['Testing', 'log_'+yvar]], predictS,\n",
    "                  left_index=True, right_index=True, how='left')\n",
    "##Convert to original scale\n",
    "tempDF[yvar] = np.e**tempDF['log_'+yvar]\n",
    "tempDF[yvar_model] = np.e**tempDF['log_'+yvar_model]\n",
    "display(tempDF)\n",
    "\n",
    "#Save the cleaned prediction DF\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-Ridge_'\n",
    "fileName = yvar_model.replace('Base', '')+'-Male.tsv'\n",
    "tempDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "protBMI_M_bcoefs = bcoefDF\n",
    "protBMI_M_intercept = interceptDF\n",
    "protBMI_M_R2 = scoreL\n",
    "protBMI_M = tempDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Both sex model\n",
    "tempDF1 = protDF_B#Standardized independent variables\n",
    "tempDF2 = bmiDF_B#Not-standardized dependent variable and info about testing set\n",
    "yvar = 'BaseBMI'\n",
    "nmodels = len(tempDF2['Testing'].unique())#Reset just in case\n",
    "ncvs = 10\n",
    "model = RidgeCV(alphas=np.geomspace(1, 1000, num=200), fit_intercept=True,\n",
    "                normalize=False, gcv_mode='auto', cv=ncvs)\n",
    "yvar_model = 'BaseProtBMI'\n",
    "\n",
    "#Perform Ridge\n",
    "bcoefDF = pd.DataFrame(index=tempDF1.columns).astype('float64')#For beta-coefficients\n",
    "bcoefDF.index.rename('Variable', inplace=True)\n",
    "interceptDF = pd.DataFrame(index=['Intercept']).astype('float64')#For intercept\n",
    "interceptDF.index.rename('Variable', inplace=True)\n",
    "scoreL = []#For the coefficient of determination R2\n",
    "predictS = pd.Series(name='log_'+yvar_model).astype('float64')#For predictions\n",
    "t_start = time.time()\n",
    "for model_k in range(nmodels):\n",
    "    #Prepare training and testing (hold-out) datasets in model k\n",
    "    model_n = 'Model_'+str(model_k+1).zfill(2)\n",
    "    yDF_train = tempDF2.loc[tempDF2['Testing']!=model_n]\n",
    "    yDF_test = tempDF2.loc[tempDF2['Testing']==model_n]\n",
    "    xDF_train = tempDF1.loc[yDF_train.index.tolist()]\n",
    "    xDF_test = tempDF1.loc[yDF_test.index.tolist()]\n",
    "    #Retrieve the dependent variable\n",
    "    yDF_train = pd.DataFrame(yDF_train['log_'+yvar])#Not Series but DF\n",
    "    yDF_test = pd.DataFrame(yDF_test['log_'+yvar])#Not Series but DF\n",
    "    #Fitting model with cross-validation using training dataset (i.e., internal training/validation datasets)\n",
    "    model.fit(xDF_train, yDF_train)\n",
    "    #Check the penalization amount decided by cross validation\n",
    "    print(model_n+': alpha =', model.alpha_)\n",
    "    #Save parameters\n",
    "    bcoefDF[model_n] = model.coef_.flatten()#weighted vector\n",
    "    interceptDF[model_n] = model.intercept_#Independent term in decision function\n",
    "    #Evaluation with testing (hold-out) dataset\n",
    "    scoreL.append(model.score(xDF_test, yDF_test))\n",
    "    #Prediction for testing dataset using the fitted model k\n",
    "    tempS = pd.Series(model.predict(xDF_test).flatten(),\n",
    "                      index=xDF_test.index, name='log_'+yvar_model)\n",
    "    predictS = pd.concat([predictS, tempS], axis=0)\n",
    "t_elapsed = time.time() - t_start\n",
    "print('Elapsed time for '+str(nmodels)+' models of '+str(ncvs)+'-fold CV Ridge:',\n",
    "      round(t_elapsed//60), 'min', round(t_elapsed%60, 1), 'sec')\n",
    "\n",
    "#Clean prediction DF\n",
    "tempDF = pd.merge(tempDF2[['Testing', 'log_'+yvar]], predictS,\n",
    "                  left_index=True, right_index=True, how='left')\n",
    "##Convert to original scale\n",
    "tempDF[yvar] = np.e**tempDF['log_'+yvar]\n",
    "tempDF[yvar_model] = np.e**tempDF['log_'+yvar_model]\n",
    "display(tempDF)\n",
    "\n",
    "#Save the cleaned prediction DF\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-Ridge_'\n",
    "fileName = yvar_model.replace('Base', '')+'-BothSex.tsv'\n",
    "tempDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "protBMI_B_bcoefs = bcoefDF\n",
    "protBMI_B_intercept = interceptDF\n",
    "protBMI_B_R2 = scoreL\n",
    "protBMI_B = tempDF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3-3. Prediction accuracy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Summary\n",
    "tempD1 = {'Female':protBMI_F_R2, 'Male':protBMI_M_R2, 'Both sex':protBMI_B_R2}\n",
    "tempD2 = {'Female':protBMI_F, 'Male':protBMI_M, 'Both sex':protBMI_B}\n",
    "yvar = 'BaseBMI'\n",
    "yvar_model = 'BaseProtBMI'\n",
    "\n",
    "for sex in tempD1.keys():\n",
    "    tempL = tempD1[sex]\n",
    "    print(sex+' 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",
    "    tempDF = tempD2[sex]\n",
    "    print(' - Observed vs. predicted log_'+yvar+': (Pearson\\'s r, P) =',\n",
    "          stats.pearsonr(tempDF['log_'+yvar], tempDF['log_'+yvar_model]))\n",
    "    print(' - Observed vs. predicted '+yvar+': (Pearson\\'s r, P) =',\n",
    "          stats.pearsonr(tempDF[yvar], tempDF[yvar_model]))\n",
    "    display(tempDF.describe())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Check difference between sex-specific and sex-mixed models for now. Note that this is a rough comparison because the sample size is different.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Prepare DF\n",
    "tempDF = pd.DataFrame({'Female':protBMI_F_R2, 'Male':protBMI_M_R2, 'Both sex':protBMI_B_R2})\n",
    "tempDF = tempDF.melt(var_name='Cohort', value_name='R2', value_vars=tempDF.columns.tolist())\n",
    "\n",
    "#Plot\n",
    "sns.set(style='ticks', font='Arial', context='notebook')\n",
    "plt.figure(figsize=(3, 1))\n",
    "sns.barplot(data=tempDF, y='Cohort', x='R2', hue='Cohort', palette='Set1', dodge=False, edgecolor='black',\n",
    "            ci=95, capsize=0.4, errwidth=1.5, errcolor='black')\n",
    "sns.stripplot(data=tempDF, y='Cohort', x='R2', hue='Cohort', dodge=False, size=8, edgecolor='black',\n",
    "              linewidth=1, alpha=0.4, palette={'Female':'gray', 'Male':'gray', 'Both sex':'gray'})\n",
    "sns.despine()\n",
    "plt.xlabel('Out-of-sample '+r'$R^2$'+' in 10 Ridge models\\n[mean with 95% CI]')\n",
    "plt.ylabel('')\n",
    "plt.legend('', frameon=False)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempD1 = {'Female':protBMI_F, 'Male':protBMI_M, 'Both sex':protBMI_B}\n",
    "yvar = 'BaseBMI'\n",
    "yvar_model = 'BaseProtBMI'\n",
    "range_min = np.min([df[var].min() for df in tempD1.values() for var in [yvar, yvar_model]])\n",
    "range_max = np.max([df[var].max() for df in tempD1.values() for var in [yvar, yvar_model]])\n",
    "axis_label = 'BMI [kg m'+r'$^{-2}$'+']'\n",
    "\n",
    "#Plot measured vs. predicted per model\n",
    "tempD2 = {'Female':'tab:red', 'Male':'tab:blue', 'Both sex':'tab:green'}\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "fig, axes = plt.subplots(nrows=1, ncols=len(tempD2), figsize=(10, 3), sharex=True, sharey=True,\n",
    "                         gridspec_kw={'hspace':0.1, 'wspace':0.1})\n",
    "for ax_i, ax in enumerate(axes.flat):\n",
    "    sex = list(tempD2.keys())[ax_i]\n",
    "    tempDF = tempD1[sex]\n",
    "    #Y=X as reference\n",
    "    ax.plot([range_min, range_max], [range_min, range_max], color='black', linestyle=(0, (1, 2)))\n",
    "    #Regplot\n",
    "    sns.regplot(data=tempDF, x=yvar_model, y=yvar, color=tempD2[sex],\n",
    "                scatter=True, fit_reg=True, ci=95, truncate=False, marker='o',\n",
    "                scatter_kws={'alpha':0.2, 'edgecolor':'k', 's':30}, ax=ax)\n",
    "    if ax_i%len(tempD2)==0:\n",
    "        plt.setp(ax, xlabel='', ylabel='Measured '+axis_label)\n",
    "    else:\n",
    "        plt.setp(ax, xlabel='', ylabel='')\n",
    "        plt.setp(ax.get_yticklabels(), visible=False)\n",
    "    ##Annotate Pearson's correlation\n",
    "    pearson_r, pval = stats.pearsonr(tempDF[yvar_model], tempDF[yvar])\n",
    "    r_text = str(Decimal(str(pearson_r)).quantize(Decimal('0.001'), rounding=ROUND_HALF_UP))\n",
    "    if pval==1.0:\n",
    "        pval_text = '1.0'\n",
    "    elif pval==0.0:\n",
    "        pval_text = '0.0'\n",
    "    else:\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",
    "    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='x-small', color='k')\n",
    "    ##Facet label\n",
    "    ax.set_title(sex, {'fontsize':'medium'})\n",
    "    #Save position to generate facet and legend later\n",
    "    if ax_i ==0:\n",
    "        ax0_pos = ax.get_position().bounds\n",
    "    elif ax_i==2:\n",
    "        ax2_pos = ax.get_position().bounds\n",
    "sns.despine()\n",
    "fig.text(x=(ax0_pos[0]+(ax2_pos[0]+ax2_pos[2]))/2, y=ax0_pos[1]-ax0_pos[3]*0.2,#Minor manual adjustment\n",
    "         s='Predicted '+axis_label, fontsize='medium',\n",
    "         verticalalignment='top', horizontalalignment='center', rotation='horizontal')\n",
    "plt.show()\n",
    "\n",
    "#Plot difference b/w sex-specific and sex-mixed models\n",
    "tempD2 = {'Female':'tab:red', 'Male':'tab:blue'}\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "fig, axes = plt.subplots(nrows=1, ncols=len(tempD2), figsize=(6.5, 3), sharex=True, sharey=True,\n",
    "                         gridspec_kw={'hspace':0.1, 'wspace':0.1})\n",
    "for ax_i, ax in enumerate(axes.flat):\n",
    "    sex = list(tempD2.keys())[ax_i]\n",
    "    #Prepare DF\n",
    "    tempS1 = tempD1[sex][yvar_model]\n",
    "    tempS1.name = 'Sex-specific'\n",
    "    tempS2 = tempD1['Both sex'][yvar_model]\n",
    "    tempS2.name = 'Sex-mixed'\n",
    "    tempDF = pd.merge(tempS1, tempS2, left_index=True, right_index=True, how='inner')\n",
    "    #Y=X as reference\n",
    "    ax.plot([range_min, range_max], [range_min, range_max], color='black', linestyle=(0, (1, 2)))\n",
    "    #Regplot\n",
    "    sns.regplot(data=tempDF, x='Sex-specific', y='Sex-mixed', color=tempD2[sex],\n",
    "                scatter=True, fit_reg=True, ci=95, truncate=False, marker='o',\n",
    "                scatter_kws={'alpha':0.2, 'edgecolor':'k', 's':30}, ax=ax)\n",
    "    if ax_i%len(tempD2)==0:\n",
    "        plt.setp(ax, xlabel='', ylabel='Sex-mixed b'+axis_label)\n",
    "    else:\n",
    "        plt.setp(ax, xlabel='', ylabel='')\n",
    "        plt.setp(ax.get_yticklabels(), visible=False)\n",
    "    ##Annotate Pearson's correlation\n",
    "    pearson_r, pval = stats.pearsonr(tempDF['Sex-specific'], tempDF['Sex-mixed'])\n",
    "    r_text = str(Decimal(str(pearson_r)).quantize(Decimal('0.001'), rounding=ROUND_HALF_UP))\n",
    "    if pval==1.0:\n",
    "        pval_text = '1.0'\n",
    "    elif pval==0.0:\n",
    "        pval_text = '0.0'\n",
    "    else:\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",
    "    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='x-small', color='k')\n",
    "    ##Facet label\n",
    "    ax.set_title(sex, {'fontsize':'medium'})\n",
    "    #Save position to generate facet and legend later\n",
    "    if ax_i ==0:\n",
    "        ax0_pos = ax.get_position().bounds\n",
    "    elif ax_i==1:\n",
    "        ax1_pos = ax.get_position().bounds\n",
    "sns.despine()\n",
    "fig.text(x=(ax0_pos[0]+(ax1_pos[0]+ax1_pos[2]))/2, y=ax0_pos[1]-ax0_pos[3]*0.2,#Minor manual adjustment\n",
    "         s='Sex-specific b'+axis_label, fontsize='medium',\n",
    "         verticalalignment='top', horizontalalignment='center', rotation='horizontal')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3-4. Clean beta-coefficient dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempD1 = {'Female':protBMI_F_bcoefs, 'Male':protBMI_M_bcoefs, 'BothSex':protBMI_B_bcoefs}\n",
    "tempD2 = {'Female':protBMI_F_intercept, 'Male':protBMI_M_intercept, 'BothSex':protBMI_B_intercept}\n",
    "yvar_model = 'ProtBMI'\n",
    "model_method = 'Ridge'\n",
    "tempD3 = {}\n",
    "for sex in tempD1.keys():\n",
    "    #Combine variables and intercept\n",
    "    tempDF1 = tempD1[sex]\n",
    "    tempDF2 = tempD2[sex]\n",
    "    tempDF = pd.concat([tempDF1, tempDF2], axis=0)\n",
    "    #Summarize\n",
    "    tempL1 = []\n",
    "    tempL2 = []\n",
    "    tempL3 = []\n",
    "    for row_n in tempDF.index.tolist():\n",
    "        tempL1.append(tempDF.loc[row_n].mean())\n",
    "        tempL2.append(tempDF.loc[row_n].std(ddof=1))#Sample standard deviation\n",
    "        tempL3.append((tempDF.loc[row_n]==0.0).astype('int64').sum())\n",
    "    tempDF['Mean'] = tempL1\n",
    "    tempDF['SD'] = tempL2\n",
    "    tempDF['nZeros'] = tempL3\n",
    "    #Save the cleaned DF\n",
    "    fileDir = './ExportData/'\n",
    "    ipynbName = '220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-Ridge_'\n",
    "    fileName = yvar_model+'-'+sex+'-'+model_method+'bcoefs.tsv'\n",
    "    tempDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "    \n",
    "    #Check\n",
    "    tempDF2 = tempDF.loc[tempDF.index.isin(['Intercept'])]#Retrieve as pd.DataFrame\n",
    "    tempDF = tempDF.drop(index=['Intercept'])\n",
    "    print(sex+':')\n",
    "    print(' - Variables:', len(tempDF))\n",
    "    tempDF1 = tempDF.loc[tempDF['nZeros']!=10]\n",
    "    print(' - Variables with non-zero beta-coefficient:', len(tempDF1),\n",
    "          '(', len(tempDF1)/len(tempDF)*100, '%)')\n",
    "    tempDF1 = tempDF.loc[tempDF['nZeros']==0]\n",
    "    print(' - Variables with non-zero beta-coefficient in all 10 models:', len(tempDF1),\n",
    "          '(', len(tempDF1)/len(tempDF)*100, '%)')\n",
    "    tempDF1 = tempDF1.sort_values(by='Mean', ascending=False)\n",
    "    display(tempDF1)\n",
    "    display(tempDF2)#Intercept\n",
    "    print('')\n",
    "    \n",
    "    tempD3[sex] = tempDF\n",
    "#Update for using R2 transition analysis\n",
    "protBMI_F_bcoefs = tempD3['Female']\n",
    "protBMI_M_bcoefs = tempD3['Male']\n",
    "protBMI_B_bcoefs = tempD3['BothSex']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Clinical labs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 4-1. Standardization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Note that all preprocessing, including outlier elimination, missingness filtering, imputation, and standardization, should be performed within the cross-validation fold, not across whole dataset, to minimize the potential data leakage. In this study, however, an external dataset is used for “validation\" of the findings (i.e., a testing set for the fitted models). Hence, the robustness of preprocessing is prioritized.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF = chemDF\n",
    "tempD1 = {'Female':bmiDF_F, 'Male':bmiDF_M, 'BothSex':bmiDF_B}\n",
    "tempD2 = {}\n",
    "for sex in tempD1.keys():\n",
    "    tempDF1 = tempD1[sex]\n",
    "    tempDF1 = tempDF.loc[tempDF1.index.tolist()]\n",
    "    #Z-score transformation\n",
    "    scaler = StandardScaler(copy=True, with_mean=True, with_std=True)\n",
    "    tempA = scaler.fit_transform(tempDF1)#Column direction\n",
    "    tempDF1 = pd.DataFrame(data=tempA, index=tempDF1.index, columns=tempDF1.columns)\n",
    "    tempD2[sex] = tempDF1\n",
    "    \n",
    "    #Confirmation\n",
    "    print(sex+':', tempDF1.shape)\n",
    "    display(tempDF1.describe())\n",
    "    sns.set(style='ticks', font='Arial', context='notebook')\n",
    "    plt.figure(figsize=(4, 3))\n",
    "    for col_i in range(3):\n",
    "        sns.distplot(tempDF1.iloc[:, col_i], label=tempDF1.columns[col_i])\n",
    "    sns.despine()\n",
    "    plt.xlabel(r'$Z$'+'-score')\n",
    "    plt.ylabel('Density')\n",
    "    plt.legend(bbox_to_anchor=(1, 0.5), loc='center left', borderaxespad=1)\n",
    "    plt.show()\n",
    "    print('')\n",
    "\n",
    "chemDF_F = tempD2['Female']\n",
    "chemDF_M = tempD2['Male']\n",
    "chemDF_B = tempD2['BothSex']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4-2. Ridge with cross-validation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Note that, in contrast to LassoCV and similarly to LinearRegression, the current RidgeCV seems to produce 2D ndarray although target y is 1D.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Female model\n",
    "tempDF1 = chemDF_F#Standardized independent variables\n",
    "tempDF2 = bmiDF_F#Not-standardized dependent variable and info about testing set\n",
    "yvar = 'BaseBMI'\n",
    "nmodels = len(tempDF2['Testing'].unique())#Reset just in case\n",
    "ncvs = 10\n",
    "model = RidgeCV(alphas=np.geomspace(1, 1000, num=200), fit_intercept=True,\n",
    "                normalize=False, gcv_mode='auto', cv=ncvs)\n",
    "yvar_model = 'BaseChemBMI'\n",
    "\n",
    "#Perform Ridge\n",
    "bcoefDF = pd.DataFrame(index=tempDF1.columns).astype('float64')#For beta-coefficients\n",
    "bcoefDF.index.rename('Variable', inplace=True)\n",
    "interceptDF = pd.DataFrame(index=['Intercept']).astype('float64')#For intercept\n",
    "interceptDF.index.rename('Variable', inplace=True)\n",
    "scoreL = []#For the coefficient of determination R2\n",
    "predictS = pd.Series(name='log_'+yvar_model).astype('float64')#For predictions\n",
    "t_start = time.time()\n",
    "for model_k in range(nmodels):\n",
    "    #Prepare training and testing (hold-out) datasets in model k\n",
    "    model_n = 'Model_'+str(model_k+1).zfill(2)\n",
    "    yDF_train = tempDF2.loc[tempDF2['Testing']!=model_n]\n",
    "    yDF_test = tempDF2.loc[tempDF2['Testing']==model_n]\n",
    "    xDF_train = tempDF1.loc[yDF_train.index.tolist()]\n",
    "    xDF_test = tempDF1.loc[yDF_test.index.tolist()]\n",
    "    #Retrieve the dependent variable\n",
    "    yDF_train = pd.DataFrame(yDF_train['log_'+yvar])#Not Series but DF\n",
    "    yDF_test = pd.DataFrame(yDF_test['log_'+yvar])#Not Series but DF\n",
    "    #Fitting model with cross-validation using training dataset (i.e., internal training/validation datasets)\n",
    "    model.fit(xDF_train, yDF_train)\n",
    "    #Check the penalization amount decided by cross validation\n",
    "    print(model_n+': alpha =', model.alpha_)\n",
    "    #Save parameters\n",
    "    bcoefDF[model_n] = model.coef_.flatten()#weighted vector\n",
    "    interceptDF[model_n] = model.intercept_#Independent term in decision function\n",
    "    #Evaluation with testing (hold-out) dataset\n",
    "    scoreL.append(model.score(xDF_test, yDF_test))\n",
    "    #Prediction for testing dataset using the fitted model k\n",
    "    tempS = pd.Series(model.predict(xDF_test).flatten(),\n",
    "                      index=xDF_test.index, name='log_'+yvar_model)\n",
    "    predictS = pd.concat([predictS, tempS], axis=0)\n",
    "t_elapsed = time.time() - t_start\n",
    "print('Elapsed time for '+str(nmodels)+' models of '+str(ncvs)+'-fold CV Ridge:',\n",
    "      round(t_elapsed//60), 'min', round(t_elapsed%60, 1), 'sec')\n",
    "\n",
    "#Clean prediction DF\n",
    "tempDF = pd.merge(tempDF2[['Testing', 'log_'+yvar]], predictS,\n",
    "                  left_index=True, right_index=True, how='left')\n",
    "##Convert to original scale\n",
    "tempDF[yvar] = np.e**tempDF['log_'+yvar]\n",
    "tempDF[yvar_model] = np.e**tempDF['log_'+yvar_model]\n",
    "display(tempDF)\n",
    "\n",
    "#Save the cleaned prediction DF\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-Ridge_'\n",
    "fileName = yvar_model.replace('Base', '')+'-Female.tsv'\n",
    "tempDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "chemBMI_F_bcoefs = bcoefDF\n",
    "chemBMI_F_intercept = interceptDF\n",
    "chemBMI_F_R2 = scoreL\n",
    "chemBMI_F = tempDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Male model\n",
    "tempDF1 = chemDF_M#Standardized independent variables\n",
    "tempDF2 = bmiDF_M#Not-standardized dependent variable and info about testing set\n",
    "yvar = 'BaseBMI'\n",
    "nmodels = len(tempDF2['Testing'].unique())#Reset just in case\n",
    "ncvs = 10\n",
    "model = RidgeCV(alphas=np.geomspace(1, 1000, num=200), fit_intercept=True,\n",
    "                normalize=False, gcv_mode='auto', cv=ncvs)\n",
    "yvar_model = 'BaseChemBMI'\n",
    "\n",
    "#Perform Ridge\n",
    "bcoefDF = pd.DataFrame(index=tempDF1.columns).astype('float64')#For beta-coefficients\n",
    "bcoefDF.index.rename('Variable', inplace=True)\n",
    "interceptDF = pd.DataFrame(index=['Intercept']).astype('float64')#For intercept\n",
    "interceptDF.index.rename('Variable', inplace=True)\n",
    "scoreL = []#For the coefficient of determination R2\n",
    "predictS = pd.Series(name='log_'+yvar_model).astype('float64')#For predictions\n",
    "t_start = time.time()\n",
    "for model_k in range(nmodels):\n",
    "    #Prepare training and testing (hold-out) datasets in model k\n",
    "    model_n = 'Model_'+str(model_k+1).zfill(2)\n",
    "    yDF_train = tempDF2.loc[tempDF2['Testing']!=model_n]\n",
    "    yDF_test = tempDF2.loc[tempDF2['Testing']==model_n]\n",
    "    xDF_train = tempDF1.loc[yDF_train.index.tolist()]\n",
    "    xDF_test = tempDF1.loc[yDF_test.index.tolist()]\n",
    "    #Retrieve the dependent variable\n",
    "    yDF_train = pd.DataFrame(yDF_train['log_'+yvar])#Not Series but DF\n",
    "    yDF_test = pd.DataFrame(yDF_test['log_'+yvar])#Not Series but DF\n",
    "    #Fitting model with cross-validation using training dataset (i.e., internal training/validation datasets)\n",
    "    model.fit(xDF_train, yDF_train)\n",
    "    #Check the penalization amount decided by cross validation\n",
    "    print(model_n+': alpha =', model.alpha_)\n",
    "    #Save parameters\n",
    "    bcoefDF[model_n] = model.coef_.flatten()#weighted vector\n",
    "    interceptDF[model_n] = model.intercept_#Independent term in decision function\n",
    "    #Evaluation with testing (hold-out) dataset\n",
    "    scoreL.append(model.score(xDF_test, yDF_test))\n",
    "    #Prediction for testing dataset using the fitted model k\n",
    "    tempS = pd.Series(model.predict(xDF_test).flatten(),\n",
    "                      index=xDF_test.index, name='log_'+yvar_model)\n",
    "    predictS = pd.concat([predictS, tempS], axis=0)\n",
    "t_elapsed = time.time() - t_start\n",
    "print('Elapsed time for '+str(nmodels)+' models of '+str(ncvs)+'-fold CV Ridge:',\n",
    "      round(t_elapsed//60), 'min', round(t_elapsed%60, 1), 'sec')\n",
    "\n",
    "#Clean prediction DF\n",
    "tempDF = pd.merge(tempDF2[['Testing', 'log_'+yvar]], predictS,\n",
    "                  left_index=True, right_index=True, how='left')\n",
    "##Convert to original scale\n",
    "tempDF[yvar] = np.e**tempDF['log_'+yvar]\n",
    "tempDF[yvar_model] = np.e**tempDF['log_'+yvar_model]\n",
    "display(tempDF)\n",
    "\n",
    "#Save the cleaned prediction DF\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-Ridge_'\n",
    "fileName = yvar_model.replace('Base', '')+'-Male.tsv'\n",
    "tempDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "chemBMI_M_bcoefs = bcoefDF\n",
    "chemBMI_M_intercept = interceptDF\n",
    "chemBMI_M_R2 = scoreL\n",
    "chemBMI_M = tempDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Both sex model\n",
    "tempDF1 = chemDF_B#Standardized independent variables\n",
    "tempDF2 = bmiDF_B#Not-standardized dependent variable and info about testing set\n",
    "yvar = 'BaseBMI'\n",
    "nmodels = len(tempDF2['Testing'].unique())#Reset just in case\n",
    "ncvs = 10\n",
    "model = RidgeCV(alphas=np.geomspace(1, 1000, num=200), fit_intercept=True,\n",
    "                normalize=False, gcv_mode='auto', cv=ncvs)\n",
    "yvar_model = 'BaseChemBMI'\n",
    "\n",
    "#Perform Ridge\n",
    "bcoefDF = pd.DataFrame(index=tempDF1.columns).astype('float64')#For beta-coefficients\n",
    "bcoefDF.index.rename('Variable', inplace=True)\n",
    "interceptDF = pd.DataFrame(index=['Intercept']).astype('float64')#For intercept\n",
    "interceptDF.index.rename('Variable', inplace=True)\n",
    "scoreL = []#For the coefficient of determination R2\n",
    "predictS = pd.Series(name='log_'+yvar_model).astype('float64')#For predictions\n",
    "t_start = time.time()\n",
    "for model_k in range(nmodels):\n",
    "    #Prepare training and testing (hold-out) datasets in model k\n",
    "    model_n = 'Model_'+str(model_k+1).zfill(2)\n",
    "    yDF_train = tempDF2.loc[tempDF2['Testing']!=model_n]\n",
    "    yDF_test = tempDF2.loc[tempDF2['Testing']==model_n]\n",
    "    xDF_train = tempDF1.loc[yDF_train.index.tolist()]\n",
    "    xDF_test = tempDF1.loc[yDF_test.index.tolist()]\n",
    "    #Retrieve the dependent variable\n",
    "    yDF_train = pd.DataFrame(yDF_train['log_'+yvar])#Not Series but DF\n",
    "    yDF_test = pd.DataFrame(yDF_test['log_'+yvar])#Not Series but DF\n",
    "    #Fitting model with cross-validation using training dataset (i.e., internal training/validation datasets)\n",
    "    model.fit(xDF_train, yDF_train)\n",
    "    #Check the penalization amount decided by cross validation\n",
    "    print(model_n+': alpha =', model.alpha_)\n",
    "    #Save parameters\n",
    "    bcoefDF[model_n] = model.coef_.flatten()#weighted vector\n",
    "    interceptDF[model_n] = model.intercept_#Independent term in decision function\n",
    "    #Evaluation with testing (hold-out) dataset\n",
    "    scoreL.append(model.score(xDF_test, yDF_test))\n",
    "    #Prediction for testing dataset using the fitted model k\n",
    "    tempS = pd.Series(model.predict(xDF_test).flatten(),\n",
    "                      index=xDF_test.index, name='log_'+yvar_model)\n",
    "    predictS = pd.concat([predictS, tempS], axis=0)\n",
    "t_elapsed = time.time() - t_start\n",
    "print('Elapsed time for '+str(nmodels)+' models of '+str(ncvs)+'-fold CV Ridge:',\n",
    "      round(t_elapsed//60), 'min', round(t_elapsed%60, 1), 'sec')\n",
    "\n",
    "#Clean prediction DF\n",
    "tempDF = pd.merge(tempDF2[['Testing', 'log_'+yvar]], predictS,\n",
    "                  left_index=True, right_index=True, how='left')\n",
    "##Convert to original scale\n",
    "tempDF[yvar] = np.e**tempDF['log_'+yvar]\n",
    "tempDF[yvar_model] = np.e**tempDF['log_'+yvar_model]\n",
    "display(tempDF)\n",
    "\n",
    "#Save the cleaned prediction DF\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-Ridge_'\n",
    "fileName = yvar_model.replace('Base', '')+'-BothSex.tsv'\n",
    "tempDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "chemBMI_B_bcoefs = bcoefDF\n",
    "chemBMI_B_intercept = interceptDF\n",
    "chemBMI_B_R2 = scoreL\n",
    "chemBMI_B = tempDF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4-3. Prediction accuracy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Summary\n",
    "tempD1 = {'Female':chemBMI_F_R2, 'Male':chemBMI_M_R2, 'Both sex':chemBMI_B_R2}\n",
    "tempD2 = {'Female':chemBMI_F, 'Male':chemBMI_M, 'Both sex':chemBMI_B}\n",
    "yvar = 'BaseBMI'\n",
    "yvar_model = 'BaseChemBMI'\n",
    "\n",
    "for sex in tempD1.keys():\n",
    "    tempL = tempD1[sex]\n",
    "    print(sex+' 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",
    "    tempDF = tempD2[sex]\n",
    "    print(' - Observed vs. predicted log_'+yvar+': (Pearson\\'s r, P) =',\n",
    "          stats.pearsonr(tempDF['log_'+yvar], tempDF['log_'+yvar_model]))\n",
    "    print(' - Observed vs. predicted '+yvar+': (Pearson\\'s r, P) =',\n",
    "          stats.pearsonr(tempDF[yvar], tempDF[yvar_model]))\n",
    "    display(tempDF.describe())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Check difference between sex-specific and sex-mixed models for now. Note that this is a rough comparison because the sample size is different.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Prepare DF\n",
    "tempDF = pd.DataFrame({'Female':chemBMI_F_R2, 'Male':chemBMI_M_R2, 'Both sex':chemBMI_B_R2})\n",
    "tempDF = tempDF.melt(var_name='Cohort', value_name='R2', value_vars=tempDF.columns.tolist())\n",
    "\n",
    "#Plot\n",
    "sns.set(style='ticks', font='Arial', context='notebook')\n",
    "plt.figure(figsize=(3, 1))\n",
    "sns.barplot(data=tempDF, y='Cohort', x='R2', hue='Cohort', palette='Set1', dodge=False, edgecolor='black',\n",
    "            ci=95, capsize=0.4, errwidth=1.5, errcolor='black')\n",
    "sns.stripplot(data=tempDF, y='Cohort', x='R2', hue='Cohort', dodge=False, size=8, edgecolor='black',\n",
    "              linewidth=1, alpha=0.4, palette={'Female':'gray', 'Male':'gray', 'Both sex':'gray'})\n",
    "sns.despine()\n",
    "plt.xlabel('Out-of-sample '+r'$R^2$'+' in 10 Ridge models\\n[mean with 95% CI]')\n",
    "plt.ylabel('')\n",
    "plt.legend('', frameon=False)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempD1 = {'Female':chemBMI_F, 'Male':chemBMI_M, 'Both sex':chemBMI_B}\n",
    "yvar = 'BaseBMI'\n",
    "yvar_model = 'BaseChemBMI'\n",
    "range_min = np.min([df[var].min() for df in tempD1.values() for var in [yvar, yvar_model]])\n",
    "range_max = np.max([df[var].max() for df in tempD1.values() for var in [yvar, yvar_model]])\n",
    "axis_label = 'BMI [kg m'+r'$^{-2}$'+']'\n",
    "\n",
    "#Plot measured vs. predicted per model\n",
    "tempD2 = {'Female':'tab:red', 'Male':'tab:blue', 'Both sex':'tab:green'}\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "fig, axes = plt.subplots(nrows=1, ncols=len(tempD2), figsize=(10, 3), sharex=True, sharey=True,\n",
    "                         gridspec_kw={'hspace':0.1, 'wspace':0.1})\n",
    "for ax_i, ax in enumerate(axes.flat):\n",
    "    sex = list(tempD2.keys())[ax_i]\n",
    "    tempDF = tempD1[sex]\n",
    "    #Y=X as reference\n",
    "    ax.plot([range_min, range_max], [range_min, range_max], color='black', linestyle=(0, (1, 2)))\n",
    "    #Regplot\n",
    "    sns.regplot(data=tempDF, x=yvar_model, y=yvar, color=tempD2[sex],\n",
    "                scatter=True, fit_reg=True, ci=95, truncate=False, marker='o',\n",
    "                scatter_kws={'alpha':0.2, 'edgecolor':'k', 's':30}, ax=ax)\n",
    "    if ax_i%len(tempD2)==0:\n",
    "        plt.setp(ax, xlabel='', ylabel='Measured '+axis_label)\n",
    "    else:\n",
    "        plt.setp(ax, xlabel='', ylabel='')\n",
    "        plt.setp(ax.get_yticklabels(), visible=False)\n",
    "    ##Annotate Pearson's correlation\n",
    "    pearson_r, pval = stats.pearsonr(tempDF[yvar_model], tempDF[yvar])\n",
    "    r_text = str(Decimal(str(pearson_r)).quantize(Decimal('0.001'), rounding=ROUND_HALF_UP))\n",
    "    if pval==1.0:\n",
    "        pval_text = '1.0'\n",
    "    elif pval==0.0:\n",
    "        pval_text = '0.0'\n",
    "    else:\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",
    "    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='x-small', color='k')\n",
    "    ##Facet label\n",
    "    ax.set_title(sex, {'fontsize':'medium'})\n",
    "    #Save position to generate facet and legend later\n",
    "    if ax_i ==0:\n",
    "        ax0_pos = ax.get_position().bounds\n",
    "    elif ax_i==2:\n",
    "        ax2_pos = ax.get_position().bounds\n",
    "sns.despine()\n",
    "fig.text(x=(ax0_pos[0]+(ax2_pos[0]+ax2_pos[2]))/2, y=ax0_pos[1]-ax0_pos[3]*0.2,#Minor manual adjustment\n",
    "         s='Predicted '+axis_label, fontsize='medium',\n",
    "         verticalalignment='top', horizontalalignment='center', rotation='horizontal')\n",
    "plt.show()\n",
    "\n",
    "#Plot difference b/w sex-specific and sex-mixed models\n",
    "tempD2 = {'Female':'tab:red', 'Male':'tab:blue'}\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "fig, axes = plt.subplots(nrows=1, ncols=len(tempD2), figsize=(6.5, 3), sharex=True, sharey=True,\n",
    "                         gridspec_kw={'hspace':0.1, 'wspace':0.1})\n",
    "for ax_i, ax in enumerate(axes.flat):\n",
    "    sex = list(tempD2.keys())[ax_i]\n",
    "    #Prepare DF\n",
    "    tempS1 = tempD1[sex][yvar_model]\n",
    "    tempS1.name = 'Sex-specific'\n",
    "    tempS2 = tempD1['Both sex'][yvar_model]\n",
    "    tempS2.name = 'Sex-mixed'\n",
    "    tempDF = pd.merge(tempS1, tempS2, left_index=True, right_index=True, how='inner')\n",
    "    #Y=X as reference\n",
    "    ax.plot([range_min, range_max], [range_min, range_max], color='black', linestyle=(0, (1, 2)))\n",
    "    #Regplot\n",
    "    sns.regplot(data=tempDF, x='Sex-specific', y='Sex-mixed', color=tempD2[sex],\n",
    "                scatter=True, fit_reg=True, ci=95, truncate=False, marker='o',\n",
    "                scatter_kws={'alpha':0.2, 'edgecolor':'k', 's':30}, ax=ax)\n",
    "    if ax_i%len(tempD2)==0:\n",
    "        plt.setp(ax, xlabel='', ylabel='Sex-mixed b'+axis_label)\n",
    "    else:\n",
    "        plt.setp(ax, xlabel='', ylabel='')\n",
    "        plt.setp(ax.get_yticklabels(), visible=False)\n",
    "    ##Annotate Pearson's correlation\n",
    "    pearson_r, pval = stats.pearsonr(tempDF['Sex-specific'], tempDF['Sex-mixed'])\n",
    "    r_text = str(Decimal(str(pearson_r)).quantize(Decimal('0.001'), rounding=ROUND_HALF_UP))\n",
    "    if pval==1.0:\n",
    "        pval_text = '1.0'\n",
    "    elif pval==0.0:\n",
    "        pval_text = '0.0'\n",
    "    else:\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",
    "    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='x-small', color='k')\n",
    "    ##Facet label\n",
    "    ax.set_title(sex, {'fontsize':'medium'})\n",
    "    #Save position to generate facet and legend later\n",
    "    if ax_i ==0:\n",
    "        ax0_pos = ax.get_position().bounds\n",
    "    elif ax_i==1:\n",
    "        ax1_pos = ax.get_position().bounds\n",
    "sns.despine()\n",
    "fig.text(x=(ax0_pos[0]+(ax1_pos[0]+ax1_pos[2]))/2, y=ax0_pos[1]-ax0_pos[3]*0.2,#Minor manual adjustment\n",
    "         s='Sex-specific b'+axis_label, fontsize='medium',\n",
    "         verticalalignment='top', horizontalalignment='center', rotation='horizontal')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4-4. Clean beta-coefficient dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempD1 = {'Female':chemBMI_F_bcoefs, 'Male':chemBMI_M_bcoefs, 'BothSex':chemBMI_B_bcoefs}\n",
    "tempD2 = {'Female':chemBMI_F_intercept, 'Male':chemBMI_M_intercept, 'BothSex':chemBMI_B_intercept}\n",
    "yvar_model = 'ChemBMI'\n",
    "model_method = 'Ridge'\n",
    "tempD3 = {}\n",
    "for sex in tempD1.keys():\n",
    "    #Combine variables and intercept\n",
    "    tempDF1 = tempD1[sex]\n",
    "    tempDF2 = tempD2[sex]\n",
    "    tempDF = pd.concat([tempDF1, tempDF2], axis=0)\n",
    "    #Summarize\n",
    "    tempL1 = []\n",
    "    tempL2 = []\n",
    "    tempL3 = []\n",
    "    for row_n in tempDF.index.tolist():\n",
    "        tempL1.append(tempDF.loc[row_n].mean())\n",
    "        tempL2.append(tempDF.loc[row_n].std(ddof=1))#Sample standard deviation\n",
    "        tempL3.append((tempDF.loc[row_n]==0.0).astype('int64').sum())\n",
    "    tempDF['Mean'] = tempL1\n",
    "    tempDF['SD'] = tempL2\n",
    "    tempDF['nZeros'] = tempL3\n",
    "    #Save the cleaned DF\n",
    "    fileDir = './ExportData/'\n",
    "    ipynbName = '220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-Ridge_'\n",
    "    fileName = yvar_model+'-'+sex+'-'+model_method+'bcoefs.tsv'\n",
    "    tempDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "    \n",
    "    #Check\n",
    "    tempDF2 = tempDF.loc[tempDF.index.isin(['Intercept'])]#Retrieve as pd.DataFrame\n",
    "    tempDF = tempDF.drop(index=['Intercept'])\n",
    "    print(sex+':')\n",
    "    print(' - Variables:', len(tempDF))\n",
    "    tempDF1 = tempDF.loc[tempDF['nZeros']!=10]\n",
    "    print(' - Variables with non-zero beta-coefficient:', len(tempDF1),\n",
    "          '(', len(tempDF1)/len(tempDF)*100, '%)')\n",
    "    tempDF1 = tempDF.loc[tempDF['nZeros']==0]\n",
    "    print(' - Variables with non-zero beta-coefficient in all 10 models:', len(tempDF1),\n",
    "          '(', len(tempDF1)/len(tempDF)*100, '%)')\n",
    "    tempDF1 = tempDF1.sort_values(by='Mean', ascending=False)\n",
    "    display(tempDF1)\n",
    "    display(tempDF2)#Intercept\n",
    "    print('')\n",
    "    \n",
    "    tempD3[sex] = tempDF\n",
    "#Update for using R2 transition analysis\n",
    "chemBMI_F_bcoefs = tempD3['Female']\n",
    "chemBMI_M_bcoefs = tempD3['Male']\n",
    "chemBMI_B_bcoefs = tempD3['BothSex']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Metabolomics, Proteomics, and Clinical labs-combined omics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 5-1. Standardization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Note that all preprocessing, including outlier elimination, missingness filtering, imputation, and standardization, should be performed within the cross-validation fold, not across whole dataset, to minimize the potential data leakage. In this study, however, an external dataset is used for “validation\" of the findings (i.e., a testing set for the fitted models). Hence, the robustness of preprocessing is prioritized.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF = combiDF\n",
    "tempD1 = {'Female':bmiDF_F, 'Male':bmiDF_M, 'BothSex':bmiDF_B}\n",
    "tempD2 = {}\n",
    "for sex in tempD1.keys():\n",
    "    tempDF1 = tempD1[sex]\n",
    "    tempDF1 = tempDF.loc[tempDF1.index.tolist()]\n",
    "    #Z-score transformation\n",
    "    scaler = StandardScaler(copy=True, with_mean=True, with_std=True)\n",
    "    tempA = scaler.fit_transform(tempDF1)#Column direction\n",
    "    tempDF1 = pd.DataFrame(data=tempA, index=tempDF1.index, columns=tempDF1.columns)\n",
    "    tempD2[sex] = tempDF1\n",
    "    \n",
    "    #Confirmation\n",
    "    print(sex+':', tempDF1.shape)\n",
    "    display(tempDF1.describe())\n",
    "    sns.set(style='ticks', font='Arial', context='notebook')\n",
    "    plt.figure(figsize=(4, 3))\n",
    "    for col_i in range(3):\n",
    "        sns.distplot(tempDF1.iloc[:, col_i], label=tempDF1.columns[col_i])\n",
    "    sns.despine()\n",
    "    plt.xlabel(r'$Z$'+'-score')\n",
    "    plt.ylabel('Density')\n",
    "    plt.legend(bbox_to_anchor=(1, 0.5), loc='center left', borderaxespad=1)\n",
    "    plt.show()\n",
    "    print('')\n",
    "\n",
    "combiDF_F = tempD2['Female']\n",
    "combiDF_M = tempD2['Male']\n",
    "combiDF_B = tempD2['BothSex']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5-2. Ridge with cross-validation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Note that, in contrast to LassoCV and similarly to LinearRegression, the current RidgeCV seems to produce 2D ndarray although target y is 1D.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Female model\n",
    "tempDF1 = combiDF_F#Standardized independent variables\n",
    "tempDF2 = bmiDF_F#Not-standardized dependent variable and info about testing set\n",
    "yvar = 'BaseBMI'\n",
    "nmodels = len(tempDF2['Testing'].unique())#Reset just in case\n",
    "ncvs = 10\n",
    "model = RidgeCV(alphas=np.geomspace(1, 1000, num=200), fit_intercept=True,\n",
    "                normalize=False, gcv_mode='auto', cv=ncvs)\n",
    "yvar_model = 'BaseCombiBMI'\n",
    "\n",
    "#Perform Ridge\n",
    "bcoefDF = pd.DataFrame(index=tempDF1.columns).astype('float64')#For beta-coefficients\n",
    "bcoefDF.index.rename('Variable', inplace=True)\n",
    "interceptDF = pd.DataFrame(index=['Intercept']).astype('float64')#For intercept\n",
    "interceptDF.index.rename('Variable', inplace=True)\n",
    "scoreL = []#For the coefficient of determination R2\n",
    "predictS = pd.Series(name='log_'+yvar_model).astype('float64')#For predictions\n",
    "t_start = time.time()\n",
    "for model_k in range(nmodels):\n",
    "    #Prepare training and testing (hold-out) datasets in model k\n",
    "    model_n = 'Model_'+str(model_k+1).zfill(2)\n",
    "    yDF_train = tempDF2.loc[tempDF2['Testing']!=model_n]\n",
    "    yDF_test = tempDF2.loc[tempDF2['Testing']==model_n]\n",
    "    xDF_train = tempDF1.loc[yDF_train.index.tolist()]\n",
    "    xDF_test = tempDF1.loc[yDF_test.index.tolist()]\n",
    "    #Retrieve the dependent variable\n",
    "    yDF_train = pd.DataFrame(yDF_train['log_'+yvar])#Not Series but DF\n",
    "    yDF_test = pd.DataFrame(yDF_test['log_'+yvar])#Not Series but DF\n",
    "    #Fitting model with cross-validation using training dataset (i.e., internal training/validation datasets)\n",
    "    model.fit(xDF_train, yDF_train)\n",
    "    #Check the penalization amount decided by cross validation\n",
    "    print(model_n+': alpha =', model.alpha_)\n",
    "    #Save parameters\n",
    "    bcoefDF[model_n] = model.coef_.flatten()#weighted vector\n",
    "    interceptDF[model_n] = model.intercept_#Independent term in decision function\n",
    "    #Evaluation with testing (hold-out) dataset\n",
    "    scoreL.append(model.score(xDF_test, yDF_test))\n",
    "    #Prediction for testing dataset using the fitted model k\n",
    "    tempS = pd.Series(model.predict(xDF_test).flatten(),\n",
    "                      index=xDF_test.index, name='log_'+yvar_model)\n",
    "    predictS = pd.concat([predictS, tempS], axis=0)\n",
    "t_elapsed = time.time() - t_start\n",
    "print('Elapsed time for '+str(nmodels)+' models of '+str(ncvs)+'-fold CV Ridge:',\n",
    "      round(t_elapsed//60), 'min', round(t_elapsed%60, 1), 'sec')\n",
    "\n",
    "#Clean prediction DF\n",
    "tempDF = pd.merge(tempDF2[['Testing', 'log_'+yvar]], predictS,\n",
    "                  left_index=True, right_index=True, how='left')\n",
    "##Convert to original scale\n",
    "tempDF[yvar] = np.e**tempDF['log_'+yvar]\n",
    "tempDF[yvar_model] = np.e**tempDF['log_'+yvar_model]\n",
    "display(tempDF)\n",
    "\n",
    "#Save the cleaned prediction DF\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-Ridge_'\n",
    "fileName = yvar_model.replace('Base', '')+'-Female.tsv'\n",
    "tempDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "combiBMI_F_bcoefs = bcoefDF\n",
    "combiBMI_F_intercept = interceptDF\n",
    "combiBMI_F_R2 = scoreL\n",
    "combiBMI_F = tempDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Male model\n",
    "tempDF1 = combiDF_M#Standardized independent variables\n",
    "tempDF2 = bmiDF_M#Not-standardized dependent variable and info about testing set\n",
    "yvar = 'BaseBMI'\n",
    "nmodels = len(tempDF2['Testing'].unique())#Reset just in case\n",
    "ncvs = 10\n",
    "model = RidgeCV(alphas=np.geomspace(1, 1000, num=200), fit_intercept=True,\n",
    "                normalize=False, gcv_mode='auto', cv=ncvs)\n",
    "yvar_model = 'BaseCombiBMI'\n",
    "\n",
    "#Perform Ridge\n",
    "bcoefDF = pd.DataFrame(index=tempDF1.columns).astype('float64')#For beta-coefficients\n",
    "bcoefDF.index.rename('Variable', inplace=True)\n",
    "interceptDF = pd.DataFrame(index=['Intercept']).astype('float64')#For intercept\n",
    "interceptDF.index.rename('Variable', inplace=True)\n",
    "scoreL = []#For the coefficient of determination R2\n",
    "predictS = pd.Series(name='log_'+yvar_model).astype('float64')#For predictions\n",
    "t_start = time.time()\n",
    "for model_k in range(nmodels):\n",
    "    #Prepare training and testing (hold-out) datasets in model k\n",
    "    model_n = 'Model_'+str(model_k+1).zfill(2)\n",
    "    yDF_train = tempDF2.loc[tempDF2['Testing']!=model_n]\n",
    "    yDF_test = tempDF2.loc[tempDF2['Testing']==model_n]\n",
    "    xDF_train = tempDF1.loc[yDF_train.index.tolist()]\n",
    "    xDF_test = tempDF1.loc[yDF_test.index.tolist()]\n",
    "    #Retrieve the dependent variable\n",
    "    yDF_train = pd.DataFrame(yDF_train['log_'+yvar])#Not Series but DF\n",
    "    yDF_test = pd.DataFrame(yDF_test['log_'+yvar])#Not Series but DF\n",
    "    #Fitting model with cross-validation using training dataset (i.e., internal training/validation datasets)\n",
    "    model.fit(xDF_train, yDF_train)\n",
    "    #Check the penalization amount decided by cross validation\n",
    "    print(model_n+': alpha =', model.alpha_)\n",
    "    #Save parameters\n",
    "    bcoefDF[model_n] = model.coef_.flatten()#weighted vector\n",
    "    interceptDF[model_n] = model.intercept_#Independent term in decision function\n",
    "    #Evaluation with testing (hold-out) dataset\n",
    "    scoreL.append(model.score(xDF_test, yDF_test))\n",
    "    #Prediction for testing dataset using the fitted model k\n",
    "    tempS = pd.Series(model.predict(xDF_test).flatten(),\n",
    "                      index=xDF_test.index, name='log_'+yvar_model)\n",
    "    predictS = pd.concat([predictS, tempS], axis=0)\n",
    "t_elapsed = time.time() - t_start\n",
    "print('Elapsed time for '+str(nmodels)+' models of '+str(ncvs)+'-fold CV Ridge:',\n",
    "      round(t_elapsed//60), 'min', round(t_elapsed%60, 1), 'sec')\n",
    "\n",
    "#Clean prediction DF\n",
    "tempDF = pd.merge(tempDF2[['Testing', 'log_'+yvar]], predictS,\n",
    "                  left_index=True, right_index=True, how='left')\n",
    "##Convert to original scale\n",
    "tempDF[yvar] = np.e**tempDF['log_'+yvar]\n",
    "tempDF[yvar_model] = np.e**tempDF['log_'+yvar_model]\n",
    "display(tempDF)\n",
    "\n",
    "#Save the cleaned prediction DF\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-Ridge_'\n",
    "fileName = yvar_model.replace('Base', '')+'-Male.tsv'\n",
    "tempDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "combiBMI_M_bcoefs = bcoefDF\n",
    "combiBMI_M_intercept = interceptDF\n",
    "combiBMI_M_R2 = scoreL\n",
    "combiBMI_M = tempDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Both sex model\n",
    "tempDF1 = combiDF_B#Standardized independent variables\n",
    "tempDF2 = bmiDF_B#Not-standardized dependent variable and info about testing set\n",
    "yvar = 'BaseBMI'\n",
    "nmodels = len(tempDF2['Testing'].unique())#Reset just in case\n",
    "ncvs = 10\n",
    "model = RidgeCV(alphas=np.geomspace(1, 1000, num=200), fit_intercept=True,\n",
    "                normalize=False, gcv_mode='auto', cv=ncvs)\n",
    "yvar_model = 'BaseCombiBMI'\n",
    "\n",
    "#Perform Ridge\n",
    "bcoefDF = pd.DataFrame(index=tempDF1.columns).astype('float64')#For beta-coefficients\n",
    "bcoefDF.index.rename('Variable', inplace=True)\n",
    "interceptDF = pd.DataFrame(index=['Intercept']).astype('float64')#For intercept\n",
    "interceptDF.index.rename('Variable', inplace=True)\n",
    "scoreL = []#For the coefficient of determination R2\n",
    "predictS = pd.Series(name='log_'+yvar_model).astype('float64')#For predictions\n",
    "t_start = time.time()\n",
    "for model_k in range(nmodels):\n",
    "    #Prepare training and testing (hold-out) datasets in model k\n",
    "    model_n = 'Model_'+str(model_k+1).zfill(2)\n",
    "    yDF_train = tempDF2.loc[tempDF2['Testing']!=model_n]\n",
    "    yDF_test = tempDF2.loc[tempDF2['Testing']==model_n]\n",
    "    xDF_train = tempDF1.loc[yDF_train.index.tolist()]\n",
    "    xDF_test = tempDF1.loc[yDF_test.index.tolist()]\n",
    "    #Retrieve the dependent variable\n",
    "    yDF_train = pd.DataFrame(yDF_train['log_'+yvar])#Not Series but DF\n",
    "    yDF_test = pd.DataFrame(yDF_test['log_'+yvar])#Not Series but DF\n",
    "    #Fitting model with cross-validation using training dataset (i.e., internal training/validation datasets)\n",
    "    model.fit(xDF_train, yDF_train)\n",
    "    #Check the penalization amount decided by cross validation\n",
    "    print(model_n+': alpha =', model.alpha_)\n",
    "    #Save parameters\n",
    "    bcoefDF[model_n] = model.coef_.flatten()#weighted vector\n",
    "    interceptDF[model_n] = model.intercept_#Independent term in decision function\n",
    "    #Evaluation with testing (hold-out) dataset\n",
    "    scoreL.append(model.score(xDF_test, yDF_test))\n",
    "    #Prediction for testing dataset using the fitted model k\n",
    "    tempS = pd.Series(model.predict(xDF_test).flatten(),\n",
    "                      index=xDF_test.index, name='log_'+yvar_model)\n",
    "    predictS = pd.concat([predictS, tempS], axis=0)\n",
    "t_elapsed = time.time() - t_start\n",
    "print('Elapsed time for '+str(nmodels)+' models of '+str(ncvs)+'-fold CV Ridge:',\n",
    "      round(t_elapsed//60), 'min', round(t_elapsed%60, 1), 'sec')\n",
    "\n",
    "#Clean prediction DF\n",
    "tempDF = pd.merge(tempDF2[['Testing', 'log_'+yvar]], predictS,\n",
    "                  left_index=True, right_index=True, how='left')\n",
    "##Convert to original scale\n",
    "tempDF[yvar] = np.e**tempDF['log_'+yvar]\n",
    "tempDF[yvar_model] = np.e**tempDF['log_'+yvar_model]\n",
    "display(tempDF)\n",
    "\n",
    "#Save the cleaned prediction DF\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-Ridge_'\n",
    "fileName = yvar_model.replace('Base', '')+'-BothSex.tsv'\n",
    "tempDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "combiBMI_B_bcoefs = bcoefDF\n",
    "combiBMI_B_intercept = interceptDF\n",
    "combiBMI_B_R2 = scoreL\n",
    "combiBMI_B = tempDF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5-3. Prediction accuracy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Summary\n",
    "tempD1 = {'Female':combiBMI_F_R2, 'Male':combiBMI_M_R2, 'Both sex':combiBMI_B_R2}\n",
    "tempD2 = {'Female':combiBMI_F, 'Male':combiBMI_M, 'Both sex':combiBMI_B}\n",
    "yvar = 'BaseBMI'\n",
    "yvar_model = 'BaseCombiBMI'\n",
    "\n",
    "for sex in tempD1.keys():\n",
    "    tempL = tempD1[sex]\n",
    "    print(sex+' 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",
    "    tempDF = tempD2[sex]\n",
    "    print(' - Observed vs. predicted log_'+yvar+': (Pearson\\'s r, P) =',\n",
    "          stats.pearsonr(tempDF['log_'+yvar], tempDF['log_'+yvar_model]))\n",
    "    print(' - Observed vs. predicted '+yvar+': (Pearson\\'s r, P) =',\n",
    "          stats.pearsonr(tempDF[yvar], tempDF[yvar_model]))\n",
    "    display(tempDF.describe())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Check difference between sex-specific and sex-mixed models for now. Note that this is a rough comparison because the sample size is different.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Prepare DF\n",
    "tempDF = pd.DataFrame({'Female':combiBMI_F_R2, 'Male':combiBMI_M_R2, 'Both sex':combiBMI_B_R2})\n",
    "tempDF = tempDF.melt(var_name='Cohort', value_name='R2', value_vars=tempDF.columns.tolist())\n",
    "\n",
    "#Plot\n",
    "sns.set(style='ticks', font='Arial', context='notebook')\n",
    "plt.figure(figsize=(3, 1))\n",
    "sns.barplot(data=tempDF, y='Cohort', x='R2', hue='Cohort', palette='Set1', dodge=False, edgecolor='black',\n",
    "            ci=95, capsize=0.4, errwidth=1.5, errcolor='black')\n",
    "sns.stripplot(data=tempDF, y='Cohort', x='R2', hue='Cohort', dodge=False, size=8, edgecolor='black',\n",
    "              linewidth=1, alpha=0.4, palette={'Female':'gray', 'Male':'gray', 'Both sex':'gray'})\n",
    "sns.despine()\n",
    "plt.xlabel('Out-of-sample '+r'$R^2$'+' in 10 Ridge models\\n[mean with 95% CI]')\n",
    "plt.ylabel('')\n",
    "plt.legend('', frameon=False)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempD1 = {'Female':combiBMI_F, 'Male':combiBMI_M, 'Both sex':combiBMI_B}\n",
    "yvar = 'BaseBMI'\n",
    "yvar_model = 'BaseCombiBMI'\n",
    "range_min = np.min([df[var].min() for df in tempD1.values() for var in [yvar, yvar_model]])\n",
    "range_max = np.max([df[var].max() for df in tempD1.values() for var in [yvar, yvar_model]])\n",
    "axis_label = 'BMI [kg m'+r'$^{-2}$'+']'\n",
    "\n",
    "#Plot measured vs. predicted per model\n",
    "tempD2 = {'Female':'tab:red', 'Male':'tab:blue', 'Both sex':'tab:green'}\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "fig, axes = plt.subplots(nrows=1, ncols=len(tempD2), figsize=(10, 3), sharex=True, sharey=True,\n",
    "                         gridspec_kw={'hspace':0.1, 'wspace':0.1})\n",
    "for ax_i, ax in enumerate(axes.flat):\n",
    "    sex = list(tempD2.keys())[ax_i]\n",
    "    tempDF = tempD1[sex]\n",
    "    #Y=X as reference\n",
    "    ax.plot([range_min, range_max], [range_min, range_max], color='black', linestyle=(0, (1, 2)))\n",
    "    #Regplot\n",
    "    sns.regplot(data=tempDF, x=yvar_model, y=yvar, color=tempD2[sex],\n",
    "                scatter=True, fit_reg=True, ci=95, truncate=False, marker='o',\n",
    "                scatter_kws={'alpha':0.2, 'edgecolor':'k', 's':30}, ax=ax)\n",
    "    if ax_i%len(tempD2)==0:\n",
    "        plt.setp(ax, xlabel='', ylabel='Measured '+axis_label)\n",
    "    else:\n",
    "        plt.setp(ax, xlabel='', ylabel='')\n",
    "        plt.setp(ax.get_yticklabels(), visible=False)\n",
    "    ##Annotate Pearson's correlation\n",
    "    pearson_r, pval = stats.pearsonr(tempDF[yvar_model], tempDF[yvar])\n",
    "    r_text = str(Decimal(str(pearson_r)).quantize(Decimal('0.001'), rounding=ROUND_HALF_UP))\n",
    "    if pval==1.0:\n",
    "        pval_text = '1.0'\n",
    "    elif pval==0.0:\n",
    "        pval_text = '0.0'\n",
    "    else:\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",
    "    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='x-small', color='k')\n",
    "    ##Facet label\n",
    "    ax.set_title(sex, {'fontsize':'medium'})\n",
    "    #Save position to generate facet and legend later\n",
    "    if ax_i ==0:\n",
    "        ax0_pos = ax.get_position().bounds\n",
    "    elif ax_i==2:\n",
    "        ax2_pos = ax.get_position().bounds\n",
    "sns.despine()\n",
    "fig.text(x=(ax0_pos[0]+(ax2_pos[0]+ax2_pos[2]))/2, y=ax0_pos[1]-ax0_pos[3]*0.2,#Minor manual adjustment\n",
    "         s='Predicted '+axis_label, fontsize='medium',\n",
    "         verticalalignment='top', horizontalalignment='center', rotation='horizontal')\n",
    "plt.show()\n",
    "\n",
    "#Plot difference b/w sex-specific and sex-mixed models\n",
    "tempD2 = {'Female':'tab:red', 'Male':'tab:blue'}\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "fig, axes = plt.subplots(nrows=1, ncols=len(tempD2), figsize=(6.5, 3), sharex=True, sharey=True,\n",
    "                         gridspec_kw={'hspace':0.1, 'wspace':0.1})\n",
    "for ax_i, ax in enumerate(axes.flat):\n",
    "    sex = list(tempD2.keys())[ax_i]\n",
    "    #Prepare DF\n",
    "    tempS1 = tempD1[sex][yvar_model]\n",
    "    tempS1.name = 'Sex-specific'\n",
    "    tempS2 = tempD1['Both sex'][yvar_model]\n",
    "    tempS2.name = 'Sex-mixed'\n",
    "    tempDF = pd.merge(tempS1, tempS2, left_index=True, right_index=True, how='inner')\n",
    "    #Y=X as reference\n",
    "    ax.plot([range_min, range_max], [range_min, range_max], color='black', linestyle=(0, (1, 2)))\n",
    "    #Regplot\n",
    "    sns.regplot(data=tempDF, x='Sex-specific', y='Sex-mixed', color=tempD2[sex],\n",
    "                scatter=True, fit_reg=True, ci=95, truncate=False, marker='o',\n",
    "                scatter_kws={'alpha':0.2, 'edgecolor':'k', 's':30}, ax=ax)\n",
    "    if ax_i%len(tempD2)==0:\n",
    "        plt.setp(ax, xlabel='', ylabel='Sex-mixed b'+axis_label)\n",
    "    else:\n",
    "        plt.setp(ax, xlabel='', ylabel='')\n",
    "        plt.setp(ax.get_yticklabels(), visible=False)\n",
    "    ##Annotate Pearson's correlation\n",
    "    pearson_r, pval = stats.pearsonr(tempDF['Sex-specific'], tempDF['Sex-mixed'])\n",
    "    r_text = str(Decimal(str(pearson_r)).quantize(Decimal('0.001'), rounding=ROUND_HALF_UP))\n",
    "    if pval==1.0:\n",
    "        pval_text = '1.0'\n",
    "    elif pval==0.0:\n",
    "        pval_text = '0.0'\n",
    "    else:\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",
    "    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='x-small', color='k')\n",
    "    ##Facet label\n",
    "    ax.set_title(sex, {'fontsize':'medium'})\n",
    "    #Save position to generate facet and legend later\n",
    "    if ax_i ==0:\n",
    "        ax0_pos = ax.get_position().bounds\n",
    "    elif ax_i==1:\n",
    "        ax1_pos = ax.get_position().bounds\n",
    "sns.despine()\n",
    "fig.text(x=(ax0_pos[0]+(ax1_pos[0]+ax1_pos[2]))/2, y=ax0_pos[1]-ax0_pos[3]*0.2,#Minor manual adjustment\n",
    "         s='Sex-specific b'+axis_label, fontsize='medium',\n",
    "         verticalalignment='top', horizontalalignment='center', rotation='horizontal')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5-4. Clean beta-coefficient dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempD1 = {'Female':combiBMI_F_bcoefs, 'Male':combiBMI_M_bcoefs, 'BothSex':combiBMI_B_bcoefs}\n",
    "tempD2 = {'Female':combiBMI_F_intercept, 'Male':combiBMI_M_intercept, 'BothSex':combiBMI_B_intercept}\n",
    "yvar_model = 'CombiBMI'\n",
    "model_method = 'Ridge'\n",
    "tempD3 = {}\n",
    "for sex in tempD1.keys():\n",
    "    #Combine variables and intercept\n",
    "    tempDF1 = tempD1[sex]\n",
    "    tempDF2 = tempD2[sex]\n",
    "    tempDF = pd.concat([tempDF1, tempDF2], axis=0)\n",
    "    #Summarize\n",
    "    tempL1 = []\n",
    "    tempL2 = []\n",
    "    tempL3 = []\n",
    "    for row_n in tempDF.index.tolist():\n",
    "        tempL1.append(tempDF.loc[row_n].mean())\n",
    "        tempL2.append(tempDF.loc[row_n].std(ddof=1))#Sample standard deviation\n",
    "        tempL3.append((tempDF.loc[row_n]==0.0).astype('int64').sum())\n",
    "    tempDF['Mean'] = tempL1\n",
    "    tempDF['SD'] = tempL2\n",
    "    tempDF['nZeros'] = tempL3\n",
    "    #Save the cleaned DF\n",
    "    fileDir = './ExportData/'\n",
    "    ipynbName = '220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-Ridge_'\n",
    "    fileName = yvar_model+'-'+sex+'-'+model_method+'bcoefs.tsv'\n",
    "    tempDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "    \n",
    "    #Check\n",
    "    tempDF2 = tempDF.loc[tempDF.index.isin(['Intercept'])]#Retrieve as pd.DataFrame\n",
    "    tempDF = tempDF.drop(index=['Intercept'])\n",
    "    print(sex+':')\n",
    "    print(' - Variables:', len(tempDF))\n",
    "    tempDF1 = tempDF.loc[tempDF['nZeros']!=10]\n",
    "    print(' - Variables with non-zero beta-coefficient:', len(tempDF1),\n",
    "          '(', len(tempDF1)/len(tempDF)*100, '%)')\n",
    "    tempDF1 = tempDF.loc[tempDF['nZeros']==0]\n",
    "    print(' - Variables with non-zero beta-coefficient in all 10 models:', len(tempDF1),\n",
    "          '(', len(tempDF1)/len(tempDF)*100, '%)')\n",
    "    tempDF1 = tempDF1.sort_values(by='Mean', ascending=False)\n",
    "    display(tempDF1)\n",
    "    display(tempDF2)#Intercept\n",
    "    print('')\n",
    "    \n",
    "    tempD3[sex] = tempDF\n",
    "#Update for using R2 transition analysis\n",
    "combiBMI_F_bcoefs = tempD3['Female']\n",
    "combiBMI_M_bcoefs = tempD3['Male']\n",
    "combiBMI_B_bcoefs = tempD3['BothSex']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Standard clinical measures"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> A model predicting BMI is generated using ordinary least squares (OLS) linear regression with sex, age, HDL-cholesterol, LDL-cholesterol, triglycerides (TG), glucose, insulin, and HOMA-IR. To directely compare with the Ridge models, regression analysis is applied to 10 split set.  \n",
    ">\n",
    "> Of note, a previous study (Cirulli, E.T. et al. Cell Metab. 2019) generated a “conventional model\" from sex, age, HDL, LDL, TG, and total cholesterol, which explained 31% variance in BMI (probably, in-sample R2). However, their model deesn't include the information about glucose, which is broadly used for metabolic health risks (Stefan, N. et al. Cell Metab. 2017). Also, total cholesterol was collinear with LDL and TG, at least in the Arivale cohort, enough to make negative beta-coefficients for them. HOMA-IR works a kind of interaction term between glucose and insulin.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6-1. Prepare independent variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF1 = bmiDF\n",
    "tempDF2 = chemDF\n",
    "\n",
    "#Select the independent variables\n",
    "tempD = {'Sex':'Sex', 'BaseAge':'Age'}\n",
    "tempDF1 = tempDF1[tempD.keys()]\n",
    "tempDF1 = tempDF1.rename(columns=tempD)\n",
    "tempD = {'HDL CHOL DIRECT':'HDL-cholesterol',\n",
    "         'LDL-CHOL CALCULATION':'LDL-cholesterol',\n",
    "         'TRIGLYCERIDES':'Triglycerides',\n",
    "         'GLUCOSE':'Glucose',\n",
    "         'INSULIN':'Insulin',\n",
    "         'HOMA-IR':'HOMA-IR'}\n",
    "tempDF2 = tempDF2[list(tempD.keys())]\n",
    "tempDF2 = tempDF2.rename(columns=tempD)\n",
    "tempDF = pd.merge(tempDF1, tempDF2, left_index=True, right_index=True, how='inner')\n",
    "\n",
    "display(tempDF.describe(include='all'))\n",
    "\n",
    "standDF = tempDF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 6-2. Standardization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Note that all preprocessing, including outlier elimination, missingness filtering, imputation, and standardization, should be performed within the cross-validation fold, not across whole dataset, to minimize the potential data leakage. In this study, however, an external dataset is used for “validation\" of the findings (i.e., a testing set for the fitted models). Hence, the robustness of preprocessing is prioritized.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF = standDF\n",
    "tempD1 = {'Female':bmiDF_F, 'Male':bmiDF_M, 'BothSex':bmiDF_B}\n",
    "tempD2 = {}\n",
    "for sex in tempD1.keys():\n",
    "    tempDF1 = tempD1[sex]\n",
    "    tempDF1 = tempDF.loc[tempDF1.index.tolist()]\n",
    "    #Z-score transformation\n",
    "    tempDF2 = tempDF1.select_dtypes(include=[np.number])\n",
    "    scaler = StandardScaler(copy=True, with_mean=True, with_std=True)\n",
    "    tempA = scaler.fit_transform(tempDF2)#Column direction\n",
    "    tempDF2 = pd.DataFrame(data=tempA, index=tempDF2.index, columns=tempDF2.columns)\n",
    "    ##Recover categorical covariates\n",
    "    tempDF3 = tempDF1.select_dtypes(exclude=[np.number])\n",
    "    tempDF1 = pd.merge(tempDF2, tempDF3, left_index=True, right_index=True, how='left')\n",
    "    tempD2[sex] = tempDF1\n",
    "    \n",
    "    #Confirmation\n",
    "    print(sex+':', tempDF1.shape)\n",
    "    display(tempDF1.describe(include='all'))\n",
    "    tempDF1 = tempDF1.select_dtypes(include=[np.number])\n",
    "    sns.set(style='ticks', font='Arial', context='notebook')\n",
    "    plt.figure(figsize=(4, 3))\n",
    "    for col_i in range(len(tempDF1.columns)):\n",
    "        sns.distplot(tempDF1.iloc[:, col_i], label=tempDF1.columns[col_i])\n",
    "    sns.despine()\n",
    "    plt.xlabel(r'$Z$'+'-score')\n",
    "    plt.ylabel('Density')\n",
    "    plt.legend(bbox_to_anchor=(1, 0.5), loc='center left', borderaxespad=1)\n",
    "    plt.show()\n",
    "    print('')\n",
    "\n",
    "standDF_F = tempD2['Female']\n",
    "standDF_M = tempD2['Male']\n",
    "standDF_B = tempD2['BothSex']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6-3. One-hot encoding for categorical variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> While statsmodels automatically recognizes categorical variables, one-hot encoding is required in scikit-learn.  \n",
    "–> In this case, only sex is the categorical variable; hence, map manually! (cf. in many cases, category_encoders is more useful than sklearn.preprocessing or pandas.get_dummies.)  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempD1 = {'Female':standDF_F, 'Male':standDF_M, 'BothSex':standDF_B}\n",
    "tempD2 = {}\n",
    "for sex in tempD1.keys():\n",
    "    tempDF = tempD1[sex]\n",
    "    #One-hot encoding for sex\n",
    "    tempD = {'F':0, 'M':1}\n",
    "    tempDF['Sex'] = tempDF['Sex'].map(tempD)\n",
    "    tempD2[sex] = tempDF\n",
    "#Update\n",
    "standDF_F = tempD2['Female']\n",
    "standDF_M = tempD2['Male']\n",
    "standDF_B = tempD2['BothSex']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 6-4. OLS linear regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> In scikit-learn, constant variables are automatically assinged zero to; hence, although sex is included in the model, there is no need to modify the code between sex-specific and sex-mixed version.  \n",
    ">\n",
    "> Note that, in contrast to LassoCV and similarly to RidgeCV, the current LinearRegression seems to produce 2D ndarray although target y is 1D.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Female model\n",
    "tempDF1 = standDF_F#Standardized independent variables\n",
    "tempDF2 = bmiDF_F#Not-standardized dependent variable and info about testing set\n",
    "yvar = 'BaseBMI'\n",
    "nmodels = len(tempDF2['Testing'].unique())#Reset just in case\n",
    "model = LinearRegression(fit_intercept=True, normalize=False)\n",
    "yvar_model = 'BaseStandBMI'\n",
    "\n",
    "#Perform OLS linear regression\n",
    "bcoefDF = pd.DataFrame(index=tempDF1.columns).astype('float64')#For beta-coefficients\n",
    "bcoefDF.index.rename('Variable', inplace=True)\n",
    "interceptDF = pd.DataFrame(index=['Intercept']).astype('float64')#For intercept\n",
    "interceptDF.index.rename('Variable', inplace=True)\n",
    "scoreL = []#For the coefficient of determination R2\n",
    "predictS = pd.Series(name='log_'+yvar_model).astype('float64')#For predictions\n",
    "t_start = time.time()\n",
    "for model_k in range(nmodels):\n",
    "    #Prepare training and testing (hold-out) datasets in model k\n",
    "    model_n = 'Model_'+str(model_k+1).zfill(2)\n",
    "    yDF_train = tempDF2.loc[tempDF2['Testing']!=model_n]\n",
    "    yDF_test = tempDF2.loc[tempDF2['Testing']==model_n]\n",
    "    xDF_train = tempDF1.loc[yDF_train.index.tolist()]\n",
    "    xDF_test = tempDF1.loc[yDF_test.index.tolist()]\n",
    "    #Retrieve the dependent variable\n",
    "    yDF_train = pd.DataFrame(yDF_train['log_'+yvar])#Not Series but DF\n",
    "    yDF_test = pd.DataFrame(yDF_test['log_'+yvar])#Not Series but DF\n",
    "    #Fitting model using training dataset\n",
    "    model.fit(xDF_train, yDF_train)\n",
    "    #Save parameters\n",
    "    bcoefDF[model_n] = model.coef_.flatten()#Estimated coefficients for the linear regression problem\n",
    "    interceptDF[model_n] = model.intercept_#Independent term in the linear model\n",
    "    #Evaluation with testing (hold-out) dataset\n",
    "    scoreL.append(model.score(xDF_test, yDF_test))\n",
    "    #Prediction for testing dataset using the fitted model k\n",
    "    tempS = pd.Series(model.predict(xDF_test).flatten(),\n",
    "                      index=xDF_test.index, name='log_'+yvar_model)\n",
    "    predictS = pd.concat([predictS, tempS], axis=0)\n",
    "t_elapsed = time.time() - t_start\n",
    "print('Elapsed time for '+str(nmodels)+' models of OLS linear regression:',\n",
    "      round(t_elapsed//60), 'min', round(t_elapsed%60, 1), 'sec')\n",
    "\n",
    "#Clean prediction DF\n",
    "tempDF = pd.merge(tempDF2[['Testing', 'log_'+yvar]], predictS,\n",
    "                  left_index=True, right_index=True, how='left')\n",
    "##Convert to original scale\n",
    "tempDF[yvar] = np.e**tempDF['log_'+yvar]\n",
    "tempDF[yvar_model] = np.e**tempDF['log_'+yvar_model]\n",
    "display(tempDF)\n",
    "\n",
    "#Save the cleaned prediction DF\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-Ridge_'\n",
    "fileName = yvar_model.replace('Base', '')+'-Female.tsv'\n",
    "tempDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "standBMI_F_bcoefs = bcoefDF\n",
    "standBMI_F_intercept = interceptDF\n",
    "standBMI_F_R2 = scoreL\n",
    "standBMI_F = tempDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Male model\n",
    "tempDF1 = standDF_M#Standardized independent variables\n",
    "tempDF2 = bmiDF_M#Not-standardized dependent variable and info about testing set\n",
    "yvar = 'BaseBMI'\n",
    "nmodels = len(tempDF2['Testing'].unique())#Reset just in case\n",
    "model = LinearRegression(fit_intercept=True, normalize=False)\n",
    "yvar_model = 'BaseStandBMI'\n",
    "\n",
    "#Perform OLS linear regression\n",
    "bcoefDF = pd.DataFrame(index=tempDF1.columns).astype('float64')#For beta-coefficients\n",
    "bcoefDF.index.rename('Variable', inplace=True)\n",
    "interceptDF = pd.DataFrame(index=['Intercept']).astype('float64')#For intercept\n",
    "interceptDF.index.rename('Variable', inplace=True)\n",
    "scoreL = []#For the coefficient of determination R2\n",
    "predictS = pd.Series(name='log_'+yvar_model).astype('float64')#For predictions\n",
    "t_start = time.time()\n",
    "for model_k in range(nmodels):\n",
    "    #Prepare training and testing (hold-out) datasets in model k\n",
    "    model_n = 'Model_'+str(model_k+1).zfill(2)\n",
    "    yDF_train = tempDF2.loc[tempDF2['Testing']!=model_n]\n",
    "    yDF_test = tempDF2.loc[tempDF2['Testing']==model_n]\n",
    "    xDF_train = tempDF1.loc[yDF_train.index.tolist()]\n",
    "    xDF_test = tempDF1.loc[yDF_test.index.tolist()]\n",
    "    #Retrieve the dependent variable\n",
    "    yDF_train = pd.DataFrame(yDF_train['log_'+yvar])#Not Series but DF\n",
    "    yDF_test = pd.DataFrame(yDF_test['log_'+yvar])#Not Series but DF\n",
    "    #Fitting model using training dataset\n",
    "    model.fit(xDF_train, yDF_train)\n",
    "    #Save parameters\n",
    "    bcoefDF[model_n] = model.coef_.flatten()#Estimated coefficients for the linear regression problem\n",
    "    interceptDF[model_n] = model.intercept_#Independent term in the linear model\n",
    "    #Evaluation with testing (hold-out) dataset\n",
    "    scoreL.append(model.score(xDF_test, yDF_test))\n",
    "    #Prediction for testing dataset using the fitted model k\n",
    "    tempS = pd.Series(model.predict(xDF_test).flatten(),\n",
    "                      index=xDF_test.index, name='log_'+yvar_model)\n",
    "    predictS = pd.concat([predictS, tempS], axis=0)\n",
    "t_elapsed = time.time() - t_start\n",
    "print('Elapsed time for '+str(nmodels)+' models of OLS linear regression:',\n",
    "      round(t_elapsed//60), 'min', round(t_elapsed%60, 1), 'sec')\n",
    "\n",
    "#Clean prediction DF\n",
    "tempDF = pd.merge(tempDF2[['Testing', 'log_'+yvar]], predictS,\n",
    "                  left_index=True, right_index=True, how='left')\n",
    "##Convert to original scale\n",
    "tempDF[yvar] = np.e**tempDF['log_'+yvar]\n",
    "tempDF[yvar_model] = np.e**tempDF['log_'+yvar_model]\n",
    "display(tempDF)\n",
    "\n",
    "#Save the cleaned prediction DF\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-Ridge_'\n",
    "fileName = yvar_model.replace('Base', '')+'-Male.tsv'\n",
    "tempDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "standBMI_M_bcoefs = bcoefDF\n",
    "standBMI_M_intercept = interceptDF\n",
    "standBMI_M_R2 = scoreL\n",
    "standBMI_M = tempDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Both sex model\n",
    "tempDF1 = standDF_B#Standardized independent variables\n",
    "tempDF2 = bmiDF_B#Not-standardized dependent variable and info about testing set\n",
    "yvar = 'BaseBMI'\n",
    "nmodels = len(tempDF2['Testing'].unique())#Reset just in case\n",
    "model = LinearRegression(fit_intercept=True, normalize=False)\n",
    "yvar_model = 'BaseStandBMI'\n",
    "\n",
    "#Perform OLS linear regression\n",
    "bcoefDF = pd.DataFrame(index=tempDF1.columns).astype('float64')#For beta-coefficients\n",
    "bcoefDF.index.rename('Variable', inplace=True)\n",
    "interceptDF = pd.DataFrame(index=['Intercept']).astype('float64')#For intercept\n",
    "interceptDF.index.rename('Variable', inplace=True)\n",
    "scoreL = []#For the coefficient of determination R2\n",
    "predictS = pd.Series(name='log_'+yvar_model).astype('float64')#For predictions\n",
    "t_start = time.time()\n",
    "for model_k in range(nmodels):\n",
    "    #Prepare training and testing (hold-out) datasets in model k\n",
    "    model_n = 'Model_'+str(model_k+1).zfill(2)\n",
    "    yDF_train = tempDF2.loc[tempDF2['Testing']!=model_n]\n",
    "    yDF_test = tempDF2.loc[tempDF2['Testing']==model_n]\n",
    "    xDF_train = tempDF1.loc[yDF_train.index.tolist()]\n",
    "    xDF_test = tempDF1.loc[yDF_test.index.tolist()]\n",
    "    #Retrieve the dependent variable\n",
    "    yDF_train = pd.DataFrame(yDF_train['log_'+yvar])#Not Series but DF\n",
    "    yDF_test = pd.DataFrame(yDF_test['log_'+yvar])#Not Series but DF\n",
    "    #Fitting model using training dataset\n",
    "    model.fit(xDF_train, yDF_train)\n",
    "    #Save parameters\n",
    "    bcoefDF[model_n] = model.coef_.flatten()#Estimated coefficients for the linear regression problem\n",
    "    interceptDF[model_n] = model.intercept_#Independent term in the linear model\n",
    "    #Evaluation with testing (hold-out) dataset\n",
    "    scoreL.append(model.score(xDF_test, yDF_test))\n",
    "    #Prediction for testing dataset using the fitted model k\n",
    "    tempS = pd.Series(model.predict(xDF_test).flatten(),\n",
    "                      index=xDF_test.index, name='log_'+yvar_model)\n",
    "    predictS = pd.concat([predictS, tempS], axis=0)\n",
    "t_elapsed = time.time() - t_start\n",
    "print('Elapsed time for '+str(nmodels)+' models of OLS linear regression:',\n",
    "      round(t_elapsed//60), 'min', round(t_elapsed%60, 1), 'sec')\n",
    "\n",
    "#Clean prediction DF\n",
    "tempDF = pd.merge(tempDF2[['Testing', 'log_'+yvar]], predictS,\n",
    "                  left_index=True, right_index=True, how='left')\n",
    "##Convert to original scale\n",
    "tempDF[yvar] = np.e**tempDF['log_'+yvar]\n",
    "tempDF[yvar_model] = np.e**tempDF['log_'+yvar_model]\n",
    "display(tempDF)\n",
    "\n",
    "#Save the cleaned prediction DF\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-Ridge_'\n",
    "fileName = yvar_model.replace('Base', '')+'-BothSex.tsv'\n",
    "tempDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "standBMI_B_bcoefs = bcoefDF\n",
    "standBMI_B_intercept = interceptDF\n",
    "standBMI_B_R2 = scoreL\n",
    "standBMI_B = tempDF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6-5. Prediction accuracy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Summary\n",
    "tempD1 = {'Female':standBMI_F_R2, 'Male':standBMI_M_R2, 'Both sex':standBMI_B_R2}\n",
    "tempD2 = {'Female':standBMI_F, 'Male':standBMI_M, 'Both sex':standBMI_B}\n",
    "yvar = 'BaseBMI'\n",
    "yvar_model = 'BaseStandBMI'\n",
    "\n",
    "for sex in tempD1.keys():\n",
    "    tempL = tempD1[sex]\n",
    "    print(sex+' 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",
    "    tempDF = tempD2[sex]\n",
    "    print(' - Observed vs. predicted log_'+yvar+': (Pearson\\'s r, P) =',\n",
    "          stats.pearsonr(tempDF['log_'+yvar], tempDF['log_'+yvar_model]))\n",
    "    print(' - Observed vs. predicted '+yvar+': (Pearson\\'s r, P) =',\n",
    "          stats.pearsonr(tempDF[yvar], tempDF[yvar_model]))\n",
    "    display(tempDF.describe())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Check difference between sex-specific and sex-mixed models for now. Note that this is a rough comparison because the sample size is different.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Prepare DF\n",
    "tempDF = pd.DataFrame({'Female':standBMI_F_R2, 'Male':standBMI_M_R2, 'Both sex':standBMI_B_R2})\n",
    "tempDF = tempDF.melt(var_name='Cohort', value_name='R2', value_vars=tempDF.columns.tolist())\n",
    "\n",
    "#Plot\n",
    "sns.set(style='ticks', font='Arial', context='notebook')\n",
    "plt.figure(figsize=(3, 1))\n",
    "sns.barplot(data=tempDF, y='Cohort', x='R2', hue='Cohort', palette='Set1', dodge=False, edgecolor='black',\n",
    "            ci=95, capsize=0.4, errwidth=1.5, errcolor='black')\n",
    "sns.stripplot(data=tempDF, y='Cohort', x='R2', hue='Cohort', dodge=False, size=8, edgecolor='black',\n",
    "              linewidth=1, alpha=0.4, palette={'Female':'gray', 'Male':'gray', 'Both sex':'gray'})\n",
    "sns.despine()\n",
    "plt.xlabel('Out-of-sample '+r'$R^2$'+' in 10 OLS-LR models\\n[mean with 95% CI]')\n",
    "plt.ylabel('')\n",
    "plt.legend('', frameon=False)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempD1 = {'Female':standBMI_F, 'Male':standBMI_M, 'Both sex':standBMI_B}\n",
    "yvar = 'BaseBMI'\n",
    "yvar_model = 'BaseStandBMI'\n",
    "range_min = np.min([df[var].min() for df in tempD1.values() for var in [yvar, yvar_model]])\n",
    "range_max = np.max([df[var].max() for df in tempD1.values() for var in [yvar, yvar_model]])\n",
    "axis_label = 'BMI [kg m'+r'$^{-2}$'+']'\n",
    "\n",
    "#Plot measured vs. predicted per model\n",
    "tempD2 = {'Female':'tab:red', 'Male':'tab:blue', 'Both sex':'tab:green'}\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "fig, axes = plt.subplots(nrows=1, ncols=len(tempD2), figsize=(10, 3), sharex=True, sharey=True,\n",
    "                         gridspec_kw={'hspace':0.1, 'wspace':0.1})\n",
    "for ax_i, ax in enumerate(axes.flat):\n",
    "    sex = list(tempD2.keys())[ax_i]\n",
    "    tempDF = tempD1[sex]\n",
    "    #Y=X as reference\n",
    "    ax.plot([range_min, range_max], [range_min, range_max], color='black', linestyle=(0, (1, 2)))\n",
    "    #Regplot\n",
    "    sns.regplot(data=tempDF, x=yvar_model, y=yvar, color=tempD2[sex],\n",
    "                scatter=True, fit_reg=True, ci=95, truncate=False, marker='o',\n",
    "                scatter_kws={'alpha':0.2, 'edgecolor':'k', 's':30}, ax=ax)\n",
    "    if ax_i%len(tempD2)==0:\n",
    "        plt.setp(ax, xlabel='', ylabel='Measured '+axis_label)\n",
    "    else:\n",
    "        plt.setp(ax, xlabel='', ylabel='')\n",
    "        plt.setp(ax.get_yticklabels(), visible=False)\n",
    "    ##Annotate Pearson's correlation\n",
    "    pearson_r, pval = stats.pearsonr(tempDF[yvar_model], tempDF[yvar])\n",
    "    r_text = str(Decimal(str(pearson_r)).quantize(Decimal('0.001'), rounding=ROUND_HALF_UP))\n",
    "    if pval==1.0:\n",
    "        pval_text = '1.0'\n",
    "    elif pval==0.0:\n",
    "        pval_text = '0.0'\n",
    "    else:\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",
    "    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='x-small', color='k')\n",
    "    ##Facet label\n",
    "    ax.set_title(sex, {'fontsize':'medium'})\n",
    "    #Save position to generate facet and legend later\n",
    "    if ax_i ==0:\n",
    "        ax0_pos = ax.get_position().bounds\n",
    "    elif ax_i==2:\n",
    "        ax2_pos = ax.get_position().bounds\n",
    "sns.despine()\n",
    "fig.text(x=(ax0_pos[0]+(ax2_pos[0]+ax2_pos[2]))/2, y=ax0_pos[1]-ax0_pos[3]*0.2,#Minor manual adjustment\n",
    "         s='Predicted '+axis_label, fontsize='medium',\n",
    "         verticalalignment='top', horizontalalignment='center', rotation='horizontal')\n",
    "plt.show()\n",
    "\n",
    "#Plot difference b/w sex-specific and sex-mixed models\n",
    "tempD2 = {'Female':'tab:red', 'Male':'tab:blue'}\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "fig, axes = plt.subplots(nrows=1, ncols=len(tempD2), figsize=(6.5, 3), sharex=True, sharey=True,\n",
    "                         gridspec_kw={'hspace':0.1, 'wspace':0.1})\n",
    "for ax_i, ax in enumerate(axes.flat):\n",
    "    sex = list(tempD2.keys())[ax_i]\n",
    "    #Prepare DF\n",
    "    tempS1 = tempD1[sex][yvar_model]\n",
    "    tempS1.name = 'Sex-specific'\n",
    "    tempS2 = tempD1['Both sex'][yvar_model]\n",
    "    tempS2.name = 'Sex-mixed'\n",
    "    tempDF = pd.merge(tempS1, tempS2, left_index=True, right_index=True, how='inner')\n",
    "    #Y=X as reference\n",
    "    ax.plot([range_min, range_max], [range_min, range_max], color='black', linestyle=(0, (1, 2)))\n",
    "    #Regplot\n",
    "    sns.regplot(data=tempDF, x='Sex-specific', y='Sex-mixed', color=tempD2[sex],\n",
    "                scatter=True, fit_reg=True, ci=95, truncate=False, marker='o',\n",
    "                scatter_kws={'alpha':0.2, 'edgecolor':'k', 's':30}, ax=ax)\n",
    "    if ax_i%len(tempD2)==0:\n",
    "        plt.setp(ax, xlabel='', ylabel='Sex-mixed b'+axis_label)\n",
    "    else:\n",
    "        plt.setp(ax, xlabel='', ylabel='')\n",
    "        plt.setp(ax.get_yticklabels(), visible=False)\n",
    "    ##Annotate Pearson's correlation\n",
    "    pearson_r, pval = stats.pearsonr(tempDF['Sex-specific'], tempDF['Sex-mixed'])\n",
    "    r_text = str(Decimal(str(pearson_r)).quantize(Decimal('0.001'), rounding=ROUND_HALF_UP))\n",
    "    if pval==1.0:\n",
    "        pval_text = '1.0'\n",
    "    elif pval==0.0:\n",
    "        pval_text = '0.0'\n",
    "    else:\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",
    "    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='x-small', color='k')\n",
    "    ##Facet label\n",
    "    ax.set_title(sex, {'fontsize':'medium'})\n",
    "    #Save position to generate facet and legend later\n",
    "    if ax_i ==0:\n",
    "        ax0_pos = ax.get_position().bounds\n",
    "    elif ax_i==1:\n",
    "        ax1_pos = ax.get_position().bounds\n",
    "sns.despine()\n",
    "fig.text(x=(ax0_pos[0]+(ax1_pos[0]+ax1_pos[2]))/2, y=ax0_pos[1]-ax0_pos[3]*0.2,#Minor manual adjustment\n",
    "         s='Sex-specific b'+axis_label, fontsize='medium',\n",
    "         verticalalignment='top', horizontalalignment='center', rotation='horizontal')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6-6. Clean beta-coefficient dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempD1 = {'Female':standBMI_F_bcoefs, 'Male':standBMI_M_bcoefs, 'BothSex':standBMI_B_bcoefs}\n",
    "tempD2 = {'Female':standBMI_F_intercept, 'Male':standBMI_M_intercept, 'BothSex':standBMI_B_intercept}\n",
    "yvar_model = 'StandBMI'\n",
    "model_method = 'OLS'\n",
    "for sex in tempD1.keys():\n",
    "    #Combine variables and intercept\n",
    "    tempDF1 = tempD1[sex]\n",
    "    tempDF2 = tempD2[sex]\n",
    "    tempDF = pd.concat([tempDF1, tempDF2], axis=0)\n",
    "    #Summarize\n",
    "    tempL1 = []\n",
    "    tempL2 = []\n",
    "    tempL3 = []\n",
    "    for row_n in tempDF.index.tolist():\n",
    "        tempL1.append(tempDF.loc[row_n].mean())\n",
    "        tempL2.append(tempDF.loc[row_n].std(ddof=1))#Sample standard deviation\n",
    "        tempL3.append((tempDF.loc[row_n]==0.0).astype('int64').sum())\n",
    "    tempDF['Mean'] = tempL1\n",
    "    tempDF['SD'] = tempL2\n",
    "    tempDF['nZeros'] = tempL3\n",
    "    #Save the cleaned DF\n",
    "    fileDir = './ExportData/'\n",
    "    ipynbName = '220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-Ridge_'\n",
    "    fileName = yvar_model+'-'+sex+'-'+model_method+'bcoefs.tsv'\n",
    "    tempDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "    \n",
    "    #Check\n",
    "    tempDF2 = tempDF.loc[tempDF.index.isin(['Intercept'])]#Retrieve as pd.DataFrame\n",
    "    tempDF = tempDF.drop(index=['Intercept'])\n",
    "    print(sex+':')\n",
    "    print(' - Variables:', len(tempDF))\n",
    "    tempDF1 = tempDF.loc[tempDF['nZeros']!=10]\n",
    "    print(' - Variables with non-zero beta-coefficient:', len(tempDF1),\n",
    "          '(', len(tempDF1)/len(tempDF)*100, '%)')\n",
    "    tempDF1 = tempDF.loc[tempDF['nZeros']==0]\n",
    "    print(' - Variables with non-zero beta-coefficient in all 10 models:', len(tempDF1),\n",
    "          '(', len(tempDF1)/len(tempDF)*100, '%)')\n",
    "    tempDF1 = tempDF1.sort_values(by='Mean', ascending=False)\n",
    "    display(tempDF1)\n",
    "    display(tempDF2)#Intercept\n",
    "    print('')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7. Comparison b/w models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7-1. Out-of-sample R2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Assess significance using two-sided Welch's t-test (using statsmodels library). Note that scipy.stats.ttest_ind() doesn't report degrees of freedom.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempD1 = {'Standard measures':standBMI_B_R2,\n",
    "          'Metabolomics':metBMI_B_R2, 'Proteomics':protBMI_B_R2,\n",
    "          'Clinical labs':chemBMI_B_R2, 'Combined omics':combiBMI_B_R2}\n",
    "tempD2 = {'Standard measures':'0.5',\n",
    "          'Metabolomics':'b', 'Proteomics':'r',\n",
    "          'Clinical labs':'g', 'Combined omics':'m'}\n",
    "\n",
    "#Prepare DF\n",
    "tempDF = pd.DataFrame(tempD1)\n",
    "display(tempDF.describe())\n",
    "\n",
    "#Statistical tests\n",
    "control = list(tempD1.keys())[0]\n",
    "tempL = list(tempD1.keys())[1:]\n",
    "tempDF1 = pd.DataFrame(columns=['Control', 'Contrast', 'Control_N', 'contrast_N', 'DoF', 'tStat', 'Pval'])\n",
    "for contrast in tempL:\n",
    "    tempS1 = tempDF[control]\n",
    "    tempS2 = tempDF[contrast]\n",
    "    #Two-sided Welch's t-test\n",
    "    tstat, pval, dof = weightstats.ttest_ind(tempS1, tempS2,\n",
    "                                             alternative='two-sided', usevar='unequal')\n",
    "    size1 = len(tempS1)\n",
    "    size2 = len(tempS2)\n",
    "    tempDF1.loc[contrast+'-vs-'+control] = [control, contrast, size1, size2, dof, tstat, 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('ComparisonLabel', inplace=True)\n",
    "display(tempDF1)\n",
    "##Save\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-Ridge_'\n",
    "fileName = 'R2-comparison-BothSex.tsv'\n",
    "tempDF1.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "\n",
    "#Plot\n",
    "tempDF = tempDF.melt(var_name='Category', value_name='R2', value_vars=tempDF.columns.tolist())\n",
    "axis_ymin = 0.0\n",
    "axis_ymax = 1.2\n",
    "ymin = 0.0\n",
    "ymax = 0.8\n",
    "yinter = 0.2\n",
    "aline_ymin = 0.8\n",
    "aline_yinter = 0.1\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "plt.figure(figsize=(3.5, 4))\n",
    "sns.barplot(data=tempDF, y='R2', x='Category', order=tempD2.keys(),\n",
    "            hue='Category', hue_order=tempD2.keys(), dodge=False, palette=tempD2,\n",
    "            ci=95, capsize=0.4, errwidth=1.5, errcolor='black', edgecolor='black')\n",
    "p = sns.stripplot(data=tempDF, y='R2', x='Category', order=tempD2.keys(),\n",
    "                  hue='Category', hue_order=tempD2.keys(), dodge=False, jitter=0.3,\n",
    "                  size=5, edgecolor='black', color='gray', linewidth=1, alpha=0.4)\n",
    "sns.despine()\n",
    "p.set(ylim=(axis_ymin, axis_ymax), yticks=np.arange(ymin, ymax+yinter/10, yinter))\n",
    "plt.setp(p.get_xticklabels(), rotation=70,\n",
    "         horizontalalignment='right', verticalalignment='center', rotation_mode='anchor')\n",
    "##P-value annotation\n",
    "for row_i in range(len(tempDF1)):\n",
    "    #Control\n",
    "    group_0 = tempDF1['Control'].iloc[row_i]\n",
    "    xcoord_0 = list(tempD2.keys()).index(group_0)\n",
    "    #Contrast\n",
    "    group_1 = tempDF1['Contrast'].iloc[row_i]\n",
    "    xcoord_1 = list(tempD2.keys()).index(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",
    "    p.plot([xcoord_0, xcoord_0, xcoord_1, xcoord_1],\n",
    "           [ycoord+aline_offset, ycoord+aline_length, ycoord+aline_length, ycoord+aline_offset],\n",
    "           lw=1.5, c='k')\n",
    "    #Retrieve P-value\n",
    "    pval = tempDF1['AdjPval'].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",
    "    p.annotate(label, xy=(xcoord, ycoord+text_offset),\n",
    "               horizontalalignment='center', verticalalignment='bottom',\n",
    "               fontsize=text_size, color='k')\n",
    "plt.ylabel('Out-of-sample '+r'$R^2$')\n",
    "plt.xlabel('')\n",
    "p.get_legend().remove()\n",
    "##Save\n",
    "fileDir = './ExportFigures/'\n",
    "ipynbName = '220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-Ridge_'\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": [
    "#Plot sex-specific version as a single figure\n",
    "tempD1 = {'Standard measures':standBMI_F_R2,\n",
    "          'Metabolomics':metBMI_F_R2, 'Proteomics':protBMI_F_R2,\n",
    "          'Clinical labs':chemBMI_F_R2, 'Combined omics':combiBMI_F_R2}\n",
    "tempD2 = {'Standard measures':standBMI_M_R2,\n",
    "          'Metabolomics':metBMI_M_R2, 'Proteomics':protBMI_M_R2,\n",
    "          'Clinical labs':chemBMI_M_R2, 'Combined omics':combiBMI_M_R2}\n",
    "tempD1 = {'Female':tempD1, 'Male':tempD2}\n",
    "tempD2 = {'Standard measures':'0.5',\n",
    "          'Metabolomics':'b', 'Proteomics':'r',\n",
    "          'Clinical labs':'g', 'Combined omics':'m'}\n",
    "\n",
    "#Prepare DFs\n",
    "tempD3 = {}\n",
    "tempD = {}\n",
    "for sex in tempD1.keys():\n",
    "    #Prepare DF\n",
    "    tempDF = pd.DataFrame(tempD1[sex])\n",
    "    tempD3[sex] = tempDF\n",
    "    print(sex)\n",
    "    display(tempDF.describe())\n",
    "    \n",
    "    #Statistical tests\n",
    "    control = list(tempD2.keys())[0]\n",
    "    tempL = list(tempD2.keys())[1:]\n",
    "    tempDF1 = pd.DataFrame(columns=['Control', 'Contrast', 'Control_N', 'contrast_N', 'DoF', 'tStat', 'Pval'])\n",
    "    for contrast in tempL:\n",
    "        tempS1 = tempDF[control]\n",
    "        tempS2 = tempDF[contrast]\n",
    "        #Two-sided Welch's t-test\n",
    "        tstat, pval, dof = weightstats.ttest_ind(tempS1, tempS2,\n",
    "                                                 alternative='two-sided', usevar='unequal')\n",
    "        size1 = len(tempS1)\n",
    "        size2 = len(tempS2)\n",
    "        tempDF1.loc[contrast+'-vs-'+control] = [control, contrast, size1, size2, dof, tstat, pval]\n",
    "    ##P-value adjustment (within sex) by using Benjamini–Hochberg method\n",
    "    tempDF1['AdjPval_sex'] = multi.multipletests(tempDF1['Pval'], alpha=0.05, method='fdr_bh',\n",
    "                                                 is_sorted=False, returnsorted=False)[1]\n",
    "    tempDF1['Sex'] = sex\n",
    "    tempD[sex] = 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.index.rename('ComparisonLabel', inplace=True)\n",
    "tempDF1 = tempDF1.reset_index().set_index(['Sex', 'ComparisonLabel'])\n",
    "display(tempDF1)\n",
    "##Save\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-Ridge_'\n",
    "fileName = 'R2-comparison-FemaleMale.tsv'\n",
    "tempDF1.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)\n",
    "\n",
    "#Plot\n",
    "axis_ymin = 0.0\n",
    "axis_ymax = 1.2\n",
    "ymin = 0.0\n",
    "ymax = 0.8\n",
    "yinter = 0.2\n",
    "aline_ymin = 0.8\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=(3.5*len(tempD3), 4), sharex=True, sharey=True)\n",
    "for ax_i, ax in enumerate(axes.flat):\n",
    "    sex = list(tempD3.keys())[ax_i]\n",
    "    tempDF = tempD3[sex]\n",
    "    tempDF = tempDF.melt(var_name='Category', value_name='R2', value_vars=tempDF.columns.tolist())\n",
    "    sns.barplot(data=tempDF, y='R2', x='Category', order=tempD2.keys(),\n",
    "                hue='Category', hue_order=tempD2.keys(), dodge=False, palette=tempD2,\n",
    "                ci=95, capsize=0.4, errwidth=1.5, errcolor='black', edgecolor='black', ax=ax)\n",
    "    sns.stripplot(data=tempDF, y='R2', x='Category', order=tempD2.keys(),\n",
    "                  hue='Category', hue_order=tempD2.keys(), dodge=False, jitter=0.3,\n",
    "                  size=5, edgecolor='black', color='gray', linewidth=1, alpha=0.4, ax=ax)\n",
    "    #P-value annotation\n",
    "    tempDF2 = tempDF1.loc[sex]#MultiIndex\n",
    "    for row_i in range(len(tempDF2)):\n",
    "        #Control\n",
    "        group_0 = tempDF2['Control'].iloc[row_i]\n",
    "        xcoord_0 = list(tempD2.keys()).index(group_0)\n",
    "        #Contrast\n",
    "        group_1 = tempDF2['Contrast'].iloc[row_i]\n",
    "        xcoord_1 = list(tempD2.keys()).index(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 = tempDF2['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(sex, {'fontsize':'large'})\n",
    "    #Legend\n",
    "    ax.get_legend().remove()\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 = '220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-Ridge_'\n",
    "fileName = 'R2-comparison-FemaleMale.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": [
    "### 7-2. Measured vs. Predicted"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> In this version, P-values from Pearson's correlation tests are adjusted across categories.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "tempD1 = {'Standard measures':standBMI_B,\n",
    "          'Metabolomics':metBMI_B, 'Proteomics':protBMI_B,\n",
    "          'Clinical labs':chemBMI_B, 'Combined omics':combiBMI_B}\n",
    "tempD2 = {'Standard measures':'BaseStandBMI',\n",
    "          'Metabolomics':'BaseMetBMI', 'Proteomics':'BaseProtBMI',\n",
    "          'Clinical labs':'BaseChemBMI', 'Combined omics':'BaseCombiBMI'}\n",
    "tempD3 = {'Standard measures':'0.5',\n",
    "          'Metabolomics':'b', 'Proteomics':'r',\n",
    "          'Clinical labs':'g', 'Combined omics':'m'}\n",
    "yvar='BaseBMI'\n",
    "axis_label = 'BMI [kg m'+r'$^{-2}$'+']'\n",
    "\n",
    "#Statistical tests\n",
    "tempDF1 = pd.DataFrame(columns=['N', 'DoF', 'Pearson_r', 'Pval'])\n",
    "for category in tempD1.keys():\n",
    "    tempDF = tempD1[category]\n",
    "    xvar = tempD2[category]\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[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 = '220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-Ridge_'\n",
    "fileName = 'regplot-comparison-BothSex.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=1, ncols=len(tempD3),\n",
    "                         figsize=(3.5*len(tempD3), 3.5+1), sharex=True, sharey=True)\n",
    "axis_xymin = 7.5\n",
    "axis_xymax = 62.5\n",
    "xymin = 10\n",
    "xymax = 60\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",
    "    category = list(tempD1.keys())[ax_i]\n",
    "    #Prepare DF\n",
    "    tempDF = tempD1[category]\n",
    "    xvar = tempD2[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",
    "    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",
    "    if pval==1.0:\n",
    "        pval_text = '1.0'\n",
    "    elif pval==0.0:\n",
    "        pval_text = '0.0'\n",
    "    else:\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",
    "    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",
    "    #Axis setting\n",
    "    if ax_i%len(tempD3)==0:\n",
    "        plt.setp(ax, xlabel='', ylabel='Measured '+axis_label)\n",
    "    elif ax_i==np.median(range(len(tempD3))):\n",
    "        plt.setp(ax, xlabel='Predicted '+axis_label, ylabel='')\n",
    "        plt.setp(ax.get_yticklabels(), visible=False)\n",
    "    else:\n",
    "        plt.setp(ax, xlabel='', 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",
    "#fig.text(x=0.525, y=0.0175,#Manual adjustment\n",
    "#         s='Predicted '+axis_label, fontsize='medium',\n",
    "#         verticalalignment='top', horizontalalignment='center')\n",
    "#fig.text(x=0.0125, y=0.5,#Manual adjustment\n",
    "#         s='Measured '+axis_label, fontsize='medium',\n",
    "#         verticalalignment='center', horizontalalignment='right', rotation='vertical')\n",
    "##Save\n",
    "fileDir = './ExportFigures/'\n",
    "ipynbName = '220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-Ridge_'\n",
    "fileName = 'regplot-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": [
    "## 8. Performance transition vs. critical variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> If a critical variable for the Ridge models is eliminated from the input variables, performance should be dropped.  \n",
    "> —> Repeat Ridge while dropping the top variable: the variable that was retained across 10 models and had the highest absolute value of the mean beta-coefficient. Of note, the former selection is dummy, because almost all variables are retained across 10 models in Ridge.  \n",
    "> ***–> In visualization, highlight the original top 30 variables by background color; otherwise, all dropped variables will be highlighted.***  \n",
    ">\n",
    "> Due to time consuming calculation, skip female/male model version.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 8-1. Metabolomics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Note that, in contrast to LassoCV and similarly to LinearRegression, the current RidgeCV seems to produce 2D ndarray although target y is 1D.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Both sex model\n",
    "niterations = 100\n",
    "sex = 'BothSex'\n",
    "tempDF1 = metDF_B#Standardized independent variables\n",
    "tempDF2 = bmiDF_B#Not-standardized dependent variable and info about testing set\n",
    "yvar = 'BaseBMI'\n",
    "nmodels = len(tempDF2['Testing'].unique())#Reset just in case\n",
    "ncvs = 10\n",
    "model = RidgeCV(alphas=np.geomspace(1, 1000, num=200), fit_intercept=True,\n",
    "                normalize=False, gcv_mode='auto', cv=ncvs)\n",
    "yvar_model = 'BaseMetBMI'\n",
    "\n",
    "#Perform Ridge while dropping the top retained variable\n",
    "print('Start:', time.ctime(time.time()))\n",
    "dropL = []#Initialize\n",
    "tempL = ['Model_'+str(model_k+1).zfill(2) for model_k in range(nmodels)]\n",
    "scoreDF = pd.DataFrame(columns=tempL).astype('float64')#For the coefficient of determination R2\n",
    "scoreDF.index.rename('Iteration', inplace=True)\n",
    "for iter_n in range(niterations+1):\n",
    "    print('Iteration', iter_n)\n",
    "    t_start = time.time()\n",
    "    #Drop the top variables\n",
    "    tempDF = tempDF1.drop(columns=dropL)\n",
    "    \n",
    "    #Perform Ridge\n",
    "    bcoefDF = pd.DataFrame(index=tempDF.columns).astype('float64')#For beta-coefficients\n",
    "    bcoefDF.index.rename('Variable', inplace=True)\n",
    "    scoreL = []#For the coefficient of determination R2\n",
    "    for model_k in range(nmodels):\n",
    "        #Prepare training and testing (hold-out) datasets in model k\n",
    "        model_n = 'Model_'+str(model_k+1).zfill(2)\n",
    "        yDF_train = tempDF2.loc[tempDF2['Testing']!=model_n]\n",
    "        yDF_test = tempDF2.loc[tempDF2['Testing']==model_n]\n",
    "        xDF_train = tempDF.loc[yDF_train.index.tolist()]\n",
    "        xDF_test = tempDF.loc[yDF_test.index.tolist()]\n",
    "        #Retrieve the dependent variable\n",
    "        yDF_train = pd.DataFrame(yDF_train['log_'+yvar])#Not Series but DF\n",
    "        yDF_test = pd.DataFrame(yDF_test['log_'+yvar])#Not Series but DF\n",
    "        #Fitting model with cross-validation using training dataset (i.e., internal training/validation datasets)\n",
    "        model.fit(xDF_train, yDF_train)\n",
    "        #Save parameters\n",
    "        bcoefDF[model_n] = model.coef_.flatten()#weighted vector\n",
    "        #Evaluation with testing (hold-out) dataset\n",
    "        scoreL.append(model.score(xDF_test, yDF_test))\n",
    "    \n",
    "    #Prediction accuracy\n",
    "    print(' - Out-of-sample R2 [Mean ± SEM]:',\n",
    "          np.mean(scoreL), '±', np.std(scoreL, ddof=1)/np.sqrt(len(scoreL)))\n",
    "    if np.mean(scoreL)<0:\n",
    "        print('-> Finish')\n",
    "        break\n",
    "    \n",
    "    #Add the scores to scoreDF\n",
    "    scoreDF.loc['Iteration_'+str(iter_n).zfill(3)] = scoreL\n",
    "    \n",
    "    #Clean beta-coefficient dataframe\n",
    "    count = len(bcoefDF)#Input variables\n",
    "    tempDF = bcoefDF#Just rename to reuse the same code without modification\n",
    "    ##Summarize\n",
    "    tempL1 = []\n",
    "    tempL3 = []\n",
    "    for row_n in tempDF.index.tolist():\n",
    "        tempL1.append(tempDF.loc[row_n].mean())\n",
    "        tempL3.append((tempDF.loc[row_n]==0.0).astype('int64').sum())\n",
    "    tempDF['Mean'] = tempL1\n",
    "    tempDF['nZeros'] = tempL3\n",
    "    \n",
    "    #Identify the top variable\n",
    "    ##Filter the variables with non-zero beta-coefficient in all 10 models\n",
    "    tempDF = tempDF.loc[tempDF['nZeros']==0]\n",
    "    print(' - Variables with non-zero beta-coefficient:', len(tempDF),\n",
    "          '(', len(tempDF)/count*100, '%)')\n",
    "    ##Identify the variable having the highest absolute value of the mean beta-coefficient\n",
    "    tempS = np.abs(tempDF['Mean'])\n",
    "    tempS = tempS.sort_values(ascending=False)\n",
    "    print(' - Top variable:', tempS.index[0])\n",
    "    ##Update the top variable list\n",
    "    dropL.append(tempS.index[0])\n",
    "    \n",
    "    t_elapsed = time.time() - t_start\n",
    "    print(' - Elapsed time:', round(t_elapsed//60), 'min', round(t_elapsed%60, 1), 'sec')\n",
    "print('Finish:', time.ctime(time.time()))\n",
    "\n",
    "#Save just in case for a connection timeout error\n",
    "scoreDF['TopVariable'] = dropL\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-Ridge_'\n",
    "fileName = yvar_model.replace('Base', '')+'-'+sex+'-R2transition.tsv'\n",
    "scoreDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Visualize performance transition\n",
    "yvar_model = 'MetBMI'\n",
    "tempDF1 = metBMI_B_bcoefs#All input variables with the beta-coefficient summary\n",
    "sex = 'BothSex'\n",
    "model_color = 'b'\n",
    "\n",
    "#Re-load scoreDF just in case for a connection timeout error\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-Ridge_'\n",
    "fileName = yvar_model+'-'+sex+'-R2transition.tsv'\n",
    "tempDF2 = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t').set_index('Iteration')\n",
    "\n",
    "#Prepare DF\n",
    "tempDF2['IterNum'] = list(range(len(tempDF2)))\n",
    "tempDF = tempDF2.drop(columns=['TopVariable'])\n",
    "tempDF = tempDF.reset_index().melt(var_name='Model', value_name='R2', id_vars=['Iteration', 'IterNum'])\n",
    "\n",
    "#Retrieve the original top 30 variables\n",
    "tempS = np.abs(tempDF1['Mean'])\n",
    "tempS = tempS.sort_values(ascending=False)\n",
    "tempL = tempS.iloc[:30].index.tolist()\n",
    "tempDF2 = tempDF2.loc[tempDF2['TopVariable'].isin(tempL)]\n",
    "tempDF2 = tempDF2.loc[:, ~tempDF2.columns.str.contains('Model_')]\n",
    "print('Original top 30 variables:', len(tempDF2))\n",
    "display(tempDF2)\n",
    "\n",
    "#Plot\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "plt.figure(figsize=(5, 3))\n",
    "p = sns.lineplot(data=tempDF, x='IterNum', y='R2', estimator='mean', ci=95, **{'color':'black'})\n",
    "p.set(xlim=(0, 100), ylim=(0, 0.8), yticks=np.arange(0, 0.81, 0.2))\n",
    "sns.despine()\n",
    "plt.xlabel('Iteration number')\n",
    "plt.ylabel('Out-of-sample '+r'$R^2$')\n",
    "for row_i in range(len(tempDF2)):\n",
    "    iternum = tempDF2['IterNum'].iloc[row_i]\n",
    "    p.axvspan(xmin=iternum, xmax=iternum+1, facecolor=model_color, alpha=0.4, zorder=0)\n",
    "plt.margins(0, 0, tight=True)\n",
    "##Save\n",
    "fileDir = './ExportFigures/'\n",
    "ipynbName = '220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-Ridge_'\n",
    "fileName = yvar_model+'-'+sex+'-R2transition.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": [
    "### 8-2. Proteomics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Note that, in contrast to LassoCV and similarly to LinearRegression, the current RidgeCV seems to produce 2D ndarray although target y is 1D.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Both sex model\n",
    "niterations = 100\n",
    "sex = 'BothSex'\n",
    "tempDF1 = protDF_B#Standardized independent variables\n",
    "tempDF2 = bmiDF_B#Not-standardized dependent variable and info about testing set\n",
    "yvar = 'BaseBMI'\n",
    "nmodels = len(tempDF2['Testing'].unique())#Reset just in case\n",
    "ncvs = 10\n",
    "model = RidgeCV(alphas=np.geomspace(1, 1000, num=200), fit_intercept=True,\n",
    "                normalize=False, gcv_mode='auto', cv=ncvs)\n",
    "yvar_model = 'BaseProtBMI'\n",
    "\n",
    "#Perform Ridge while dropping the top retained variable\n",
    "print('Start:', time.ctime(time.time()))\n",
    "dropL = []#Initialize\n",
    "tempL = ['Model_'+str(model_k+1).zfill(2) for model_k in range(nmodels)]\n",
    "scoreDF = pd.DataFrame(columns=tempL).astype('float64')#For the coefficient of determination R2\n",
    "scoreDF.index.rename('Iteration', inplace=True)\n",
    "for iter_n in range(niterations+1):\n",
    "    print('Iteration', iter_n)\n",
    "    t_start = time.time()\n",
    "    #Drop the top variables\n",
    "    tempDF = tempDF1.drop(columns=dropL)\n",
    "    \n",
    "    #Perform Ridge\n",
    "    bcoefDF = pd.DataFrame(index=tempDF.columns).astype('float64')#For beta-coefficients\n",
    "    bcoefDF.index.rename('Variable', inplace=True)\n",
    "    scoreL = []#For the coefficient of determination R2\n",
    "    for model_k in range(nmodels):\n",
    "        #Prepare training and testing (hold-out) datasets in model k\n",
    "        model_n = 'Model_'+str(model_k+1).zfill(2)\n",
    "        yDF_train = tempDF2.loc[tempDF2['Testing']!=model_n]\n",
    "        yDF_test = tempDF2.loc[tempDF2['Testing']==model_n]\n",
    "        xDF_train = tempDF.loc[yDF_train.index.tolist()]\n",
    "        xDF_test = tempDF.loc[yDF_test.index.tolist()]\n",
    "        #Retrieve the dependent variable\n",
    "        yDF_train = pd.DataFrame(yDF_train['log_'+yvar])#Not Series but DF\n",
    "        yDF_test = pd.DataFrame(yDF_test['log_'+yvar])#Not Series but DF\n",
    "        #Fitting model with cross-validation using training dataset (i.e., internal training/validation datasets)\n",
    "        model.fit(xDF_train, yDF_train)\n",
    "        #Save parameters\n",
    "        bcoefDF[model_n] = model.coef_.flatten()#weighted vector\n",
    "        #Evaluation with testing (hold-out) dataset\n",
    "        scoreL.append(model.score(xDF_test, yDF_test))\n",
    "    \n",
    "    #Prediction accuracy\n",
    "    print(' - Out-of-sample R2 [Mean ± SEM]:',\n",
    "          np.mean(scoreL), '±', np.std(scoreL, ddof=1)/np.sqrt(len(scoreL)))\n",
    "    if np.mean(scoreL)<0:\n",
    "        print('-> Finish')\n",
    "        break\n",
    "    \n",
    "    #Add the scores to scoreDF\n",
    "    scoreDF.loc['Iteration_'+str(iter_n).zfill(3)] = scoreL\n",
    "    \n",
    "    #Clean beta-coefficient dataframe\n",
    "    count = len(bcoefDF)#Input variables\n",
    "    tempDF = bcoefDF#Just rename to reuse the same code without modification\n",
    "    ##Summarize\n",
    "    tempL1 = []\n",
    "    tempL3 = []\n",
    "    for row_n in tempDF.index.tolist():\n",
    "        tempL1.append(tempDF.loc[row_n].mean())\n",
    "        tempL3.append((tempDF.loc[row_n]==0.0).astype('int64').sum())\n",
    "    tempDF['Mean'] = tempL1\n",
    "    tempDF['nZeros'] = tempL3\n",
    "    \n",
    "    #Identify the top variable\n",
    "    ##Filter the variables with non-zero beta-coefficient in all 10 models\n",
    "    tempDF = tempDF.loc[tempDF['nZeros']==0]\n",
    "    print(' - Variables with non-zero beta-coefficient:', len(tempDF),\n",
    "          '(', len(tempDF)/count*100, '%)')\n",
    "    ##Identify the variable having the highest absolute value of the mean beta-coefficient\n",
    "    tempS = np.abs(tempDF['Mean'])\n",
    "    tempS = tempS.sort_values(ascending=False)\n",
    "    print(' - Top variable:', tempS.index[0])\n",
    "    ##Update the top variable list\n",
    "    dropL.append(tempS.index[0])\n",
    "    \n",
    "    t_elapsed = time.time() - t_start\n",
    "    print(' - Elapsed time:', round(t_elapsed//60), 'min', round(t_elapsed%60, 1), 'sec')\n",
    "print('Finish:', time.ctime(time.time()))\n",
    "\n",
    "#Save just in case for a connection timeout error\n",
    "scoreDF['TopVariable'] = dropL\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-Ridge_'\n",
    "fileName = yvar_model.replace('Base', '')+'-'+sex+'-R2transition.tsv'\n",
    "scoreDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Visualize performance transition\n",
    "yvar_model = 'ProtBMI'\n",
    "tempDF1 = protBMI_B_bcoefs#All input variables with the beta-coefficient summary\n",
    "sex = 'BothSex'\n",
    "model_color = 'r'\n",
    "\n",
    "#Re-load scoreDF just in case for a connection timeout error\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-Ridge_'\n",
    "fileName = yvar_model+'-'+sex+'-R2transition.tsv'\n",
    "tempDF2 = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t').set_index('Iteration')\n",
    "\n",
    "#Prepare DF\n",
    "tempDF2['IterNum'] = list(range(len(tempDF2)))\n",
    "tempDF = tempDF2.drop(columns=['TopVariable'])\n",
    "tempDF = tempDF.reset_index().melt(var_name='Model', value_name='R2', id_vars=['Iteration', 'IterNum'])\n",
    "\n",
    "#Retrieve the original top 30 variables\n",
    "tempS = np.abs(tempDF1['Mean'])\n",
    "tempS = tempS.sort_values(ascending=False)\n",
    "tempL = tempS.iloc[:30].index.tolist()\n",
    "tempDF2 = tempDF2.loc[tempDF2['TopVariable'].isin(tempL)]\n",
    "tempDF2 = tempDF2.loc[:, ~tempDF2.columns.str.contains('Model_')]\n",
    "print('Original top 30 variables:', len(tempDF2))\n",
    "display(tempDF2)\n",
    "\n",
    "#Plot\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "plt.figure(figsize=(5, 3))\n",
    "p = sns.lineplot(data=tempDF, x='IterNum', y='R2', estimator='mean', ci=95, **{'color':'black'})\n",
    "p.set(xlim=(0, 100), ylim=(0, 0.8), yticks=np.arange(0, 0.81, 0.2))\n",
    "sns.despine()\n",
    "plt.xlabel('Iteration number')\n",
    "plt.ylabel('Out-of-sample '+r'$R^2$')\n",
    "for row_i in range(len(tempDF2)):\n",
    "    iternum = tempDF2['IterNum'].iloc[row_i]\n",
    "    p.axvspan(xmin=iternum, xmax=iternum+1, facecolor=model_color, alpha=0.4, zorder=0)\n",
    "plt.margins(0, 0, tight=True)\n",
    "##Save\n",
    "fileDir = './ExportFigures/'\n",
    "ipynbName = '220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-Ridge_'\n",
    "fileName = yvar_model+'-'+sex+'-R2transition.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": [
    "### 8-3. Clinical labs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Note that, in contrast to LassoCV and similarly to LinearRegression, the current RidgeCV seems to produce 2D ndarray although target y is 1D.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Both sex model\n",
    "niterations = 100\n",
    "sex = 'BothSex'\n",
    "tempDF1 = chemDF_B#Standardized independent variables\n",
    "tempDF2 = bmiDF_B#Not-standardized dependent variable and info about testing set\n",
    "yvar = 'BaseBMI'\n",
    "nmodels = len(tempDF2['Testing'].unique())#Reset just in case\n",
    "ncvs = 10\n",
    "model = RidgeCV(alphas=np.geomspace(1, 1000, num=200), fit_intercept=True,\n",
    "                normalize=False, gcv_mode='auto', cv=ncvs)\n",
    "yvar_model = 'BaseChemBMI'\n",
    "\n",
    "#Perform Ridge while dropping the top retained variable\n",
    "print('Start:', time.ctime(time.time()))\n",
    "dropL = []#Initialize\n",
    "tempL = ['Model_'+str(model_k+1).zfill(2) for model_k in range(nmodels)]\n",
    "scoreDF = pd.DataFrame(columns=tempL).astype('float64')#For the coefficient of determination R2\n",
    "scoreDF.index.rename('Iteration', inplace=True)\n",
    "for iter_n in range(niterations+1):\n",
    "    print('Iteration', iter_n)\n",
    "    t_start = time.time()\n",
    "    #Drop the top variables\n",
    "    tempDF = tempDF1.drop(columns=dropL)\n",
    "    \n",
    "    #Perform Ridge\n",
    "    bcoefDF = pd.DataFrame(index=tempDF.columns).astype('float64')#For beta-coefficients\n",
    "    bcoefDF.index.rename('Variable', inplace=True)\n",
    "    scoreL = []#For the coefficient of determination R2\n",
    "    for model_k in range(nmodels):\n",
    "        #Prepare training and testing (hold-out) datasets in model k\n",
    "        model_n = 'Model_'+str(model_k+1).zfill(2)\n",
    "        yDF_train = tempDF2.loc[tempDF2['Testing']!=model_n]\n",
    "        yDF_test = tempDF2.loc[tempDF2['Testing']==model_n]\n",
    "        xDF_train = tempDF.loc[yDF_train.index.tolist()]\n",
    "        xDF_test = tempDF.loc[yDF_test.index.tolist()]\n",
    "        #Retrieve the dependent variable\n",
    "        yDF_train = pd.DataFrame(yDF_train['log_'+yvar])#Not Series but DF\n",
    "        yDF_test = pd.DataFrame(yDF_test['log_'+yvar])#Not Series but DF\n",
    "        #Fitting model with cross-validation using training dataset (i.e., internal training/validation datasets)\n",
    "        model.fit(xDF_train, yDF_train)\n",
    "        #Save parameters\n",
    "        bcoefDF[model_n] = model.coef_.flatten()#weighted vector\n",
    "        #Evaluation with testing (hold-out) dataset\n",
    "        scoreL.append(model.score(xDF_test, yDF_test))\n",
    "    \n",
    "    #Prediction accuracy\n",
    "    print(' - Out-of-sample R2 [Mean ± SEM]:',\n",
    "          np.mean(scoreL), '±', np.std(scoreL, ddof=1)/np.sqrt(len(scoreL)))\n",
    "    if np.mean(scoreL)<0:\n",
    "        print('-> Finish')\n",
    "        break\n",
    "    \n",
    "    #Add the scores to scoreDF\n",
    "    scoreDF.loc['Iteration_'+str(iter_n).zfill(3)] = scoreL\n",
    "    \n",
    "    #Clean beta-coefficient dataframe\n",
    "    count = len(bcoefDF)#Input variables\n",
    "    tempDF = bcoefDF#Just rename to reuse the same code without modification\n",
    "    ##Summarize\n",
    "    tempL1 = []\n",
    "    tempL3 = []\n",
    "    for row_n in tempDF.index.tolist():\n",
    "        tempL1.append(tempDF.loc[row_n].mean())\n",
    "        tempL3.append((tempDF.loc[row_n]==0.0).astype('int64').sum())\n",
    "    tempDF['Mean'] = tempL1\n",
    "    tempDF['nZeros'] = tempL3\n",
    "    \n",
    "    #Identify the top variable\n",
    "    ##Filter the variables with non-zero beta-coefficient in all 10 models\n",
    "    tempDF = tempDF.loc[tempDF['nZeros']==0]\n",
    "    print(' - Variables with non-zero beta-coefficient:', len(tempDF),\n",
    "          '(', len(tempDF)/count*100, '%)')\n",
    "    ##Identify the variable having the highest absolute value of the mean beta-coefficient\n",
    "    tempS = np.abs(tempDF['Mean'])\n",
    "    tempS = tempS.sort_values(ascending=False)\n",
    "    print(' - Top variable:', tempS.index[0])\n",
    "    ##Update the top variable list\n",
    "    dropL.append(tempS.index[0])\n",
    "    \n",
    "    t_elapsed = time.time() - t_start\n",
    "    print(' - Elapsed time:', round(t_elapsed//60), 'min', round(t_elapsed%60, 1), 'sec')\n",
    "print('Finish:', time.ctime(time.time()))\n",
    "\n",
    "#Save just in case for a connection timeout error\n",
    "scoreDF['TopVariable'] = dropL\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-Ridge_'\n",
    "fileName = yvar_model.replace('Base', '')+'-'+sex+'-R2transition.tsv'\n",
    "scoreDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Visualize performance transition\n",
    "yvar_model = 'ChemBMI'\n",
    "tempDF1 = chemBMI_B_bcoefs#All input variables with the beta-coefficient summary\n",
    "sex = 'BothSex'\n",
    "model_color = 'g'\n",
    "\n",
    "#Re-load scoreDF just in case for a connection timeout error\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-Ridge_'\n",
    "fileName = yvar_model+'-'+sex+'-R2transition.tsv'\n",
    "tempDF2 = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t').set_index('Iteration')\n",
    "\n",
    "#Prepare DF\n",
    "tempDF2['IterNum'] = list(range(len(tempDF2)))\n",
    "tempDF = tempDF2.drop(columns=['TopVariable'])\n",
    "tempDF = tempDF.reset_index().melt(var_name='Model', value_name='R2', id_vars=['Iteration', 'IterNum'])\n",
    "\n",
    "#Retrieve the original top 30 variables\n",
    "tempS = np.abs(tempDF1['Mean'])\n",
    "tempS = tempS.sort_values(ascending=False)\n",
    "tempL = tempS.iloc[:30].index.tolist()\n",
    "tempDF2 = tempDF2.loc[tempDF2['TopVariable'].isin(tempL)]\n",
    "tempDF2 = tempDF2.loc[:, ~tempDF2.columns.str.contains('Model_')]\n",
    "print('Original top 30 variables:', len(tempDF2))\n",
    "display(tempDF2)\n",
    "\n",
    "#Plot\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "plt.figure(figsize=(5, 3))\n",
    "p = sns.lineplot(data=tempDF, x='IterNum', y='R2', estimator='mean', ci=95, **{'color':'black'})\n",
    "p.set(xlim=(0, 100), ylim=(0, 0.8), yticks=np.arange(0, 0.81, 0.2))\n",
    "sns.despine()\n",
    "plt.xlabel('Iteration number')\n",
    "plt.ylabel('Out-of-sample '+r'$R^2$')\n",
    "for row_i in range(len(tempDF2)):\n",
    "    iternum = tempDF2['IterNum'].iloc[row_i]\n",
    "    p.axvspan(xmin=iternum, xmax=iternum+1, facecolor=model_color, alpha=0.4, zorder=0)\n",
    "plt.margins(0, 0, tight=True)\n",
    "##Save\n",
    "fileDir = './ExportFigures/'\n",
    "ipynbName = '220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-Ridge_'\n",
    "fileName = yvar_model+'-'+sex+'-R2transition.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": [
    "### 8-4. Combined omics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Note that, in contrast to LassoCV and similarly to LinearRegression, the current RidgeCV seems to produce 2D ndarray although target y is 1D.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Both sex model\n",
    "niterations = 100\n",
    "sex = 'BothSex'\n",
    "tempDF1 = combiDF_B#Standardized independent variables\n",
    "tempDF2 = bmiDF_B#Not-standardized dependent variable and info about testing set\n",
    "yvar = 'BaseBMI'\n",
    "nmodels = len(tempDF2['Testing'].unique())#Reset just in case\n",
    "ncvs = 10\n",
    "model = RidgeCV(alphas=np.geomspace(1, 1000, num=200), fit_intercept=True,\n",
    "                normalize=False, gcv_mode='auto', cv=ncvs)\n",
    "yvar_model = 'BaseCombiBMI'\n",
    "\n",
    "#Perform Ridge while dropping the top retained variable\n",
    "print('Start:', time.ctime(time.time()))\n",
    "dropL = []#Initialize\n",
    "tempL = ['Model_'+str(model_k+1).zfill(2) for model_k in range(nmodels)]\n",
    "scoreDF = pd.DataFrame(columns=tempL).astype('float64')#For the coefficient of determination R2\n",
    "scoreDF.index.rename('Iteration', inplace=True)\n",
    "for iter_n in range(niterations+1):\n",
    "    print('Iteration', iter_n)\n",
    "    t_start = time.time()\n",
    "    #Drop the top variables\n",
    "    tempDF = tempDF1.drop(columns=dropL)\n",
    "    \n",
    "    #Perform Ridge\n",
    "    bcoefDF = pd.DataFrame(index=tempDF.columns).astype('float64')#For beta-coefficients\n",
    "    bcoefDF.index.rename('Variable', inplace=True)\n",
    "    scoreL = []#For the coefficient of determination R2\n",
    "    for model_k in range(nmodels):\n",
    "        #Prepare training and testing (hold-out) datasets in model k\n",
    "        model_n = 'Model_'+str(model_k+1).zfill(2)\n",
    "        yDF_train = tempDF2.loc[tempDF2['Testing']!=model_n]\n",
    "        yDF_test = tempDF2.loc[tempDF2['Testing']==model_n]\n",
    "        xDF_train = tempDF.loc[yDF_train.index.tolist()]\n",
    "        xDF_test = tempDF.loc[yDF_test.index.tolist()]\n",
    "        #Retrieve the dependent variable\n",
    "        yDF_train = pd.DataFrame(yDF_train['log_'+yvar])#Not Series but DF\n",
    "        yDF_test = pd.DataFrame(yDF_test['log_'+yvar])#Not Series but DF\n",
    "        #Fitting model with cross-validation using training dataset (i.e., internal training/validation datasets)\n",
    "        model.fit(xDF_train, yDF_train)\n",
    "        #Save parameters\n",
    "        bcoefDF[model_n] = model.coef_.flatten()#weighted vector\n",
    "        #Evaluation with testing (hold-out) dataset\n",
    "        scoreL.append(model.score(xDF_test, yDF_test))\n",
    "    \n",
    "    #Prediction accuracy\n",
    "    print(' - Out-of-sample R2 [Mean ± SEM]:',\n",
    "          np.mean(scoreL), '±', np.std(scoreL, ddof=1)/np.sqrt(len(scoreL)))\n",
    "    if np.mean(scoreL)<0:\n",
    "        print('-> Finish')\n",
    "        break\n",
    "    \n",
    "    #Add the scores to scoreDF\n",
    "    scoreDF.loc['Iteration_'+str(iter_n).zfill(3)] = scoreL\n",
    "    \n",
    "    #Clean beta-coefficient dataframe\n",
    "    count = len(bcoefDF)#Input variables\n",
    "    tempDF = bcoefDF#Just rename to reuse the same code without modification\n",
    "    ##Summarize\n",
    "    tempL1 = []\n",
    "    tempL3 = []\n",
    "    for row_n in tempDF.index.tolist():\n",
    "        tempL1.append(tempDF.loc[row_n].mean())\n",
    "        tempL3.append((tempDF.loc[row_n]==0.0).astype('int64').sum())\n",
    "    tempDF['Mean'] = tempL1\n",
    "    tempDF['nZeros'] = tempL3\n",
    "    \n",
    "    #Identify the top variable\n",
    "    ##Filter the variables with non-zero beta-coefficient in all 10 models\n",
    "    tempDF = tempDF.loc[tempDF['nZeros']==0]\n",
    "    print(' - Variables with non-zero beta-coefficient:', len(tempDF),\n",
    "          '(', len(tempDF)/count*100, '%)')\n",
    "    ##Identify the variable having the highest absolute value of the mean beta-coefficient\n",
    "    tempS = np.abs(tempDF['Mean'])\n",
    "    tempS = tempS.sort_values(ascending=False)\n",
    "    print(' - Top variable:', tempS.index[0])\n",
    "    ##Update the top variable list\n",
    "    dropL.append(tempS.index[0])\n",
    "    \n",
    "    t_elapsed = time.time() - t_start\n",
    "    print(' - Elapsed time:', round(t_elapsed//60), 'min', round(t_elapsed%60, 1), 'sec')\n",
    "print('Finish:', time.ctime(time.time()))\n",
    "\n",
    "#Save just in case for a connection timeout error\n",
    "scoreDF['TopVariable'] = dropL\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-Ridge_'\n",
    "fileName = yvar_model.replace('Base', '')+'-'+sex+'-R2transition.tsv'\n",
    "scoreDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Visualize performance transition\n",
    "yvar_model = 'CombiBMI'\n",
    "tempDF1 = combiBMI_B_bcoefs#All input variables with the beta-coefficient summary\n",
    "sex = 'BothSex'\n",
    "model_color = 'm'\n",
    "\n",
    "#Re-load scoreDF just in case for a connection timeout error\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-Ridge_'\n",
    "fileName = yvar_model+'-'+sex+'-R2transition.tsv'\n",
    "tempDF2 = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t').set_index('Iteration')\n",
    "\n",
    "#Prepare DF\n",
    "tempDF2['IterNum'] = list(range(len(tempDF2)))\n",
    "tempDF = tempDF2.drop(columns=['TopVariable'])\n",
    "tempDF = tempDF.reset_index().melt(var_name='Model', value_name='R2', id_vars=['Iteration', 'IterNum'])\n",
    "\n",
    "#Retrieve the original top 30 variables\n",
    "tempS = np.abs(tempDF1['Mean'])\n",
    "tempS = tempS.sort_values(ascending=False)\n",
    "tempL = tempS.iloc[:30].index.tolist()\n",
    "tempDF2 = tempDF2.loc[tempDF2['TopVariable'].isin(tempL)]\n",
    "tempDF2 = tempDF2.loc[:, ~tempDF2.columns.str.contains('Model_')]\n",
    "print('Original top 30 variables:', len(tempDF2))\n",
    "display(tempDF2)\n",
    "\n",
    "#Plot\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "plt.figure(figsize=(5, 3))\n",
    "p = sns.lineplot(data=tempDF, x='IterNum', y='R2', estimator='mean', ci=95, **{'color':'black'})\n",
    "p.set(xlim=(0, 100), ylim=(0, 0.8), yticks=np.arange(0, 0.81, 0.2))\n",
    "sns.despine()\n",
    "plt.xlabel('Iteration number')\n",
    "plt.ylabel('Out-of-sample '+r'$R^2$')\n",
    "for row_i in range(len(tempDF2)):\n",
    "    iternum = tempDF2['IterNum'].iloc[row_i]\n",
    "    p.axvspan(xmin=iternum, xmax=iternum+1, facecolor=model_color, alpha=0.4, zorder=0)\n",
    "plt.margins(0, 0, tight=True)\n",
    "##Save\n",
    "fileDir = './ExportFigures/'\n",
    "ipynbName = '220827_Multiomics-BMI-NatMed1stRevision_BMI-baseline-Ridge_'\n",
    "fileName = yvar_model+'-'+sex+'-R2transition.tif'\n",
    "plt.gcf().savefig(fileDir+ipynbName+fileName, dpi=300, bbox_inches='tight', pad_inches=0.04,\n",
    "                  pil_kwargs={'compression':'tiff_lzw'})\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## — End of this notebook —"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "arivale-py3 - Python",
   "language": "python",
   "name": "conda-env-arivale-py3-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  },
  "toc-showcode": false,
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}