Code20_Network-interaction-analysis-GLM-GEE.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multiomics BMI Paper — Network Interaction Analysis with Baseline GLM and Longitudinal GEE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***by Kengo Watanabe***  \n",
    "\n",
    "This Jupyter Notebook (with Python 3 kernel) assessed the baseline analyte–analyte interactions with MetBMI using generalized linear model (GLM), and subsequently assessed the interactions between the siginificant analyte–analyte pairs and days in program using generalized estimating equations (GEE).  \n",
    "\n",
    "Input files:  \n",
    "* Arivale baseline blood omics: 210104_Biological-BMI-paper_data-cleaning-BMI-omics_baseline-\\[combiDF/metDF/protDF/chemDF\\]-without-imputation_final-cohort.tsv  \n",
    "* Arivale baseline/longitudinal MetBMI predictions: 220805_Multiomics-BMI-NatMed1stRevision_BMI-longitudinal-LASSO_MetBMI-FemaleMale.tsv  \n",
    "* Arivale time-series BMI and blood omics: 210104_Biological-BMI-paper_data-cleaning-BMI-omics_time-series-\\[bmiDF/combiDF\\]-without-imputation_final-cohort.tsv  \n",
    "\n",
    "Output figures and tables:  \n",
    "* Figure 6b, 6c  \n",
    "* Table for Supplementary Data 7  \n",
    "\n",
    "Original notebook (memo for my future tracing):  \n",
    "* dalek:\\[JupyterLab HOME\\]/220621_Multiomics-BMI-NatMedRevision/220806_Multiomics-BMI-NatMed1stRevision_LongitudinalNetworkAnalysis-ver2.ipynb  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# packages in environment at /opt/conda/envs/arivale-py3:\n",
      "#\n",
      "# Name                    Version                   Build  Channel\n",
      "_libgcc_mutex             0.1                 conda_forge    conda-forge\n",
      "_openmp_mutex             4.5                       1_gnu    conda-forge\n",
      "analytics                 0.1                      pypi_0    pypi\n",
      "argon2-cffi               21.1.0           py39h3811e60_0    conda-forge\n",
      "arivale-data-interface    0.1.0                    pypi_0    pypi\n",
      "async_generator           1.10                       py_0    conda-forge\n",
      "atk-1.0                   2.36.0               h3371d22_4    conda-forge\n",
      "attrs                     21.2.0             pyhd8ed1ab_0    conda-forge\n",
      "backcall                  0.2.0              pyh9f0ad1d_0    conda-forge\n",
      "backports                 1.0                        py_2    conda-forge\n",
      "backports.functools_lru_cache 1.6.4              pyhd8ed1ab_0    conda-forge\n",
      "biopython                 1.79             py39h3811e60_0    conda-forge\n",
      "bleach                    4.1.0              pyhd8ed1ab_0    conda-forge\n",
      "bokeh                     2.4.1            py39hf3d152e_1    conda-forge\n",
      "boto                      2.49.0                     py_0    conda-forge\n",
      "boto3                     1.19.2                   pypi_0    pypi\n",
      "botocore                  1.22.2                   pypi_0    pypi\n",
      "brotlipy                  0.7.0           py39h3811e60_1001    conda-forge\n",
      "bzip2                     1.0.8                h7f98852_4    conda-forge\n",
      "c-ares                    1.17.2               h7f98852_0    conda-forge\n",
      "ca-certificates           2022.6.15            ha878542_0    conda-forge\n",
      "cachecontrol              0.12.6                     py_0    conda-forge\n",
      "cairo                     1.16.0            h6cf1ce9_1008    conda-forge\n",
      "certifi                   2022.6.15        py39hf3d152e_0    conda-forge\n",
      "cffi                      1.14.6           py39h4bc2ebd_1    conda-forge\n",
      "chardet                   4.0.0            py39hf3d152e_1    conda-forge\n",
      "charset-normalizer        2.0.0              pyhd8ed1ab_0    conda-forge\n",
      "click                     8.0.3            py39hf3d152e_0    conda-forge\n",
      "colorama                  0.4.4              pyh9f0ad1d_0    conda-forge\n",
      "cryptography              35.0.0           py39h95dcef6_1    conda-forge\n",
      "cycler                    0.10.0                     py_2    conda-forge\n",
      "cython                    0.29.24          py39he80948d_0    conda-forge\n",
      "debugpy                   1.4.1            py39he80948d_0    conda-forge\n",
      "decorator                 4.4.2                      py_0    conda-forge\n",
      "defusedxml                0.7.1              pyhd8ed1ab_0    conda-forge\n",
      "english                   2020.7.0           pyhd8ed1ab_0    conda-forge\n",
      "entrypoints               0.3             py39hde42818_1002    conda-forge\n",
      "et_xmlfile                1.0.1                   py_1001    conda-forge\n",
      "expat                     2.4.1                h9c3ff4c_0    conda-forge\n",
      "font-ttf-dejavu-sans-mono 2.37                 hab24e00_0    conda-forge\n",
      "font-ttf-inconsolata      3.000                h77eed37_0    conda-forge\n",
      "font-ttf-source-code-pro  2.038                h77eed37_0    conda-forge\n",
      "font-ttf-ubuntu           0.83                 hab24e00_0    conda-forge\n",
      "fontconfig                2.13.1            hba837de_1005    conda-forge\n",
      "fonts-conda-ecosystem     1                             0    conda-forge\n",
      "fonts-conda-forge         1                             0    conda-forge\n",
      "freetype                  2.10.4               h0708190_1    conda-forge\n",
      "fribidi                   1.0.10               h36c2ea0_0    conda-forge\n",
      "future                    0.18.2           py39hf3d152e_3    conda-forge\n",
      "gdk-pixbuf                2.42.6               h04a7f16_0    conda-forge\n",
      "gettext                   0.19.8.1          h73d1719_1008    conda-forge\n",
      "giflib                    5.2.1                h36c2ea0_2    conda-forge\n",
      "graphite2                 1.3.13            h58526e2_1001    conda-forge\n",
      "graphviz                  2.49.1               h85b4f2f_0    conda-forge\n",
      "gtk2                      2.24.33              h539f30e_1    conda-forge\n",
      "gts                       0.7.6                h64030ff_2    conda-forge\n",
      "harfbuzz                  3.0.0                h83ec7ef_1    conda-forge\n",
      "hdmedians                 0.14.2           py39hce5d2b2_0    conda-forge\n",
      "icu                       68.1                 h58526e2_0    conda-forge\n",
      "idna                      3.1                pyhd3deb0d_0    conda-forge\n",
      "importlib-metadata        4.8.1            py39hf3d152e_0    conda-forge\n",
      "iniconfig                 1.1.1              pyh9f0ad1d_0    conda-forge\n",
      "interchange               2021.0.4           pyhd8ed1ab_0    conda-forge\n",
      "ipykernel                 6.4.2            py39hef51801_0    conda-forge\n",
      "ipython                   7.28.0           py39hef51801_0    conda-forge\n",
      "ipython_genutils          0.2.0                      py_1    conda-forge\n",
      "ipywidgets                7.6.5              pyhd8ed1ab_0    conda-forge\n",
      "jbig                      2.1               h7f98852_2003    conda-forge\n",
      "jedi                      0.18.0           py39hf3d152e_2    conda-forge\n",
      "jinja2                    3.0.2              pyhd8ed1ab_0    conda-forge\n",
      "jmespath                  0.10.0                   pypi_0    pypi\n",
      "joblib                    1.1.0              pyhd8ed1ab_0    conda-forge\n",
      "jpeg                      9d                   h36c2ea0_0    conda-forge\n",
      "jsonschema                4.1.2              pyhd8ed1ab_0    conda-forge\n",
      "jupyter_client            7.0.6              pyhd8ed1ab_0    conda-forge\n",
      "jupyter_core              4.8.1            py39hf3d152e_0    conda-forge\n",
      "jupyterlab_pygments       0.1.2              pyh9f0ad1d_0    conda-forge\n",
      "jupyterlab_widgets        1.0.2              pyhd8ed1ab_0    conda-forge\n",
      "kiwisolver                1.3.2            py39h1a9c180_0    conda-forge\n",
      "krb5                      1.19.2               hcc1bbae_2    conda-forge\n",
      "lcms2                     2.12                 hddcbb42_0    conda-forge\n",
      "ld_impl_linux-64          2.36.1               hea4e1c9_2    conda-forge\n",
      "lerc                      2.2.1                h9c3ff4c_0    conda-forge\n",
      "libblas                   3.9.0           12_linux64_openblas    conda-forge\n",
      "libcblas                  3.9.0           12_linux64_openblas    conda-forge\n",
      "libcurl                   7.79.1               h2574ce0_1    conda-forge\n",
      "libdeflate                1.7                  h7f98852_5    conda-forge\n",
      "libedit                   3.1.20191231         he28a2e2_2    conda-forge\n",
      "libev                     4.33                 h516909a_1    conda-forge\n",
      "libffi                    3.4.2                h9c3ff4c_4    conda-forge\n",
      "libgcc-ng                 11.2.0              h1d223b6_11    conda-forge\n",
      "libgd                     2.3.3                h6ad9fb6_0    conda-forge\n",
      "libgfortran-ng            11.2.0              h69a702a_11    conda-forge\n",
      "libgfortran5              11.2.0              h5c6108e_11    conda-forge\n",
      "libglib                   2.70.0               h174f98d_1    conda-forge\n",
      "libgomp                   11.2.0              h1d223b6_11    conda-forge\n",
      "libiconv                  1.16                 h516909a_0    conda-forge\n",
      "liblapack                 3.9.0           12_linux64_openblas    conda-forge\n",
      "libnghttp2                1.43.0               h812cca2_1    conda-forge\n",
      "libopenblas               0.3.18          pthreads_h8fe5266_0    conda-forge\n",
      "libpng                    1.6.37               h21135ba_2    conda-forge\n",
      "libpq                     13.3                 hd57d9b9_1    conda-forge\n",
      "librsvg                   2.52.2               hc3c00ef_0    conda-forge\n",
      "libsodium                 1.0.18               h36c2ea0_1    conda-forge\n",
      "libssh2                   1.10.0               ha56f1ee_2    conda-forge\n",
      "libstdcxx-ng              11.2.0              he4da1e4_11    conda-forge\n",
      "libtiff                   4.3.0                hf544144_1    conda-forge\n",
      "libtool                   2.4.6             h9c3ff4c_1008    conda-forge\n",
      "libuuid                   2.32.1            h7f98852_1000    conda-forge\n",
      "libwebp                   1.2.1                h3452ae3_0    conda-forge\n",
      "libwebp-base              1.2.1                h7f98852_0    conda-forge\n",
      "libxcb                    1.13              h7f98852_1003    conda-forge\n",
      "libxml2                   2.9.12               h72842e0_0    conda-forge\n",
      "libxslt                   1.1.33               h15afd5d_2    conda-forge\n",
      "libzlib                   1.2.11            h36c2ea0_1013    conda-forge\n",
      "lockfile                  0.12.2                     py_1    conda-forge\n",
      "lxml                      4.6.3            py39h107f48f_0    conda-forge\n",
      "lz4-c                     1.9.3                h9c3ff4c_1    conda-forge\n",
      "markupsafe                2.0.1            py39h3811e60_0    conda-forge\n",
      "matplotlib-base           3.4.3            py39h2fa2bec_1    conda-forge\n",
      "matplotlib-inline         0.1.3              pyhd8ed1ab_0    conda-forge\n",
      "matplotlib-venn           0.11.7             pyhd8ed1ab_0    conda-forge\n",
      "missingpy                 0.2.0                    pypi_0    pypi\n",
      "mistune                   0.8.4           py39h3811e60_1004    conda-forge\n",
      "more-itertools            8.10.0             pyhd8ed1ab_0    conda-forge\n",
      "mscorefonts               0.0.1                         3    conda-forge\n",
      "msgpack-python            1.0.2            py39h1a9c180_1    conda-forge\n",
      "natsort                   7.1.1              pyhd8ed1ab_0    conda-forge\n",
      "nbclient                  0.5.4              pyhd8ed1ab_0    conda-forge\n",
      "nbconvert                 6.2.0            py39hf3d152e_0    conda-forge\n",
      "nbformat                  5.1.3              pyhd8ed1ab_0    conda-forge\n",
      "ncurses                   6.2                  h58526e2_4    conda-forge\n",
      "neo4j-python-driver       4.3.7              pyhd8ed1ab_0    conda-forge\n",
      "neobolt                   1.7.17           py39h3811e60_2    conda-forge\n",
      "neotime                   1.7.4                      py_0    conda-forge\n",
      "nest-asyncio              1.5.1              pyhd8ed1ab_0    conda-forge\n",
      "networkx                  2.6.3              pyhd8ed1ab_0    conda-forge\n",
      "notebook                  6.4.5              pyha770c72_0    conda-forge\n",
      "numpy                     1.21.3           py39hdbf815f_0    conda-forge\n",
      "olefile                   0.46               pyh9f0ad1d_1    conda-forge\n",
      "openjpeg                  2.4.0                hb52868f_1    conda-forge\n",
      "openpyxl                  3.0.9              pyhd8ed1ab_0    conda-forge\n",
      "openssl                   1.1.1o               h166bdaf_0    conda-forge\n",
      "packaging                 21.0               pyhd8ed1ab_0    conda-forge\n",
      "pandas                    1.3.4            py39hde0f152_0    conda-forge\n",
      "pandoc                    2.15                 h7f98852_0    conda-forge\n",
      "pandocfilters             1.5.0              pyhd8ed1ab_0    conda-forge\n",
      "pango                     1.48.10              h54213e6_2    conda-forge\n",
      "parso                     0.8.2              pyhd8ed1ab_0    conda-forge\n",
      "patsy                     0.5.2              pyhd8ed1ab_0    conda-forge\n",
      "pcre                      8.45                 h9c3ff4c_0    conda-forge\n",
      "pexpect                   4.8.0              pyh9f0ad1d_2    conda-forge\n",
      "pickleshare               0.7.5           py39hde42818_1002    conda-forge\n",
      "pillow                    8.3.2            py39ha612740_0    conda-forge\n",
      "pip                       21.3.1             pyhd8ed1ab_0    conda-forge\n",
      "pixman                    0.40.0               h36c2ea0_0    conda-forge\n",
      "plotly                    5.3.1              pyhd8ed1ab_0    conda-forge\n",
      "pluggy                    1.0.0            py39hf3d152e_1    conda-forge\n",
      "prometheus_client         0.11.0             pyhd8ed1ab_0    conda-forge\n",
      "prompt-toolkit            3.0.21             pyha770c72_0    conda-forge\n",
      "prompt_toolkit            3.0.21               hd8ed1ab_0    conda-forge\n",
      "psycopg2                  2.9.1            py39h3811e60_0    conda-forge\n",
      "pthread-stubs             0.4               h36c2ea0_1001    conda-forge\n",
      "ptyprocess                0.7.0              pyhd3deb0d_0    conda-forge\n",
      "py                        1.10.0             pyhd3deb0d_0    conda-forge\n",
      "py2neo                    2021.2.3           pyhd8ed1ab_0    conda-forge\n",
      "pycparser                 2.20               pyh9f0ad1d_2    conda-forge\n",
      "pygments                  2.10.0             pyhd8ed1ab_0    conda-forge\n",
      "pygraphviz                1.7              py39h78163bd_0    conda-forge\n",
      "pyopenssl                 21.0.0             pyhd8ed1ab_0    conda-forge\n",
      "pyparsing                 2.4.7              pyh9f0ad1d_0    conda-forge\n",
      "pyrsistent                0.17.3           py39h3811e60_2    conda-forge\n",
      "pysam                     0.17.0           py39h051187c_0    bioconda\n",
      "pysocks                   1.7.1            py39hf3d152e_3    conda-forge\n",
      "pytabix                   0.1              py39h98c8e45_1    bioconda\n",
      "pytest                    6.2.5            py39hf3d152e_0    conda-forge\n",
      "python                    3.9.7           hb7a2778_3_cpython    conda-forge\n",
      "python-dateutil           2.8.2              pyhd8ed1ab_0    conda-forge\n",
      "python-graphviz           0.17               pyhaef67bd_0    conda-forge\n",
      "python_abi                3.9                      2_cp39    conda-forge\n",
      "pytz                      2021.3             pyhd8ed1ab_0    conda-forge\n",
      "pyvcf                     0.6.8           py39hde42818_1002    conda-forge\n",
      "pyyaml                    6.0              py39h3811e60_0    conda-forge\n",
      "pyzmq                     22.3.0           py39h37b5a0c_0    conda-forge\n",
      "readline                  8.1                  h46c0cb4_0    conda-forge\n",
      "requests                  2.26.0             pyhd8ed1ab_0    conda-forge\n",
      "s3transfer                0.5.0                    pypi_0    pypi\n",
      "scikit-bio                0.5.6            py39h16ac069_4    conda-forge\n",
      "scikit-learn              1.0.1            py39h7c5d8c9_0    conda-forge\n",
      "scipy                     1.7.1            py39hee8e79c_0    conda-forge\n",
      "seaborn                   0.11.2               hd8ed1ab_0    conda-forge\n",
      "seaborn-base              0.11.2             pyhd8ed1ab_0    conda-forge\n",
      "send2trash                1.8.0              pyhd8ed1ab_0    conda-forge\n",
      "setuptools                58.2.0           py39hf3d152e_0    conda-forge\n",
      "six                       1.16.0             pyh6c4a22f_0    conda-forge\n",
      "sqlite                    3.36.0               h9cd32fc_2    conda-forge\n",
      "statsmodels               0.13.0           py39hce5d2b2_0    conda-forge\n",
      "tenacity                  8.0.1              pyhd8ed1ab_0    conda-forge\n",
      "terminado                 0.12.1           py39hf3d152e_0    conda-forge\n",
      "testpath                  0.5.0              pyhd8ed1ab_0    conda-forge\n",
      "threadpoolctl             3.0.0              pyh8a188c0_0    conda-forge\n",
      "tk                        8.6.11               h27826a3_1    conda-forge\n",
      "toml                      0.10.2             pyhd8ed1ab_0    conda-forge\n",
      "tornado                   6.1              py39h3811e60_1    conda-forge\n",
      "traitlets                 5.1.0              pyhd8ed1ab_0    conda-forge\n",
      "typing_extensions         3.10.0.2           pyha770c72_0    conda-forge\n",
      "tzdata                    2021e                he74cb21_0    conda-forge\n",
      "urllib3                   1.26.7             pyhd8ed1ab_0    conda-forge\n",
      "venn                      0.1.3                    pypi_0    pypi\n",
      "wcwidth                   0.2.5              pyh9f0ad1d_2    conda-forge\n",
      "webencodings              0.5.1                      py_1    conda-forge\n",
      "wheel                     0.37.0             pyhd8ed1ab_1    conda-forge\n",
      "widgetsnbextension        3.5.1            py39hf3d152e_4    conda-forge\n",
      "wordcloud                 1.8.1                    pypi_0    pypi\n",
      "xlrd                      2.0.1              pyhd8ed1ab_3    conda-forge\n",
      "xorg-kbproto              1.0.7             h7f98852_1002    conda-forge\n",
      "xorg-libice               1.0.10               h7f98852_0    conda-forge\n",
      "xorg-libsm                1.2.3             hd9c2040_1000    conda-forge\n",
      "xorg-libx11               1.7.2                h7f98852_0    conda-forge\n",
      "xorg-libxau               1.0.9                h7f98852_0    conda-forge\n",
      "xorg-libxdmcp             1.1.3                h7f98852_0    conda-forge\n",
      "xorg-libxext              1.3.4                h7f98852_1    conda-forge\n",
      "xorg-libxrender           0.9.10            h7f98852_1003    conda-forge\n",
      "xorg-renderproto          0.11.1            h7f98852_1002    conda-forge\n",
      "xorg-xextproto            7.3.0             h7f98852_1002    conda-forge\n",
      "xorg-xproto               7.0.31            h7f98852_1007    conda-forge\n",
      "xz                        5.2.5                h516909a_1    conda-forge\n",
      "yaml                      0.2.5                h516909a_0    conda-forge\n",
      "zeromq                    4.3.4                h9c3ff4c_1    conda-forge\n",
      "zipp                      3.6.0              pyhd8ed1ab_0    conda-forge\n",
      "zlib                      1.2.11            h36c2ea0_1013    conda-forge\n",
      "zstd                      1.5.0                ha95c52a_0    conda-forge\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scipy.stats as stats\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import seaborn as sns\n",
    "#For Arial font\n",
    "#!conda install -c conda-forge -y mscorefonts\n",
    "##-> The below was also needed in matplotlib 3.4.2\n",
    "#import shutil\n",
    "#import matplotlib\n",
    "#shutil.rmtree(matplotlib.get_cachedir())\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "from IPython.display import display\n",
    "import time\n",
    "\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "from statsmodels.stats import multitest as multi\n",
    "import matplotlib.lines as mlines\n",
    "\n",
    "!conda list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Preprocessing for GLM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1-1. Prepare dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Import the baseline all omics without imputation\n",
    "fileDir = '../210104_Biological-BMI-paper/ExportData/'\n",
    "ipynbName = '210104_Biological-BMI-paper_data-cleaning-BMI-omics_'\n",
    "fileName = 'baseline-combiDF-without-imputation_final-cohort.tsv'\n",
    "tempDF = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t', dtype={'public_client_id': str})\n",
    "tempDF = tempDF.set_index('public_client_id')\n",
    "\n",
    "#Retrieve the baseline MetBMI-based class from the sex-stratified model\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220805_Multiomics-BMI-NatMed1stRevision_BMI-longitudinal-LASSO_'\n",
    "fileName = 'MetBMI-FemaleMale.tsv'\n",
    "tempDF1 = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t', dtype={'public_client_id': str})\n",
    "tempDF1 = tempDF1.sort_values(by=['public_client_id', 'days_in_program'], ascending=True)\n",
    "tempDF1 = tempDF1.drop_duplicates('public_client_id', keep='first')\n",
    "tempDF1 = tempDF1.set_index('public_client_id')\n",
    "tempDF1['log_BaseMetBMI'] = np.log(tempDF1['BaseMetBMI'])\n",
    "tempDF1 = tempDF1[['log_BaseMetBMI', 'BaseMetBMI_class']]\n",
    "tempDF = pd.merge(tempDF, tempDF1, left_index=True, right_index=True, how='left')\n",
    "\n",
    "#To apply the skewness reduction pipeline, revert the log-transformed metabolite values back to orignal scale\n",
    "covarL = ['log_BaseBMI', 'Sex', 'BaseAge', 'PC1', 'PC2', 'PC3', 'PC4', 'PC5', 'Race']\n",
    "fileDir = '../210104_Biological-BMI-paper/ExportData/'\n",
    "ipynbName = '210104_Biological-BMI-paper_data-cleaning-BMI-omics_'\n",
    "fileName = 'baseline-metDF-without-imputation_final-cohort.tsv'\n",
    "tempDF1 = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t', dtype={'public_client_id': str})\n",
    "tempDF1 = tempDF1.set_index('public_client_id')\n",
    "for col_n in tempDF1.drop(columns=covarL).columns.tolist():\n",
    "    tempDF[col_n] = np.e**tempDF[col_n]\n",
    "\n",
    "display(tempDF)\n",
    "#Update\n",
    "combiDF_base = tempDF\n",
    "covarL = ['log_BaseBMI', 'Sex', 'BaseAge', 'PC1', 'PC2', 'PC3', 'PC4', 'PC5', 'Race',\n",
    "          'log_BaseMetBMI', 'BaseMetBMI_class']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF = combiDF_base#Not copy but just rename\n",
    "print('Unique ID before elimination:', len(tempDF))\n",
    "\n",
    "#Select the cohort used in the longitudinal analysis\n",
    "month_threshold = 18\n",
    "##Select participants who have 2 or more measuremnts in BMI\n",
    "fileDir = '../210104_Biological-BMI-paper/ExportData/'\n",
    "ipynbName = '210104_Biological-BMI-paper_data-cleaning-BMI-omics_'\n",
    "fileName = 'time-series-bmiDF-without-imputation_final-cohort.tsv'\n",
    "tempDF1 = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t', dtype={'public_client_id': str})\n",
    "tempDF1 = tempDF1.set_index('KeyIndex')\n",
    "tempDF1 = tempDF1.loc[tempDF1['days_in_program'] <= 365.25/12*month_threshold]\n",
    "tempS = tempDF1['public_client_id'].value_counts()\n",
    "tempL1 = tempS.loc[tempS>=2].index.tolist()\n",
    "##Select participants who have 2 or more measuremnts in omics\n",
    "fileDir = '../210104_Biological-BMI-paper/ExportData/'\n",
    "ipynbName = '210104_Biological-BMI-paper_data-cleaning-BMI-omics_'\n",
    "fileName = 'time-series-combiDF-without-imputation_final-cohort.tsv'\n",
    "tempDF1 = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t', dtype={'public_client_id': str})\n",
    "tempDF1 = tempDF1.set_index('KeyIndex')\n",
    "tempDF1 = tempDF1.loc[tempDF1['days_in_program'] <= 365.25/12*month_threshold]\n",
    "tempS = tempDF1['public_client_id'].value_counts()\n",
    "tempL2 = tempS.loc[tempS>=2].index.tolist()\n",
    "##Select participants who have 2 or more measuremnts in both BMI and omics\n",
    "tempL = list(set(tempL1) & set(tempL2))\n",
    "print('Participants who have more than 2 measurements in both BMI and omics:', len(tempL))\n",
    "tempDF = tempDF.loc[tempDF.index.isin(tempL)]\n",
    "print('Unique ID after elimination:', len(tempDF))\n",
    "\n",
    "display(tempDF)\n",
    "#Update\n",
    "combiDF_base = tempDF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1-2. Drop pseudo numeric analytes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Some analytes in clinical labs have almost same values between participants.  \n",
    "> –> Because becoming invariant values after dropping outliers, these analytes will raise errors in GLM.  \n",
    "> –> Eliminate these in advance.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF = combiDF_base#Not copy but just rename\n",
    "print('nAnalytes before the step: ', len(tempDF.drop(columns=covarL).columns))\n",
    "t_start = time.time()\n",
    "for col_n in tempDF.drop(columns=covarL).columns.tolist():\n",
    "    #Drop category level <15 (cf.) The category level of 'race' is 14.\n",
    "    if len(tempDF[col_n].unique())<15:\n",
    "        tempDF = tempDF.drop(columns=[col_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('nAnalytes after the step: ', len(tempDF.drop(columns=covarL).columns))\n",
    "combiDF_base = tempDF#Update the original"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1-3. Replace zero/negative value with NaN"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Some analytes were measured as meaningful zeros in clinical labs. Also, analytes in proteomics could be negative values after the batch correction.  \n",
    "> –> Skip this step.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Check zero or negative values\n",
    "tempDF = combiDF_base.drop(columns=covarL)\n",
    "print('Analytes having zero or negative values')\n",
    "tempS = tempDF.min(axis=0)\n",
    "display(tempDF.describe(include='all').loc[:, tempS.loc[tempS<=0].index])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1-4. Drop outliers (>3SD from MEAN)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Check\n",
    "print('Before the step:')\n",
    "tempDF = combiDF_base.drop(columns=covarL)\n",
    "display(tempDF.describe())\n",
    "#Distribution of examples\n",
    "sns.set(style='ticks', font='Arial', context='paper')\n",
    "fig, axes = plt.subplots(nrows=2, ncols=5, figsize=(15, 6), sharex=False, sharey=False)\n",
    "for col_i, ax in enumerate(axes.flat):\n",
    "    sns.distplot(tempDF.iloc[:, col_i].dropna(), color='g', ax=ax)\n",
    "    ax.set_title(tempDF.columns[col_i], fontsize='large')\n",
    "sns.despine()\n",
    "plt.setp(axes, xlabel='', ylabel='')\n",
    "plt.setp(axes[1, :], xlabel='Raw value')\n",
    "plt.setp(axes[:, 0], ylabel='Density')\n",
    "fig.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF = combiDF_base#Not copy but just rename\n",
    "t_start = time.time()\n",
    "#Replace more than 3SD from MEAN with NaN\n",
    "for col_n in tempDF.drop(columns=covarL).columns.tolist():\n",
    "    col_mean = tempDF[col_n].mean()\n",
    "    col_sd = tempDF[col_n].std(ddof=1)#Sample standard deviation: devided with N-1\n",
    "    tempL = []\n",
    "    for row_n in tempDF.index.tolist():\n",
    "        value = tempDF.at[row_n, col_n]\n",
    "        if (value > (col_mean - 3*col_sd)) & (value < (col_mean + 3*col_sd)):\n",
    "            tempL.append(value)\n",
    "        else:#outlier\n",
    "            tempL.append(np.nan)\n",
    "    tempDF[col_n] = tempL\n",
    "t_elapsed = time.time() - t_start\n",
    "print('Elapsed time:', round(t_elapsed//60), 'min', round(t_elapsed%60, 1), 'sec')\n",
    "combiDF_base = tempDF#Update the original"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Check\n",
    "print('After the step:')\n",
    "tempDF = combiDF_base.drop(columns=covarL)\n",
    "display(tempDF.describe())\n",
    "#Distribution of examples\n",
    "sns.set(style='ticks', font='Arial', context='paper')\n",
    "fig, axes = plt.subplots(nrows=2, ncols=5, figsize=(15, 6), sharex=False, sharey=False)\n",
    "for col_i, ax in enumerate(axes.flat):\n",
    "    sns.distplot(tempDF.iloc[:, col_i].dropna(), color='g', ax=ax)\n",
    "    ax.set_title(tempDF.columns[col_i], fontsize='large')\n",
    "sns.despine()\n",
    "plt.setp(axes, xlabel='', ylabel='')\n",
    "plt.setp(axes[1, :], xlabel='Raw value')\n",
    "plt.setp(axes[:, 0], ylabel='Density')\n",
    "fig.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1-5. Drop analytes with small sample size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF = combiDF_base#Not copy but just rename\n",
    "print('nAnalytes before the step: ', len(tempDF.drop(columns=covarL).columns))\n",
    "t_start = time.time()\n",
    "for col_n in tempDF.drop(columns=covarL).columns.tolist():\n",
    "    tempS = tempDF[col_n].dropna()\n",
    "    if len(tempS) < 100:\n",
    "        #Drop this analyte as one with small sample size\n",
    "        tempDF = tempDF.drop(columns=[col_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('nAnalytes after the step: ', len(tempDF.drop(columns=covarL).columns))\n",
    "combiDF_base = tempDF#Update the original"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1-6. Transform to less skewed distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Check\n",
    "print('Before the step:')\n",
    "tempDF = combiDF_base.drop(columns=covarL)\n",
    "print('Confirm total number of NaNs in DF:', tempDF.isnull().to_numpy().sum(axis=None))\n",
    "#Skewness\n",
    "tempL = []\n",
    "for col_i in range(len(tempDF.columns)):\n",
    "    tempL.append(stats.skew(tempDF.iloc[:, col_i].dropna()))\n",
    "tempS = pd.Series(tempL, index=tempDF.columns, name='Skewness')\n",
    "display(tempS.describe())\n",
    "#Distribution and probability plot of examples\n",
    "sns.set(style='ticks', font='Arial', context='paper')\n",
    "fig, axes = plt.subplots(nrows=4, ncols=5, figsize=(15, 12), sharex=False, sharey=False)\n",
    "for ax_i, ax in enumerate(axes.flat):\n",
    "    if (ax_i//5)%2 == 0:\n",
    "        col_i = (ax_i//10)*5 + ax_i%5\n",
    "        sns.set(style='ticks', font='Arial', context='paper')\n",
    "        sns.distplot(tempDF.iloc[:, col_i].dropna(), color='g', ax=ax)\n",
    "        sns.despine()\n",
    "        ax.set_title(tempDF.columns[col_i], fontsize='large')\n",
    "        ax.set_xlabel('Raw value')\n",
    "    else:\n",
    "        sns.set(style='whitegrid', font='Arial', context='paper')\n",
    "        col_i = (ax_i//10)*5 + ax_i%5\n",
    "        stats.probplot(tempDF.iloc[:, col_i].dropna(), plot=ax)\n",
    "        ax.get_lines()[0].set_markerfacecolor('g')\n",
    "        ax.get_lines()[0].set_markeredgecolor('g')\n",
    "        ax.get_lines()[1].set_color('k')\n",
    "        ax.get_lines()[1].set_linewidth(3)\n",
    "        skewness = tempS.loc[tempDF.columns[col_i]]\n",
    "        ax.set_title('Skewness: '+str(round(skewness, 3)), fontsize='large')\n",
    "fig.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF = combiDF_base#Not copy but just rename\n",
    "tempL = tempDF.drop(columns=covarL).columns.tolist()\n",
    "\n",
    "#Prepare analyte category list to check transformation per each omics\n",
    "fileDir = '../210104_Biological-BMI-paper/ExportData/'\n",
    "ipynbName = '210104_Biological-BMI-paper_data-cleaning-BMI-omics_'\n",
    "fileName = 'baseline-metDF-without-imputation_final-cohort.tsv'\n",
    "tempDF1 = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t', dtype={'public_client_id': str})\n",
    "tempDF1 = tempDF1.set_index('public_client_id')\n",
    "tempL1 = tempDF1.loc[:, ~tempDF1.columns.isin(covarL)].columns.tolist()\n",
    "tempL1 = list(set(tempL) & set(tempL1))\n",
    "fileName = 'baseline-protDF-without-imputation_final-cohort.tsv'\n",
    "tempDF1 = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t', dtype={'public_client_id': str})\n",
    "tempDF1 = tempDF1.set_index('public_client_id')\n",
    "tempL2 = tempDF1.loc[:, ~tempDF1.columns.isin(covarL)].columns.tolist()\n",
    "tempL2 = list(set(tempL) & set(tempL2))\n",
    "fileName = 'baseline-chemDF-without-imputation_final-cohort.tsv'\n",
    "tempDF1 = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t', dtype={'public_client_id': str})\n",
    "tempDF1 = tempDF1.set_index('public_client_id')\n",
    "tempL3 = tempDF1.loc[:, ~tempDF1.columns.isin(covarL)].columns.tolist()\n",
    "tempL3 = list(set(tempL) & set(tempL3))\n",
    "tempD = {'Metabolomics':tempL1, 'Proteomics':tempL2, 'Clinical labs':tempL3}\n",
    "\n",
    "for omics in tempD.keys():\n",
    "    t_start = time.time()\n",
    "    print(omics)\n",
    "    #Transform to less skewed distribution\n",
    "    counter_r1 = 0\n",
    "    counter_r2 = 0\n",
    "    counter_r3 = 0\n",
    "    counter_l1 = 0\n",
    "    counter_l2 = 0\n",
    "    counter_l3 = 0\n",
    "    counter_l4 = 0\n",
    "    counter_l5 = 0\n",
    "    counter_l6 = 0\n",
    "    for col_n in tempD[omics]:\n",
    "        tempS = tempDF[col_n].copy(deep=True)\n",
    "        col_skew_0 = stats.skew(tempS.dropna())\n",
    "        ##Shift the distribution to avoid NaN transformation for zero or negative values\n",
    "        tempS = tempS - tempS.min() + 1\n",
    "        if col_skew_0 > 0:#Right/positive skewed distribution\n",
    "            #Transformation\n",
    "            tempS1 = np.log(tempS)#Assumed for extremely skewed\n",
    "            tempS2 = np.cbrt(tempS)#Assumed for strongly skewed\n",
    "            tempS3 = np.sqrt(tempS)#Assumed for skewed\n",
    "            #Evaluation\n",
    "            col_skew_1 = stats.skew(tempS1.dropna())\n",
    "            col_skew_2 = stats.skew(tempS2.dropna())\n",
    "            col_skew_3 = stats.skew(tempS3.dropna())\n",
    "            col_skew_min = np.min([abs(col_skew_0), abs(col_skew_1), abs(col_skew_2), abs(col_skew_3)])\n",
    "            if col_skew_min == abs(col_skew_1):\n",
    "                tempDF[col_n] = tempS1\n",
    "                counter_r1 += 1\n",
    "            elif col_skew_min == abs(col_skew_2):\n",
    "                tempDF[col_n] = tempS2\n",
    "                counter_r2 += 1\n",
    "            elif col_skew_min == abs(col_skew_3):\n",
    "                tempDF[col_n] = tempS3\n",
    "                counter_r3 += 1\n",
    "        elif col_skew_0 < 0:#Left/negative skewed distribution\n",
    "            #Version1\n",
    "            ##Transformation\n",
    "            tempS1 = np.e**tempS#Assumed for extremely skewed\n",
    "            tempS2 = tempS**3#Assumed for strongly skewed\n",
    "            tempS3 = tempS**2#Assumed for skewed\n",
    "            #Version2\n",
    "            ##Flip horizontally to transform to the right skewed\n",
    "            col_max = tempS.max()\n",
    "            col_min = tempS.min()\n",
    "            tempS = (col_max + col_min) - tempS\n",
    "            ##Transformation\n",
    "            tempS4 = np.log(tempS)#Assumed for extremely skewed\n",
    "            tempS5 = np.cbrt(tempS)#Assumed for strongly skewed\n",
    "            tempS6 = np.sqrt(tempS)#Assumed for skewed\n",
    "            #Evaluation for both version 1 and 2\n",
    "            col_skew_1 = stats.skew(tempS1.dropna())\n",
    "            col_skew_2 = stats.skew(tempS2.dropna())\n",
    "            col_skew_3 = stats.skew(tempS3.dropna())\n",
    "            col_skew_4 = stats.skew(tempS4.dropna())\n",
    "            col_skew_5 = stats.skew(tempS5.dropna())\n",
    "            col_skew_6 = stats.skew(tempS6.dropna())\n",
    "            col_skew_min = np.min([abs(col_skew_0), abs(col_skew_1), abs(col_skew_2), abs(col_skew_3),\n",
    "                                   abs(col_skew_4), abs(col_skew_5), abs(col_skew_6)])\n",
    "            if col_skew_min == abs(col_skew_1):\n",
    "                tempDF[col_n] = tempS1\n",
    "                counter_l1 += 1\n",
    "            elif col_skew_min == abs(col_skew_2):\n",
    "                tempDF[col_n] = tempS2\n",
    "                counter_l2 += 1\n",
    "            elif col_skew_min == abs(col_skew_3):\n",
    "                tempDF[col_n] = tempS3\n",
    "                counter_l3 += 1\n",
    "            elif col_skew_min == abs(col_skew_4):\n",
    "                #Flip horizontally again to maintain the direction of relationship\n",
    "                col_max = tempS4.max()\n",
    "                col_min = tempS4.min()\n",
    "                tempDF[col_n] = (col_max + col_min) - tempS4\n",
    "                counter_l4 += 1\n",
    "            elif col_skew_min == abs(col_skew_5):\n",
    "                #Flip horizontally again to maintain the direction of relationship\n",
    "                col_max = tempS5.max()\n",
    "                col_min = tempS5.min()\n",
    "                tempDF[col_n] = (col_max + col_min) - tempS5\n",
    "                counter_l5 += 1\n",
    "            elif col_skew_min == abs(col_skew_6):\n",
    "                #Flip horizontally again to maintain the direction of relationship\n",
    "                col_max = tempS6.max()\n",
    "                col_min = tempS6.min()\n",
    "                tempDF[col_n] = (col_max + col_min) - tempS6\n",
    "                counter_l6 += 1\n",
    "    print(' nAnalyte:', len(tempD[omics]))\n",
    "    print(' - log transformation for right skewed:', counter_r1)\n",
    "    print(' - cbrt transformation for right skewed:', counter_r2)\n",
    "    print(' - sqrt transformation for right skewed:', counter_r3)\n",
    "    print(' - exponential transformation for left skewed:', counter_l1)\n",
    "    print(' - cube transformation for left skewed:', counter_l2)\n",
    "    print(' - square transformation for left skewed:', counter_l3)\n",
    "    print(' - log transformation for left skewed with mirroring:', counter_l4)\n",
    "    print(' - cbrt transformation for left skewed with mirroring:', counter_l5)\n",
    "    print(' - sqrt transformation for left skewed with mirroring:', counter_l6)\n",
    "    print(' - no transformation:', len(tempD[omics]) - counter_r1 - counter_r2 - counter_r3 -\n",
    "          counter_l1 - counter_l2 - counter_l3 - counter_l4 - counter_l5 -counter_l6)\n",
    "    t_elapsed = time.time() - t_start\n",
    "    print(' Elapsed time:', round(t_elapsed//60), 'min', round(t_elapsed%60, 1), 'sec')\n",
    "    print('')\n",
    "combiDF_base = tempDF#Update the original"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Check\n",
    "print('After the step:')\n",
    "tempDF = combiDF_base.drop(columns=covarL)\n",
    "print('Confirm total number of NaNs in DF:', tempDF.isnull().to_numpy().sum(axis=None))\n",
    "#Skewness\n",
    "tempL = []\n",
    "for col_i in range(len(tempDF.columns)):\n",
    "    tempL.append(stats.skew(tempDF.iloc[:, col_i].dropna()))\n",
    "tempS = pd.Series(tempL, index=tempDF.columns, name='Skewness')\n",
    "display(tempS.describe())\n",
    "#Distribution and probability plot of examples\n",
    "sns.set(style='ticks', font='Arial', context='paper')\n",
    "fig, axes = plt.subplots(nrows=4, ncols=5, figsize=(15, 12), sharex=False, sharey=False)\n",
    "for ax_i, ax in enumerate(axes.flat):\n",
    "    if (ax_i//5)%2 == 0:\n",
    "        col_i = (ax_i//10)*5 + ax_i%5\n",
    "        sns.set(style='ticks', font='Arial', context='paper')\n",
    "        sns.distplot(tempDF.iloc[:, col_i].dropna(), color='g', ax=ax)\n",
    "        sns.despine()\n",
    "        ax.set_title(tempDF.columns[col_i], fontsize='large')\n",
    "        ax.set_xlabel('Raw or transformed value')\n",
    "    else:\n",
    "        sns.set(style='whitegrid', font='Arial', context='paper')\n",
    "        col_i = (ax_i//10)*5 + ax_i%5\n",
    "        stats.probplot(tempDF.iloc[:, col_i].dropna(), plot=ax)\n",
    "        ax.get_lines()[0].set_markerfacecolor('g')\n",
    "        ax.get_lines()[0].set_markeredgecolor('g')\n",
    "        ax.get_lines()[1].set_color('k')\n",
    "        ax.get_lines()[1].set_linewidth(3)\n",
    "        skewness = tempS.loc[tempDF.columns[col_i]]\n",
    "        ax.set_title('Skewness: '+str(round(skewness, 3)), fontsize='large')\n",
    "fig.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1-7. Prepare pairwise combinations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Prepare combination\n",
    "tempDF = combiDF_base.drop(columns=covarL)\n",
    "#Compute correlation matrix\n",
    "tempDF1 = tempDF.corr(method='pearson')\n",
    "print('nCombinations:', int(len(tempDF1)*(len(tempDF1)-1)/2))\n",
    "#Extract lower triangle matrix\n",
    "tempDF1 = tempDF1.where(np.tril(np.ones(tempDF1.shape), k=-1).astype(np.bool), other=np.nan)\n",
    "tempDF1.index.set_names('Variable1', inplace=True)\n",
    "tempDF1 = tempDF1.reset_index().melt(var_name='Variable2', value_name='Pearson_r', id_vars=['Variable1'])\n",
    "tempDF1 = tempDF1.dropna()\n",
    "print('nrows:', len(tempDF1))\n",
    "print('|Pearson\\'s r| > 0.8:', len(tempDF1.loc[abs(tempDF1['Pearson_r'])>0.8]))\n",
    "\n",
    "#Clean\n",
    "tempDF1['PairLabel'] = tempDF1['Variable1'].str.cat(tempDF1['Variable2'], sep=' vs. ')\n",
    "pairDF = tempDF1.set_index('PairLabel')\n",
    "display(pairDF)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1-8. Split DF for each model, and then drop NaN"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Sample size is different b/w analyte pairs.  \n",
    "> ***–> To maximize the sample size for each model, dropping NaN is performed after splitting.***  \n",
    "\n",
    "> ***The large number of combinations seems to halt this step due to lack of memory.***  \n",
    "> –> The subsequent preprocessing steps such as standardization should be performed at each modeling.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Baseline analyte–analyte interactions with MetBMI"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2-1. GLM: Analyte1 ~ b0 + b1\\*Analyte2 + b2\\*MetBMI + b3\\*Analyte2:MetBMI + b4\\*Sex + b5\\*Age + b6\\*AncestryPCs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Because skewness was reduced in preprocessing, use Gaussian family with an identity link function.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempL1 = []#For sample size\n",
    "tempL2 = []#For beta-coefficient\n",
    "tempL3 = []#For SE of beta-coefficient\n",
    "tempL4 = []#For t-statistic\n",
    "tempL5 = []#For residual degrees of freedom\n",
    "tempL6 = []#For P-value\n",
    "t_start = time.time()\n",
    "for pairs, var1, var2 in zip(pairDF.index, pairDF['Variable1'], pairDF['Variable2']):\n",
    "    #Split DF for each model, and then drop NaN\n",
    "    tempL = list(set(covarL) - set(['Race']))#Race can have NaN in this study and must be eliminated for dropna()\n",
    "    tempL = [variable for sublist in [[var1, var2], tempL] for variable in sublist]\n",
    "    tempDF = combiDF_base[tempL]\n",
    "    tempDF = tempDF.dropna()\n",
    "    \n",
    "    #Eliminate features with small sample size\n",
    "    ##-> Able to skip because analytes are already filtered with 10% missingness.\n",
    "    \n",
    "    #Standardization of continuous variables\n",
    "    tempDF1 = tempDF.select_dtypes(include=[np.number])\n",
    "    ##Z-score transformation\n",
    "    scaler = StandardScaler(copy=True, with_mean=True, with_std=True)\n",
    "    tempA = scaler.fit_transform(tempDF1)\n",
    "    tempDF1 = pd.DataFrame(data=tempA, index=tempDF1.index, columns=tempDF1.columns)\n",
    "    ##Recover the categorical covariates\n",
    "    tempDF2 = tempDF.select_dtypes(exclude=[np.number])\n",
    "    tempDF = pd.merge(tempDF1, tempDF2, left_index=True, right_index=True, how='left')\n",
    "    \n",
    "    #One-hot encoding for categorical covariates\n",
    "    ##-> In statsmodels, categorical variables are automatically recognized!\n",
    "    \n",
    "    #Generate GLM with Gaussian family\n",
    "    ##Rename the variables\n",
    "    tempDF = tempDF.rename(columns={var1:'Variable1', var2:'Variable2'})\n",
    "    ##Fit the model using Gaussian with an identity link function\n",
    "    formula = 'Variable1 ~ Variable2 * log_BaseMetBMI + C(Sex) + BaseAge + PC1 + PC2 + PC3 + PC4 + PC5'\n",
    "    model = smf.glm(formula=formula, data=tempDF,\n",
    "                    family=sm.families.Gaussian(sm.families.links.identity())).fit()\n",
    "    \n",
    "    #Save sample size\n",
    "    tempL1.append(len(tempDF))\n",
    "    #Save beta-coefficient of the interaction term\n",
    "    tempL2.append(model.params['Variable2:log_BaseMetBMI'])\n",
    "    tempL3.append(model.bse['Variable2:log_BaseMetBMI'])\n",
    "    #Save t-statistic of the interaction term\n",
    "    tempL4.append(model.tvalues['Variable2:log_BaseMetBMI'])\n",
    "    #Save residual degrees of freedom\n",
    "    tempL5.append(int(model.df_resid))\n",
    "    #Save P-value of the interaction term\n",
    "    tempL6.append(model.pvalues['Variable2:log_BaseMetBMI'])\n",
    "t_elapsed = time.time() - t_start\n",
    "print('Elapsed time for', len(pairDF), 'GLMs: ', round(t_elapsed//60), 'min', round(t_elapsed%60, 1), 'sec')\n",
    "\n",
    "#Clean the results\n",
    "pairDF['GLM_N'] = tempL1\n",
    "pairDF['GLM_Bcoef'] = tempL2\n",
    "pairDF['GLM_BcoefSE'] = tempL3\n",
    "pairDF['GLM_tStat'] = tempL4\n",
    "pairDF['GLM_DoF'] = tempL5\n",
    "pairDF['GLM_Pval'] = tempL6\n",
    "##P-value adjustment by using Benjamini–Hochberg method\n",
    "pairDF['GLM_AdjPval'] = multi.multipletests(pairDF['GLM_Pval'], alpha=0.05, method='fdr_bh',\n",
    "                                           is_sorted=False, returnsorted=False)[1]\n",
    "\n",
    "pairDF = pairDF.sort_values(by=['GLM_Pval'], ascending=True)\n",
    "display(pairDF)\n",
    "\n",
    "#Save\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220806_Multiomics-BMI-NatMed1stRevision_LongitudinalNetworkAnalysis-ver2_'\n",
    "fileName = 'GLM-interaction.tsv'\n",
    "pairDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2-2. Significant pairs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Import (just in case of time-out)\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220806_Multiomics-BMI-NatMed1stRevision_LongitudinalNetworkAnalysis-ver2_'\n",
    "fileName = 'GLM-interaction.tsv'\n",
    "pairDF = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t').set_index('PairLabel')\n",
    "pairDF = pairDF.loc[pairDF['GLM_AdjPval']<0.05]\n",
    "tempS = set(pairDF['Variable1']) | set(pairDF['Variable2'])\n",
    "print('Significant pairwise analyte interactions with MetBMI (FDR < 0.05):', len(pairDF))\n",
    "print('-> nAnalytes:', len(tempS))\n",
    "display(pairDF.iloc[0:30])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Example plots\n",
    "\n",
    "for row_i in range(5):\n",
    "    pairs = pairDF.index.tolist()[row_i]\n",
    "    print(pairs)\n",
    "    var1 = pairDF.iloc[row_i]['Variable1']\n",
    "    var2 = pairDF.iloc[row_i]['Variable2']\n",
    "    tempDF = combiDF_base[[var1, var2]]\n",
    "    tempDF = tempDF.dropna()\n",
    "    print('n =', len(tempDF))\n",
    "    #Z-score transformation\n",
    "    scaler = StandardScaler(copy=True, with_mean=True, with_std=True)\n",
    "    tempA = scaler.fit_transform(tempDF)\n",
    "    tempDF = pd.DataFrame(data=tempA, index=tempDF.index, columns=tempDF.columns)\n",
    "    #Add BMI class\n",
    "    tempDF = pd.merge(tempDF, combiDF_base['BaseMetBMI_class'],\n",
    "                      left_index=True, right_index=True, how='inner')\n",
    "    #Plot\n",
    "    tempD = {'Underweight':'blue', 'Normal':'green', 'Overweight':'orange', 'Obese':'red'}\n",
    "    sns.set(style='ticks', font='Arial', context='notebook')\n",
    "    p = sns.lmplot(data=tempDF, x=var1, y=var2,\n",
    "                   hue='BaseMetBMI_class', hue_order=tempD.keys(), palette=tempD,\n",
    "                   scatter_kws={'alpha':0.5, 'edgecolor':'0.3', 's':25}, height=3, aspect=1)\n",
    "    p.set(xlim=(-3.5, 3.5), xticks=np.arange(-2, 2.1, 2),\n",
    "          ylim=(-3.5, 3.5), yticks=np.arange(-2, 2.1, 2))\n",
    "    plt.xlabel(var1+'\\n('+r'$Z$'+'-score)')\n",
    "    plt.ylabel(var2+'\\n('+r'$Z$'+'-score)')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Preprocessing for GEE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3-1. Prepare dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Import the time-series all omics without imputation\n",
    "fileDir = '../210104_Biological-BMI-paper/ExportData/'\n",
    "ipynbName = '210104_Biological-BMI-paper_data-cleaning-BMI-omics_'\n",
    "fileName = 'time-series-combiDF-without-imputation_final-cohort.tsv'\n",
    "tempDF = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t', dtype={'public_client_id': str})\n",
    "tempDF = tempDF.set_index('KeyIndex')\n",
    "\n",
    "#Retrieve the baseline MetBMI-based class from the sex-stratified model\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220805_Multiomics-BMI-NatMed1stRevision_BMI-longitudinal-LASSO_'\n",
    "fileName = 'MetBMI-FemaleMale.tsv'\n",
    "tempDF1 = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t', dtype={'public_client_id': str})\n",
    "tempDF1 = tempDF1.sort_values(by=['public_client_id', 'days_in_program'], ascending=True)\n",
    "tempDF1 = tempDF1.drop_duplicates('public_client_id', keep='first')\n",
    "tempDF1 = tempDF1.set_index('public_client_id')\n",
    "tempDF1['log_BaseMetBMI'] = np.log(tempDF1['BaseMetBMI'])\n",
    "tempDF1 = tempDF1[['log_BaseMetBMI', 'BaseMetBMI_class']]\n",
    "tempDF = pd.merge(tempDF, tempDF1, left_on='public_client_id', right_index=True, how='left')\n",
    "\n",
    "#To apply the skewness reduction pipeline, revert the log-transformed metabolite values back to orignal scale\n",
    "covarL = ['log_BaseBMI', 'Sex', 'BaseAge', 'PC1', 'PC2', 'PC3', 'PC4', 'PC5', 'Race']\n",
    "fileDir = '../210104_Biological-BMI-paper/ExportData/'\n",
    "ipynbName = '210104_Biological-BMI-paper_data-cleaning-BMI-omics_'\n",
    "fileName = 'baseline-metDF-without-imputation_final-cohort.tsv'\n",
    "tempDF1 = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t', dtype={'public_client_id': str})\n",
    "tempDF1 = tempDF1.set_index('public_client_id')\n",
    "for col_n in tempDF1.drop(columns=covarL).columns.tolist():\n",
    "    tempDF[col_n] = np.e**tempDF[col_n]\n",
    "\n",
    "display(tempDF)\n",
    "#Update\n",
    "combiDF_ts = tempDF\n",
    "covarL = ['public_client_id', 'days_in_program', 'Season',\n",
    "          'log_BaseBMI', 'Sex', 'BaseAge', 'PC1', 'PC2', 'PC3', 'PC4', 'PC5', 'Race',\n",
    "          'log_BaseMetBMI', 'BaseMetBMI_class']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF = combiDF_ts#Not copy but just rename\n",
    "print('Time-series DF nrows before elimination:', len(tempDF))\n",
    "print(' -> Unique ID:', len(tempDF['public_client_id'].unique()))\n",
    "\n",
    "#Select measurements used in the longitudinal analysis\n",
    "month_threshold = 18\n",
    "tempDF = tempDF.loc[tempDF['days_in_program'] <= 365.25/12*month_threshold]\n",
    "print('Time-series DF nrows after selecting measurements during the taget period:', len(tempDF))\n",
    "print(' -> Unique ID:', len(tempDF['public_client_id'].unique()))\n",
    "\n",
    "#Select the cohort used in the longitudinal analysis\n",
    "##Select participants who have 2 or more measuremnts in BMI\n",
    "fileDir = '../210104_Biological-BMI-paper/ExportData/'\n",
    "ipynbName = '210104_Biological-BMI-paper_data-cleaning-BMI-omics_'\n",
    "fileName = 'time-series-bmiDF-without-imputation_final-cohort.tsv'\n",
    "tempDF1 = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t', dtype={'public_client_id': str})\n",
    "tempDF1 = tempDF1.set_index('KeyIndex')\n",
    "tempDF1 = tempDF1.loc[tempDF1['days_in_program'] <= 365.25/12*month_threshold]\n",
    "tempS = tempDF1['public_client_id'].value_counts()\n",
    "tempL1 = tempS.loc[tempS>=2].index.tolist()\n",
    "##Select participants who have 2 or more measuremnts in omics\n",
    "fileDir = '../210104_Biological-BMI-paper/ExportData/'\n",
    "ipynbName = '210104_Biological-BMI-paper_data-cleaning-BMI-omics_'\n",
    "fileName = 'time-series-combiDF-without-imputation_final-cohort.tsv'\n",
    "tempDF1 = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t', dtype={'public_client_id': str})\n",
    "tempDF1 = tempDF1.set_index('KeyIndex')\n",
    "tempDF1 = tempDF1.loc[tempDF1['days_in_program'] <= 365.25/12*month_threshold]\n",
    "tempS = tempDF1['public_client_id'].value_counts()\n",
    "tempL2 = tempS.loc[tempS>=2].index.tolist()\n",
    "##Select participants who have 2 or more measuremnts in both BMI and omics\n",
    "tempL = list(set(tempL1) & set(tempL2))\n",
    "print('Participants who have more than 2 measurements in both BMI and omics:', len(tempL))\n",
    "tempDF = tempDF.loc[tempDF['public_client_id'].isin(tempL)]\n",
    "print('Time-series DF nrows after selecting participants who have ≥2 measuremnts in both BMI and omics:',\n",
    "      len(tempDF))\n",
    "print(' -> Unique ID:', len(tempDF['public_client_id'].unique()))\n",
    "\n",
    "#Use only the measurements during period presented in the final figures of LMM\n",
    "tempDF = tempDF.loc[tempDF['days_in_program'] <= 365.25/12*12]\n",
    "print('Time-series DF nrows after selecting measurements during the final taget period:', len(tempDF))\n",
    "print(' -> Unique ID:', len(tempDF['public_client_id'].unique()))\n",
    "print(' -> nParticipants in metabolically obese group:',\n",
    "      len(tempDF.loc[tempDF['BaseMetBMI_class']=='Obese']['public_client_id'].unique()))\n",
    "\n",
    "display(tempDF)\n",
    "#Update\n",
    "combiDF_ts = tempDF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3-2. Drop pseudo numeric analytes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Some analytes in clinical labs have almost same values between participants.  \n",
    "> –> Because becoming invariant values after dropping outliers, these analytes will raise errors in GEE.  \n",
    "> –> Eliminate these in advance.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF = combiDF_ts#Not copy but just rename\n",
    "print('nAnalytes before the step: ', len(tempDF.drop(columns=covarL).columns))\n",
    "t_start = time.time()\n",
    "for col_n in tempDF.drop(columns=covarL).columns.tolist():\n",
    "    #Drop category level <15 (cf.) The category level of 'race' is 14.\n",
    "    if len(tempDF[col_n].unique())<15:\n",
    "        tempDF = tempDF.drop(columns=[col_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('nAnalytes after the step: ', len(tempDF.drop(columns=covarL).columns))\n",
    "combiDF_ts = tempDF#Update the original"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3-3. Replace zero/negative value with NaN"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Some analytes were measured as meaningful zeros in clinical labs. Also, analytes in proteomics could be negative values after the batch correction.  \n",
    "> –> Skip this step.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Check zero or negative values\n",
    "tempDF = combiDF_ts.drop(columns=covarL)\n",
    "print('Analytes having zero or negative values')\n",
    "tempS = tempDF.min(axis=0)\n",
    "display(tempDF.describe(include='all').loc[:, tempS.loc[tempS<=0].index])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3-4. Drop outliers (>3SD from MEAN)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Outliers are defined based on the baseline distribution.  \n",
    "> –> Because only using metabolically obese group in GEE, rather consider sample size to maximize as much as possible.  \n",
    "> –> By assuming that outliers occurred not participant-dependently but stochastically, use the whole distribution.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Check\n",
    "print('Before the step:')\n",
    "tempDF = combiDF_ts.drop(columns=covarL)\n",
    "display(tempDF.describe())\n",
    "#Distribution of examples\n",
    "sns.set(style='ticks', font='Arial', context='paper')\n",
    "fig, axes = plt.subplots(nrows=2, ncols=5, figsize=(15, 6), sharex=False, sharey=False)\n",
    "for col_i, ax in enumerate(axes.flat):\n",
    "    sns.distplot(tempDF.iloc[:, col_i].dropna(), color='g', ax=ax)\n",
    "    ax.set_title(tempDF.columns[col_i], fontsize='large')\n",
    "sns.despine()\n",
    "plt.setp(axes, xlabel='', ylabel='')\n",
    "plt.setp(axes[1, :], xlabel='Raw value')\n",
    "plt.setp(axes[:, 0], ylabel='Density')\n",
    "fig.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF = combiDF_ts#Not copy but just rename\n",
    "#Prepare the baseline DF\n",
    "#tempDF1 = tempDF.sort_values(by=['public_client_id', 'days_in_program'], ascending=True)\n",
    "#tempDF1 = tempDF1.drop_duplicates('public_client_id', keep='first')\n",
    "\n",
    "t_start = time.time()\n",
    "#Replace more than 3SD from MEAN with NaN based on the baseline distribution\n",
    "#-> Based on the whole distribution\n",
    "for col_n in tempDF.drop(columns=covarL).columns.tolist():\n",
    "    #col_mean = tempDF1[col_n].mean()\n",
    "    #col_sd = tempDF1[col_n].std(ddof=1)#Sample standard deviation: devided with N-1\n",
    "    col_mean = tempDF[col_n].mean()\n",
    "    col_sd = tempDF[col_n].std(ddof=1)#Sample standard deviation: devided with N-1\n",
    "    tempL = []\n",
    "    for row_n in tempDF.index.tolist():\n",
    "        value = tempDF.at[row_n, col_n]\n",
    "        if (value > (col_mean - 3*col_sd)) & (value < (col_mean + 3*col_sd)):\n",
    "            tempL.append(value)\n",
    "        else:#outlier\n",
    "            tempL.append(np.nan)\n",
    "    tempDF[col_n] = tempL\n",
    "t_elapsed = time.time() - t_start\n",
    "print('Elapsed time:', round(t_elapsed//60), 'min', round(t_elapsed%60, 1), 'sec')\n",
    "combiDF_ts = tempDF#Update the original"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Check\n",
    "print('After the step:')\n",
    "tempDF = combiDF_ts.drop(columns=covarL)\n",
    "display(tempDF.describe())\n",
    "#Distribution of examples\n",
    "sns.set(style='ticks', font='Arial', context='paper')\n",
    "fig, axes = plt.subplots(nrows=2, ncols=5, figsize=(15, 6), sharex=False, sharey=False)\n",
    "for col_i, ax in enumerate(axes.flat):\n",
    "    sns.distplot(tempDF.iloc[:, col_i].dropna(), color='g', ax=ax)\n",
    "    ax.set_title(tempDF.columns[col_i], fontsize='large')\n",
    "sns.despine()\n",
    "plt.setp(axes, xlabel='', ylabel='')\n",
    "plt.setp(axes[1, :], xlabel='Raw value')\n",
    "plt.setp(axes[:, 0], ylabel='Density')\n",
    "fig.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3-5. Drop analytes with small sample size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF = combiDF_ts#Not copy but just rename\n",
    "print('nAnalytes before the step: ', len(tempDF.drop(columns=covarL).columns))\n",
    "t_start = time.time()\n",
    "for col_n in tempDF.drop(columns=covarL).columns.tolist():\n",
    "    tempS = tempDF[col_n].dropna()\n",
    "    if len(tempS) < 100:\n",
    "        #Drop this analyte as one with small sample size\n",
    "        tempDF = tempDF.drop(columns=[col_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('nAnalytes after the step: ', len(tempDF.drop(columns=covarL).columns))\n",
    "combiDF_ts = tempDF#Update the original"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3-6. Transform to less skewed distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Check\n",
    "print('Before the step:')\n",
    "tempDF = combiDF_ts.drop(columns=covarL)\n",
    "print('Confirm total number of NaNs in DF:', tempDF.isnull().to_numpy().sum(axis=None))\n",
    "#Skewness\n",
    "tempL = []\n",
    "for col_i in range(len(tempDF.columns)):\n",
    "    tempL.append(stats.skew(tempDF.iloc[:, col_i].dropna()))\n",
    "tempS = pd.Series(tempL, index=tempDF.columns, name='Skewness')\n",
    "display(tempS.describe())\n",
    "#Distribution and probability plot of examples\n",
    "sns.set(style='ticks', font='Arial', context='paper')\n",
    "fig, axes = plt.subplots(nrows=4, ncols=5, figsize=(15, 12), sharex=False, sharey=False)\n",
    "for ax_i, ax in enumerate(axes.flat):\n",
    "    if (ax_i//5)%2 == 0:\n",
    "        col_i = (ax_i//10)*5 + ax_i%5\n",
    "        sns.set(style='ticks', font='Arial', context='paper')\n",
    "        sns.distplot(tempDF.iloc[:, col_i].dropna(), color='g', ax=ax)\n",
    "        sns.despine()\n",
    "        ax.set_title(tempDF.columns[col_i], fontsize='large')\n",
    "        ax.set_xlabel('Raw value')\n",
    "    else:\n",
    "        sns.set(style='whitegrid', font='Arial', context='paper')\n",
    "        col_i = (ax_i//10)*5 + ax_i%5\n",
    "        stats.probplot(tempDF.iloc[:, col_i].dropna(), plot=ax)\n",
    "        ax.get_lines()[0].set_markerfacecolor('g')\n",
    "        ax.get_lines()[0].set_markeredgecolor('g')\n",
    "        ax.get_lines()[1].set_color('k')\n",
    "        ax.get_lines()[1].set_linewidth(3)\n",
    "        skewness = tempS.loc[tempDF.columns[col_i]]\n",
    "        ax.set_title('Skewness: '+str(round(skewness, 3)), fontsize='large')\n",
    "fig.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tempDF = combiDF_ts#Not copy but just rename\n",
    "tempL = tempDF.drop(columns=covarL).columns.tolist()\n",
    "\n",
    "#Prepare analyte category list to check transformation per each omics\n",
    "fileDir = '../210104_Biological-BMI-paper/ExportData/'\n",
    "ipynbName = '210104_Biological-BMI-paper_data-cleaning-BMI-omics_'\n",
    "fileName = 'baseline-metDF-without-imputation_final-cohort.tsv'\n",
    "tempDF1 = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t', dtype={'public_client_id': str})\n",
    "tempDF1 = tempDF1.set_index('public_client_id')\n",
    "tempL1 = tempDF1.loc[:, ~tempDF1.columns.isin(covarL)].columns.tolist()\n",
    "tempL1 = list(set(tempL) & set(tempL1))\n",
    "fileName = 'baseline-protDF-without-imputation_final-cohort.tsv'\n",
    "tempDF1 = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t', dtype={'public_client_id': str})\n",
    "tempDF1 = tempDF1.set_index('public_client_id')\n",
    "tempL2 = tempDF1.loc[:, ~tempDF1.columns.isin(covarL)].columns.tolist()\n",
    "tempL2 = list(set(tempL) & set(tempL2))\n",
    "fileName = 'baseline-chemDF-without-imputation_final-cohort.tsv'\n",
    "tempDF1 = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t', dtype={'public_client_id': str})\n",
    "tempDF1 = tempDF1.set_index('public_client_id')\n",
    "tempL3 = tempDF1.loc[:, ~tempDF1.columns.isin(covarL)].columns.tolist()\n",
    "tempL3 = list(set(tempL) & set(tempL3))\n",
    "tempD = {'Metabolomics':tempL1, 'Proteomics':tempL2, 'Clinical labs':tempL3}\n",
    "\n",
    "for omics in tempD.keys():\n",
    "    t_start = time.time()\n",
    "    print(omics)\n",
    "    #Transform to less skewed distribution\n",
    "    counter_r1 = 0\n",
    "    counter_r2 = 0\n",
    "    counter_r3 = 0\n",
    "    counter_l1 = 0\n",
    "    counter_l2 = 0\n",
    "    counter_l3 = 0\n",
    "    counter_l4 = 0\n",
    "    counter_l5 = 0\n",
    "    counter_l6 = 0\n",
    "    for col_n in tempD[omics]:\n",
    "        tempS = tempDF[col_n].copy(deep=True)\n",
    "        col_skew_0 = stats.skew(tempS.dropna())\n",
    "        ##Shift the distribution to avoid NaN transformation for zero or negative values\n",
    "        tempS = tempS - tempS.min() + 1\n",
    "        if col_skew_0 > 0:#Right/positive skewed distribution\n",
    "            #Transformation\n",
    "            tempS1 = np.log(tempS)#Assumed for extremely skewed\n",
    "            tempS2 = np.cbrt(tempS)#Assumed for strongly skewed\n",
    "            tempS3 = np.sqrt(tempS)#Assumed for skewed\n",
    "            #Evaluation\n",
    "            col_skew_1 = stats.skew(tempS1.dropna())\n",
    "            col_skew_2 = stats.skew(tempS2.dropna())\n",
    "            col_skew_3 = stats.skew(tempS3.dropna())\n",
    "            col_skew_min = np.min([abs(col_skew_0), abs(col_skew_1), abs(col_skew_2), abs(col_skew_3)])\n",
    "            if col_skew_min == abs(col_skew_1):\n",
    "                tempDF[col_n] = tempS1\n",
    "                counter_r1 += 1\n",
    "            elif col_skew_min == abs(col_skew_2):\n",
    "                tempDF[col_n] = tempS2\n",
    "                counter_r2 += 1\n",
    "            elif col_skew_min == abs(col_skew_3):\n",
    "                tempDF[col_n] = tempS3\n",
    "                counter_r3 += 1\n",
    "        elif col_skew_0 < 0:#Left/negative skewed distribution\n",
    "            #Version1\n",
    "            ##Transformation\n",
    "            tempS1 = np.e**tempS#Assumed for extremely skewed\n",
    "            tempS2 = tempS**3#Assumed for strongly skewed\n",
    "            tempS3 = tempS**2#Assumed for skewed\n",
    "            #Version2\n",
    "            ##Flip horizontally to transform to the right skewed\n",
    "            col_max = tempS.max()\n",
    "            col_min = tempS.min()\n",
    "            tempS = (col_max + col_min) - tempS\n",
    "            ##Transformation\n",
    "            tempS4 = np.log(tempS)#Assumed for extremely skewed\n",
    "            tempS5 = np.cbrt(tempS)#Assumed for strongly skewed\n",
    "            tempS6 = np.sqrt(tempS)#Assumed for skewed\n",
    "            #Evaluation for both version 1 and 2\n",
    "            col_skew_1 = stats.skew(tempS1.dropna())\n",
    "            col_skew_2 = stats.skew(tempS2.dropna())\n",
    "            col_skew_3 = stats.skew(tempS3.dropna())\n",
    "            col_skew_4 = stats.skew(tempS4.dropna())\n",
    "            col_skew_5 = stats.skew(tempS5.dropna())\n",
    "            col_skew_6 = stats.skew(tempS6.dropna())\n",
    "            col_skew_min = np.min([abs(col_skew_0), abs(col_skew_1), abs(col_skew_2), abs(col_skew_3),\n",
    "                                   abs(col_skew_4), abs(col_skew_5), abs(col_skew_6)])\n",
    "            if col_skew_min == abs(col_skew_1):\n",
    "                tempDF[col_n] = tempS1\n",
    "                counter_l1 += 1\n",
    "            elif col_skew_min == abs(col_skew_2):\n",
    "                tempDF[col_n] = tempS2\n",
    "                counter_l2 += 1\n",
    "            elif col_skew_min == abs(col_skew_3):\n",
    "                tempDF[col_n] = tempS3\n",
    "                counter_l3 += 1\n",
    "            elif col_skew_min == abs(col_skew_4):\n",
    "                #Flip horizontally again to maintain the direction of relationship\n",
    "                col_max = tempS4.max()\n",
    "                col_min = tempS4.min()\n",
    "                tempDF[col_n] = (col_max + col_min) - tempS4\n",
    "                counter_l4 += 1\n",
    "            elif col_skew_min == abs(col_skew_5):\n",
    "                #Flip horizontally again to maintain the direction of relationship\n",
    "                col_max = tempS5.max()\n",
    "                col_min = tempS5.min()\n",
    "                tempDF[col_n] = (col_max + col_min) - tempS5\n",
    "                counter_l5 += 1\n",
    "            elif col_skew_min == abs(col_skew_6):\n",
    "                #Flip horizontally again to maintain the direction of relationship\n",
    "                col_max = tempS6.max()\n",
    "                col_min = tempS6.min()\n",
    "                tempDF[col_n] = (col_max + col_min) - tempS6\n",
    "                counter_l6 += 1\n",
    "    print(' nAnalyte:', len(tempD[omics]))\n",
    "    print(' - log transformation for right skewed:', counter_r1)\n",
    "    print(' - cbrt transformation for right skewed:', counter_r2)\n",
    "    print(' - sqrt transformation for right skewed:', counter_r3)\n",
    "    print(' - exponential transformation for left skewed:', counter_l1)\n",
    "    print(' - cube transformation for left skewed:', counter_l2)\n",
    "    print(' - square transformation for left skewed:', counter_l3)\n",
    "    print(' - log transformation for left skewed with mirroring:', counter_l4)\n",
    "    print(' - cbrt transformation for left skewed with mirroring:', counter_l5)\n",
    "    print(' - sqrt transformation for left skewed with mirroring:', counter_l6)\n",
    "    print(' - no transformation:', len(tempD[omics]) - counter_r1 - counter_r2 - counter_r3 -\n",
    "          counter_l1 - counter_l2 - counter_l3 - counter_l4 - counter_l5 -counter_l6)\n",
    "    t_elapsed = time.time() - t_start\n",
    "    print(' Elapsed time:', round(t_elapsed//60), 'min', round(t_elapsed%60, 1), 'sec')\n",
    "    print('')\n",
    "combiDF_ts = tempDF#Update the original"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Check\n",
    "print('After the step:')\n",
    "tempDF = combiDF_ts.drop(columns=covarL)\n",
    "print('Confirm total number of NaNs in DF:', tempDF.isnull().to_numpy().sum(axis=None))\n",
    "#Skewness\n",
    "tempL = []\n",
    "for col_i in range(len(tempDF.columns)):\n",
    "    tempL.append(stats.skew(tempDF.iloc[:, col_i].dropna()))\n",
    "tempS = pd.Series(tempL, index=tempDF.columns, name='Skewness')\n",
    "display(tempS.describe())\n",
    "#Distribution and probability plot of examples\n",
    "sns.set(style='ticks', font='Arial', context='paper')\n",
    "fig, axes = plt.subplots(nrows=4, ncols=5, figsize=(15, 12), sharex=False, sharey=False)\n",
    "for ax_i, ax in enumerate(axes.flat):\n",
    "    if (ax_i//5)%2 == 0:\n",
    "        col_i = (ax_i//10)*5 + ax_i%5\n",
    "        sns.set(style='ticks', font='Arial', context='paper')\n",
    "        sns.distplot(tempDF.iloc[:, col_i].dropna(), color='g', ax=ax)\n",
    "        sns.despine()\n",
    "        ax.set_title(tempDF.columns[col_i], fontsize='large')\n",
    "        ax.set_xlabel('Raw or transformed value')\n",
    "    else:\n",
    "        sns.set(style='whitegrid', font='Arial', context='paper')\n",
    "        col_i = (ax_i//10)*5 + ax_i%5\n",
    "        stats.probplot(tempDF.iloc[:, col_i].dropna(), plot=ax)\n",
    "        ax.get_lines()[0].set_markerfacecolor('g')\n",
    "        ax.get_lines()[0].set_markeredgecolor('g')\n",
    "        ax.get_lines()[1].set_color('k')\n",
    "        ax.get_lines()[1].set_linewidth(3)\n",
    "        skewness = tempS.loc[tempDF.columns[col_i]]\n",
    "        ax.set_title('Skewness: '+str(round(skewness, 3)), fontsize='large')\n",
    "fig.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3-7. Prepare the significant pairwise combinations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Import (just in case of time-out)\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220806_Multiomics-BMI-NatMed1stRevision_LongitudinalNetworkAnalysis-ver2_'\n",
    "fileName = 'GLM-interaction.tsv'\n",
    "pairDF = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t').set_index('PairLabel')\n",
    "pairDF = pairDF.loc[pairDF['GLM_AdjPval']<0.05]\n",
    "tempS = set(pairDF['Variable1']) | set(pairDF['Variable2'])\n",
    "print('Significant pairwise analyte interactions with MetBMI (FDR < 0.05):', len(pairDF))\n",
    "print('-> nAnalytes:', len(tempS))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3-8. Split DF for each model, and then drop NaN"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Sample size is different b/w analyte pairs.  \n",
    "> ***–> To maximize the sample size for each model, dropping NaN is performed after splitting.***  \n",
    "\n",
    "> ***The large number of combinations seems to halt this step due to lack of memory.***  \n",
    "> –> The subsequent preprocessing steps such as standardization should be performed at each modeling.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Analyte–analyte interactions with days in program in metabolically obese group"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4-1. GEE: Analyte1 ~ b0 + b1\\*Analyte2 + b2\\*Days + b3\\*Analyte2:Days + b4\\*Sex + b5\\*Age + b6\\*Season + b7\\*AncestryPCs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Because of multiple measurements per participant, use GEE with exchangeable covariance structure. As well as the above GLM, use Gaussian family with an identity link function.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "base_class = 'Obese'\n",
    "tempL1 = []#For sample size (observations)\n",
    "tempL2 = []#For the number of groups (= unique participants)\n",
    "tempL3 = []#For beta-coefficient\n",
    "tempL4 = []#For SE of beta-coefficient\n",
    "tempL5 = []#For t-statistic\n",
    "tempL6 = []#For P-value\n",
    "t_start = time.time()\n",
    "for pairs, var1, var2 in zip(pairDF.index, pairDF['Variable1'], pairDF['Variable2']):\n",
    "    #Split DF for each model, and then drop NaN\n",
    "    tempL = list(set(covarL) - set(['Race']))#Race can have NaN in this study and must be eliminated for dropna()\n",
    "    tempL = [variable for sublist in [[var1, var2], tempL] for variable in sublist]\n",
    "    tempDF = combiDF_ts[tempL]\n",
    "    tempDF = tempDF.dropna()\n",
    "    \n",
    "    #Select BMI class\n",
    "    tempDF = tempDF.loc[tempDF['BaseMetBMI_class']==base_class]\n",
    "    \n",
    "    #Standardization of continuous variables based on the whole distribution\n",
    "    tempDF1 = tempDF.select_dtypes(include=[np.number])\n",
    "    ##Z-score transformation\n",
    "    scaler = StandardScaler(copy=True, with_mean=True, with_std=True)\n",
    "    tempA = scaler.fit_transform(tempDF1)\n",
    "    tempDF1 = pd.DataFrame(data=tempA, index=tempDF1.index, columns=tempDF1.columns)\n",
    "    ##Recover the categorical covariates\n",
    "    tempDF2 = tempDF.select_dtypes(exclude=[np.number])\n",
    "    tempDF = pd.merge(tempDF1, tempDF2, left_index=True, right_index=True, how='left')\n",
    "    \n",
    "    #One-hot encoding for categorical covariates\n",
    "    ##-> In statsmodels, categorical variables are automatically recognized!\n",
    "    \n",
    "    #Generate GEE with exchangeable covariance structure\n",
    "    ##Rename the variables\n",
    "    tempDF = tempDF.rename(columns={var1:'Variable1', var2:'Variable2'})\n",
    "    ##Fit the model using exchangeable covariance structure\n",
    "    formula = 'Variable1 ~ Variable2 * days_in_program + C(Sex) + BaseAge + C(Season) + PC1 + PC2 + PC3 + PC4 + PC5'\n",
    "    model = sm.GEE.from_formula(formula=formula, data=tempDF, groups='public_client_id',\n",
    "                                family=sm.families.Gaussian(sm.families.links.identity()),\n",
    "                                cov_struct=sm.cov_struct.Exchangeable()).fit()\n",
    "    \n",
    "    #Save\n",
    "    #Save sample size\n",
    "    tempL1.append(len(tempDF))\n",
    "    #Save the number of groups\n",
    "    tempL2.append(len(tempDF['public_client_id'].unique()))\n",
    "    #Save beta-coefficient of the interaction term\n",
    "    tempL3.append(model.params['Variable2:days_in_program'])\n",
    "    tempL4.append(model.bse['Variable2:days_in_program'])\n",
    "    #Save t-statistic of the interaction term\n",
    "    tempL5.append(model.tvalues['Variable2:days_in_program'])\n",
    "    #Save P-value of the interaction term\n",
    "    tempL6.append(model.pvalues['Variable2:days_in_program'])\n",
    "t_elapsed = time.time() - t_start\n",
    "print('Elapsed time for', len(pairDF), 'GEEs: ', round(t_elapsed//60), 'min', round(t_elapsed%60, 1), 'sec')\n",
    "\n",
    "#Clean the results\n",
    "pairDF['GEE_nObs'] = tempL1\n",
    "pairDF['GEE_nGroups'] = tempL2\n",
    "pairDF['GEE_Bcoef'] = tempL3\n",
    "pairDF['GEE_BcoefSE'] = tempL4\n",
    "pairDF['GEE_tStat'] = tempL5\n",
    "pairDF['GEE_Pval'] = tempL6\n",
    "##P-value adjustment by using Benjamini–Hochberg method\n",
    "pairDF['GEE_AdjPval'] = multi.multipletests(pairDF['GEE_Pval'], alpha=0.05, method='fdr_bh',\n",
    "                                            is_sorted=False, returnsorted=False)[1]\n",
    "\n",
    "pairDF = pairDF.sort_values(by=['GEE_Pval'], ascending=True)\n",
    "display(pairDF)\n",
    "\n",
    "#Save\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220806_Multiomics-BMI-NatMed1stRevision_LongitudinalNetworkAnalysis-ver2_'\n",
    "fileName = 'GEE-interaction.tsv'\n",
    "pairDF.to_csv(fileDir+ipynbName+fileName, sep='\\t', index=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4-2. Significant pairs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Import (just in case of time-out)\n",
    "fileDir = './ExportData/'\n",
    "ipynbName = '220806_Multiomics-BMI-NatMed1stRevision_LongitudinalNetworkAnalysis-ver2_'\n",
    "fileName = 'GEE-interaction.tsv'\n",
    "pairDF = pd.read_csv(fileDir+ipynbName+fileName, sep='\\t').set_index('PairLabel')\n",
    "pairDF = pairDF.loc[pairDF['GEE_AdjPval']<0.05]\n",
    "tempS = set(pairDF['Variable1']) | set(pairDF['Variable2'])\n",
    "print('Significant pairwise analyte interactions with days in program (FDR < 0.05):', len(pairDF))\n",
    "print('-> nAnalytes:', len(tempS))\n",
    "display(pairDF)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> If Pearson's r \\* beta-coefficient in GLM > 0, the analyte correlation is positively modified by MetBMI.  \n",
    "> If Pearson's r \\* beta-coefficient in GEE > 0, the analyte correlation in obese group is positively modified by days.  \n",
    "> –> If beta-coefficient in GLM \\* beta-coefficient in GEE < 0, the association direction is concordant; that is, the analyte correlation in metabolically obese group is changed to that in normal group during the program.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Concordant pairs\n",
    "pairDF = pairDF.loc[pairDF['GLM_Bcoef']*pairDF['GEE_Bcoef'] < 0]\n",
    "tempS = set(pairDF['Variable1']) | set(pairDF['Variable2'])\n",
    "print('Concordant pairs:', len(pairDF))\n",
    "print('-> nAnalytes:', len(tempS))\n",
    "display(pairDF)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 4-2-1. Categorize with measurement number"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Example plots\n",
    "\n",
    "#Convert days_in_program value to categorical value\n",
    "##-> Although different meaning from GEE, try measurement number.\n",
    "##Ensure order\n",
    "combiDF_ts = combiDF_ts.sort_values(by=['public_client_id', 'days_in_program'], ascending=True)\n",
    "tempL = []\n",
    "counter = 1\n",
    "row_n_before = 'InitialDummyName'\n",
    "for row_n in combiDF_ts['public_client_id'].tolist():\n",
    "    if row_n == row_n_before:\n",
    "        counter += 1\n",
    "        tempL.append(counter)\n",
    "        row_n_before = row_n\n",
    "    else:\n",
    "        counter = 1\n",
    "        tempL.append(counter)\n",
    "        row_n_before = row_n\n",
    "combiDF_ts['MeasureNum'] = tempL\n",
    "\n",
    "for row_i in range(np.min([15, len(pairDF)])):\n",
    "    pairs = pairDF.index.tolist()[row_i]\n",
    "    print(pairs)\n",
    "    var1 = pairDF.iloc[row_i]['Variable1']\n",
    "    var2 = pairDF.iloc[row_i]['Variable2']\n",
    "    \n",
    "    #GLM result\n",
    "    tempDF = combiDF_base[[var1, var2]]\n",
    "    tempDF = tempDF.dropna()\n",
    "    print('GLM n =', len(tempDF))\n",
    "    ##Z-score transformation\n",
    "    scaler = StandardScaler(copy=True, with_mean=True, with_std=True)\n",
    "    tempA = scaler.fit_transform(tempDF)\n",
    "    tempDF = pd.DataFrame(data=tempA, index=tempDF.index, columns=tempDF.columns)\n",
    "    ##Add BMI class\n",
    "    tempDF = pd.merge(tempDF, combiDF_base['BaseMetBMI_class'],\n",
    "                      left_index=True, right_index=True, how='inner')\n",
    "    ##Plot\n",
    "    tempD = {'Underweight':'blue', 'Normal':'green', 'Overweight':'orange', 'Obese':'red'}\n",
    "    sns.set(style='ticks', font='Arial', context='notebook')\n",
    "    p = sns.lmplot(data=tempDF, x=var1, y=var2,\n",
    "                   hue='BaseMetBMI_class', hue_order=tempD.keys(), palette=tempD,\n",
    "                   scatter_kws={'alpha':0.5, 'edgecolor':'0.3', 's':25}, height=3, aspect=1)\n",
    "    p.set(xlim=(-3.5, 3.5), xticks=np.arange(-2, 2.1, 2),\n",
    "          ylim=(-3.5, 3.5), yticks=np.arange(-2, 2.1, 2))\n",
    "    plt.xlabel(var1+'\\n('+r'$Z$'+'-score)')\n",
    "    plt.ylabel(var2+'\\n('+r'$Z$'+'-score)')\n",
    "    plt.show()\n",
    "    \n",
    "    #GEE result\n",
    "    tempDF = combiDF_ts[[var1, var2, 'BaseMetBMI_class', 'MeasureNum']]\n",
    "    tempDF = tempDF.dropna()\n",
    "    ##Z-score transformation based on the baseline distribution\n",
    "    tempDF1 = tempDF.loc[tempDF['MeasureNum']==1].drop(columns=['BaseMetBMI_class', 'MeasureNum'])#baseline\n",
    "    tempDF2 = tempDF.drop(columns=['BaseMetBMI_class', 'MeasureNum'])#time-series\n",
    "    scaler = StandardScaler(copy=True, with_mean=True, with_std=True)\n",
    "    scaler.fit(tempDF1)\n",
    "    tempA = scaler.transform(tempDF2)\n",
    "    tempDF2 = pd.DataFrame(data=tempA, index=tempDF2.index, columns=tempDF2.columns)\n",
    "    ##Add category label\n",
    "    tempDF1 = tempDF[['BaseMetBMI_class', 'MeasureNum']]\n",
    "    tempDF = pd.merge(tempDF1, tempDF2, left_index=True, right_index=True, how='left')\n",
    "    tempD = {1:'1st', 2:'2nd', 3:'3rd and more', 4:'3rd and more', 5:'3rd and more'}\n",
    "    tempDF['Measurement'] = tempDF['MeasureNum'].map(tempD)\n",
    "    ##Select obese group\n",
    "    tempDF = tempDF.loc[tempDF['BaseMetBMI_class']=='Obese']\n",
    "    print('GEE n =', len(tempDF))\n",
    "    #Plot\n",
    "    tempD = {'1st':'red', '2nd':'orange', '3rd and more':'green'}\n",
    "    sns.set(style='ticks', font='Arial', context='notebook')\n",
    "    p = sns.lmplot(data=tempDF, x=var1, y=var2, hue='Measurement', hue_order=tempD.keys(), palette=tempD,\n",
    "                   scatter_kws={'alpha':0.5, 'edgecolor':'0.3', 's':25}, height=3, aspect=1)\n",
    "    p.set(xlim=(-3.5, 3.5), xticks=np.arange(-2, 2.1, 2),\n",
    "          ylim=(-3.5, 3.5), yticks=np.arange(-2, 2.1, 2))\n",
    "    plt.xlabel(var1+'\\n('+r'$Z$'+'-score)')\n",
    "    plt.ylabel(var2+'\\n('+r'$Z$'+'-score)')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 4-2-2. Categorize with binned days"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Example plots\n",
    "\n",
    "#Convert days_in_program value to categorical value\n",
    "tempL = []\n",
    "for value in combiDF_ts['days_in_program'].tolist():\n",
    "    if value < 365.25/12*4:\n",
    "        tempL.append('0–4 month')\n",
    "    elif value < 365.25/12*8:\n",
    "        tempL.append('4–8 month')\n",
    "    elif value <= 365.25/12*12:\n",
    "        tempL.append('8–12 month')\n",
    "    else:#Just in case\n",
    "        tempL.append('Error?')\n",
    "combiDF_ts['Period'] = tempL\n",
    "\n",
    "for row_i in range(np.min([15, len(pairDF)])):\n",
    "    pairs = pairDF.index.tolist()[row_i]\n",
    "    print(pairs)\n",
    "    var1 = pairDF.iloc[row_i]['Variable1']\n",
    "    var2 = pairDF.iloc[row_i]['Variable2']\n",
    "    \n",
    "    #GLM result\n",
    "    tempDF = combiDF_base[[var1, var2]]\n",
    "    tempDF = tempDF.dropna()\n",
    "    print('GLM n =', len(tempDF))\n",
    "    ##Z-score transformation\n",
    "    scaler = StandardScaler(copy=True, with_mean=True, with_std=True)\n",
    "    tempA = scaler.fit_transform(tempDF)\n",
    "    tempDF = pd.DataFrame(data=tempA, index=tempDF.index, columns=tempDF.columns)\n",
    "    ##Add BMI class\n",
    "    tempDF = pd.merge(tempDF, combiDF_base['BaseMetBMI_class'],\n",
    "                      left_index=True, right_index=True, how='inner')\n",
    "    ##Plot\n",
    "    tempD = {'Underweight':'blue', 'Normal':'green', 'Overweight':'orange', 'Obese':'red'}\n",
    "    sns.set(style='ticks', font='Arial', context='notebook')\n",
    "    p = sns.lmplot(data=tempDF, x=var1, y=var2,\n",
    "                   hue='BaseMetBMI_class', hue_order=tempD.keys(), palette=tempD,\n",
    "                   scatter_kws={'alpha':0.5, 'edgecolor':'0.3', 's':25}, height=3, aspect=1)\n",
    "    p.set(xlim=(-3.5, 3.5), xticks=np.arange(-2, 2.1, 2),\n",
    "          ylim=(-3.5, 3.5), yticks=np.arange(-2, 2.1, 2))\n",
    "    plt.xlabel(var1+'\\n('+r'$Z$'+'-score)')\n",
    "    plt.ylabel(var2+'\\n('+r'$Z$'+'-score)')\n",
    "    plt.show()\n",
    "    \n",
    "    #GEE result\n",
    "    tempDF = combiDF_ts[[var1, var2, 'BaseMetBMI_class', 'MeasureNum', 'Period']]\n",
    "    tempDF = tempDF.dropna()\n",
    "    ##Z-score transformation based on the baseline distribution\n",
    "    tempDF1 = tempDF.loc[tempDF['MeasureNum']==1].drop(columns=['BaseMetBMI_class', 'MeasureNum', 'Period'])#baseline\n",
    "    tempDF2 = tempDF.drop(columns=['BaseMetBMI_class', 'MeasureNum', 'Period'])#time-series\n",
    "    scaler = StandardScaler(copy=True, with_mean=True, with_std=True)\n",
    "    scaler.fit(tempDF1)\n",
    "    tempA = scaler.transform(tempDF2)\n",
    "    tempDF2 = pd.DataFrame(data=tempA, index=tempDF2.index, columns=tempDF2.columns)\n",
    "    ##Add category label\n",
    "    tempDF1 = tempDF[['BaseMetBMI_class', 'Period']]\n",
    "    tempDF = pd.merge(tempDF1, tempDF2, left_index=True, right_index=True, how='left')\n",
    "    ##Select obese group\n",
    "    tempDF = tempDF.loc[tempDF['BaseMetBMI_class']=='Obese']\n",
    "    print('GEE n =', len(tempDF))\n",
    "    #Plot\n",
    "    tempD = {'0–4 month':'red', '4–8 month':'orange', '8–12 month':'green'}\n",
    "    sns.set(style='ticks', font='Arial', context='notebook')\n",
    "    p = sns.lmplot(data=tempDF, x=var1, y=var2, hue='Period', hue_order=tempD.keys(), palette=tempD,\n",
    "                   scatter_kws={'alpha':0.5, 'edgecolor':'0.3', 's':25}, height=3, aspect=1)\n",
    "    p.set(xlim=(-3.5, 3.5), xticks=np.arange(-2, 2.1, 2),\n",
    "          ylim=(-3.5, 3.5), yticks=np.arange(-2, 2.1, 2))\n",
    "    plt.xlabel(var1+'\\n('+r'$Z$'+'-score)')\n",
    "    plt.ylabel(var2+'\\n('+r'$Z$'+'-score)')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Save representative plots for presentation\n",
    "tempL = ['2-oxoarginine* vs. alpha-hydroxyisocaproate',\n",
    "         'homoarginine vs. phenyllactate (PLA)',\n",
    "         '1-stearoyl-2-docosahexaenoyl-GPE (18:0/22:6)* vs. ALBUMIN']\n",
    "\n",
    "#GLM\n",
    "tempDF1 = pd.DataFrame()\n",
    "for pairs in tempL:\n",
    "    var1 = pairDF.at[pairs, 'Variable1']\n",
    "    var2 = pairDF.at[pairs, 'Variable2']\n",
    "    tempDF = combiDF_base[[var1, var2]]\n",
    "    tempDF = tempDF.dropna()\n",
    "    print(pairs, 'GLM: n =', len(tempDF))\n",
    "    #Z-score transformation\n",
    "    scaler = StandardScaler(copy=True, with_mean=True, with_std=True)\n",
    "    tempA = scaler.fit_transform(tempDF)\n",
    "    tempDF = pd.DataFrame(data=tempA, index=tempDF.index, columns=tempDF.columns)\n",
    "    #Add BMI class\n",
    "    tempDF = pd.merge(tempDF, combiDF_base['BaseMetBMI_class'],\n",
    "                      left_index=True, right_index=True, how='inner')\n",
    "    #Merge\n",
    "    tempDF['PairLabel'] = np.repeat(pairs, len(tempDF))\n",
    "    tempDF['Variable1'] = np.repeat(var1, len(tempDF))\n",
    "    tempDF['Variable2'] = np.repeat(var2, len(tempDF))\n",
    "    tempDF = tempDF.rename(columns={var1:'Zscore1', var2:'Zscore2'})\n",
    "    tempDF1 = pd.concat([tempDF1, tempDF], axis=0)\n",
    "#Plot\n",
    "tempD = {'Normal':'green', 'Overweight':'orange', 'Obese':'red'}\n",
    "tempD1 = {tempL[0]:'Intra-metabolomics', tempL[1]:'Intra-metabolomics', tempL[2]:'Inter-omics'}\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "fig, axes = plt.subplots(nrows=1, ncols=len(tempL),\n",
    "                         figsize=(3.75*len(tempL), 3.75+0.75), sharex=False, sharey=False)\n",
    "axis_xymin = -3.5\n",
    "axis_xymax = 3.5\n",
    "xymin = -2.0\n",
    "xymax = 2.0\n",
    "xyinter = 2.0\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",
    "    pairs = tempL[ax_i]\n",
    "    #Prepare DF\n",
    "    tempDF = tempDF1.loc[tempDF1['PairLabel']==pairs]\n",
    "    #Scatterplot with regression line\n",
    "    for bbmi_class in tempD.keys():\n",
    "        tempDF2 = tempDF.loc[tempDF['BaseMetBMI_class']==bbmi_class]\n",
    "        sns.regplot(data=tempDF2, x='Zscore1', y='Zscore2', color=tempD[bbmi_class],\n",
    "                    scatter=True, fit_reg=True, ci=95, truncate=False, marker='o',\n",
    "                    scatter_kws={'alpha':0.5, 'edgecolor':'0.3', 's':25}, ax=ax)\n",
    "    #Axis label\n",
    "    tempDF = tempDF.reset_index()[['PairLabel', 'Variable1', 'Variable2']]\n",
    "    tempDF = tempDF.drop_duplicates('PairLabel', keep='first').set_index('PairLabel')\n",
    "    var1 = tempDF.at[pairs, 'Variable1']\n",
    "    var2 = tempDF.at[pairs, 'Variable2']\n",
    "    ax.set_xlabel(var1.replace('-GPE', '-\\nGPE')+'\\n('+r'$Z$'+'-score)')\n",
    "    ax.set_ylabel(var2.replace('-GPE', '-\\nGPE')+'\\n('+r'$Z$'+'-score)')\n",
    "    #Facet title\n",
    "    ax.set_title(tempD1[pairs], {'fontsize':'large'})\n",
    "sns.despine()\n",
    "plt.tight_layout()\n",
    "#Generate legend; default point in lmplot is tiny\n",
    "tempL1 = []\n",
    "for bmi_class in tempD.keys():\n",
    "    tempL1.append(mlines.Line2D([], [], color=tempD[bmi_class], label=bmi_class, linewidth=4))\n",
    "plt.legend(handles=tempL1, fontsize='medium', title='Baseline MetBMI class', title_fontsize='large',\n",
    "           bbox_to_anchor=(1, 0.5), loc='center left', borderaxespad=3, frameon=True)\n",
    "##Save\n",
    "fileDir = './ExportFigures/'\n",
    "ipynbName = '220806_Multiomics-BMI-NatMed1stRevision_LongitudinalNetworkAnalysis-ver2_'\n",
    "fileName = 'GLM-interaction-examples.tif'\n",
    "plt.gcf().savefig(fileDir+ipynbName+fileName, dpi=300, bbox_inches='tight', pad_inches=0.04,\n",
    "                  pil_kwargs={'compression':'tiff_lzw'})\n",
    "plt.show()\n",
    "\n",
    "#GEE result\n",
    "tempDF1 = pd.DataFrame()\n",
    "for pairs in tempL:\n",
    "    var1 = pairDF.at[pairs, 'Variable1']\n",
    "    var2 = pairDF.at[pairs, 'Variable2']\n",
    "    tempDF = combiDF_ts[[var1, var2, 'BaseMetBMI_class', 'MeasureNum', 'Period']]\n",
    "    tempDF = tempDF.dropna()\n",
    "    #Z-score transformation based on the baseline distribution\n",
    "    tempDF2 = tempDF.loc[tempDF['MeasureNum']==1].drop(columns=['BaseMetBMI_class', 'MeasureNum', 'Period'])#baseline\n",
    "    tempDF3 = tempDF.drop(columns=['BaseMetBMI_class', 'MeasureNum', 'Period'])#time-series\n",
    "    scaler = StandardScaler(copy=True, with_mean=True, with_std=True)\n",
    "    scaler.fit(tempDF2)\n",
    "    tempA = scaler.transform(tempDF3)\n",
    "    tempDF3 = pd.DataFrame(data=tempA, index=tempDF3.index, columns=tempDF3.columns)\n",
    "    #Add category label\n",
    "    tempDF2 = tempDF[['BaseMetBMI_class', 'Period']]\n",
    "    tempDF = pd.merge(tempDF2, tempDF3, left_index=True, right_index=True, how='left')\n",
    "    #Select obese group\n",
    "    tempDF = tempDF.loc[tempDF['BaseMetBMI_class']=='Obese']\n",
    "    print(pairs, 'GEE: n =', len(tempDF))\n",
    "    #Merge\n",
    "    tempDF['PairLabel'] = np.repeat(pairs, len(tempDF))\n",
    "    tempDF['Variable1'] = np.repeat(var1, len(tempDF))\n",
    "    tempDF['Variable2'] = np.repeat(var2, len(tempDF))\n",
    "    tempDF = tempDF.rename(columns={var1:'Zscore1', var2:'Zscore2'})\n",
    "    tempDF1 = pd.concat([tempDF1, tempDF], axis=0)\n",
    "#Plot\n",
    "tempD = {'0–4 month':'red', '4–8 month':'orange', '8–12 month':'green'}\n",
    "tempD1 = {tempL[0]:'Intra-metabolomics', tempL[1]:'Intra-metabolomics', tempL[2]:'Inter-omics'}\n",
    "sns.set(style='ticks', font='Arial', context='talk')\n",
    "fig, axes = plt.subplots(nrows=1, ncols=len(tempL),\n",
    "                         figsize=(3.75*len(tempL), 3.75+0.75), sharex=False, sharey=False)\n",
    "axis_xymin = -3.5\n",
    "axis_xymax = 3.5\n",
    "xymin = -2.0\n",
    "xymax = 2.0\n",
    "xyinter = 2.0\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",
    "    pairs = tempL[ax_i]\n",
    "    #Prepare DF\n",
    "    tempDF = tempDF1.loc[tempDF1['PairLabel']==pairs]\n",
    "    #Scatterplot with regression line\n",
    "    for period in tempD.keys():\n",
    "        tempDF2 = tempDF.loc[tempDF['Period']==period]\n",
    "        sns.regplot(data=tempDF2, x='Zscore1', y='Zscore2', color=tempD[period],\n",
    "                    scatter=True, fit_reg=True, ci=95, truncate=False, marker='o',\n",
    "                    scatter_kws={'alpha':0.5, 'edgecolor':'0.3', 's':25}, ax=ax)\n",
    "    #Axis label\n",
    "    tempDF = tempDF.reset_index()[['PairLabel', 'Variable1', 'Variable2']]\n",
    "    tempDF = tempDF.drop_duplicates('PairLabel', keep='first').set_index('PairLabel')\n",
    "    var1 = tempDF.at[pairs, 'Variable1']\n",
    "    var2 = tempDF.at[pairs, 'Variable2']\n",
    "    ax.set_xlabel(var1.replace('-GPE', '-\\nGPE')+'\\n('+r'$Z$'+'-score)')\n",
    "    ax.set_ylabel(var2.replace('-GPE', '-\\nGPE')+'\\n('+r'$Z$'+'-score)')\n",
    "    #Facet title\n",
    "    ax.set_title(tempD1[pairs], {'fontsize':'large'})\n",
    "sns.despine()\n",
    "plt.tight_layout()\n",
    "#Generate legend; default point in lmplot is tiny\n",
    "tempL1 = []\n",
    "for bmi_class in tempD.keys():\n",
    "    tempL1.append(mlines.Line2D([], [], color=tempD[bmi_class], label=bmi_class, linewidth=4))\n",
    "plt.legend(handles=tempL1, fontsize='medium', title='Period in program', title_fontsize='large',\n",
    "           bbox_to_anchor=(1, 0.5), loc='center left', borderaxespad=3, frameon=True)\n",
    "##Save\n",
    "fileDir = './ExportFigures/'\n",
    "ipynbName = '220806_Multiomics-BMI-NatMed1stRevision_LongitudinalNetworkAnalysis-ver2_'\n",
    "fileName = 'GEE-interaction-examples.tif'\n",
    "plt.gcf().savefig(fileDir+ipynbName+fileName, dpi=300, bbox_inches='tight', pad_inches=0.04,\n",
    "                  pil_kwargs={'compression':'tiff_lzw'})\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# — End of notebook —"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "arivale-py3 - Python",
   "language": "python",
   "name": "conda-env-arivale-py3-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  },
  "toc-showcode": false,
  "toc-showmarkdowntxt": false,
  "toc-showtags": false,
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}