Univariate Analysis-NMR.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "OkiuAzj7SMMq"
      },
      "source": [
        "# Univariate Analysis of NMR Data"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ljaFjrKzSMMt"
      },
      "source": [
        "In this notebook we will perform a *univariate* analysis of 4 racks from the NMR profiling dataset of the *Dementia research cohort*, as an alternative to the mulivariate PLS-DA analysis. For details of the study see the Metabolights Study [MTBLS719](https://www.ebi.ac.uk/metabolights/MTBLS719).\n",
        "\n",
        "This analysis is comprised of the following steps:\n",
        "\n",
        "1) For each signal in a given chemical shift position (ppm) we fit a linear regression model explaining variance in signal intensity based on a set of fixed covariates (age and gender). This model will be used to test for significant differences in signal intensity due between gender.\n",
        "\n",
        "2) After obtaining the *p-values* and regression coefficient estimates from step 1), multiple testing correction is applied to control the number of false positives.\n",
        "\n",
        "3) As an extra step to assess robustness of findings, we also perform bootstrap resampling of the models fitted in 1).\n",
        "\n",
        "4) Visualization of the regions of the spectrum associated with each covariate of interest.\n",
        "\n",
        "In this examples we use linear regression to fit a model and then perform statistical hypothesis tests, but these analyses can be performed also by using *t-tests*, one-way and two-way ANOVA models directly, or even their non-parametric equivalent."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dkabRMsvSMMu"
      },
      "source": [
        "## Package import and environment setup"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "TAGbPtB_SMMu"
      },
      "source": [
        "Import all the packages and configure notebook plotting mode."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "y_5T7r2TSMMu"
      },
      "outputs": [],
      "source": [
        "!pip install -q ipympl\n",
        "!pip install -q kneed"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "zGFO-piFSMMv"
      },
      "outputs": [],
      "source": [
        "!git clone https://github.com/IPTC-DataAnalysisCourse/chemometrics-tutorials.git\n",
        "%cd chemometrics-tutorials"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "QSJdL2mvSMMv"
      },
      "outputs": [],
      "source": [
        "# Import the required python packages including\n",
        "# the custom Chemometric Model objects\n",
        "import numpy as np\n",
        "from sklearn import preprocessing\n",
        "import pandas as pds\n",
        "import matplotlib.pyplot as plt\n",
        "import matplotlib as mpl\n",
        "from matplotlib import colormaps\n",
        "import statsmodels.formula.api as smf\n",
        "import numpy as np\n",
        "\n",
        "from plotly.offline import init_notebook_mode, iplot\n",
        "init_notebook_mode(connected=True)\n",
        "\n",
        "from statsmodels.stats.multitest import multipletests\n",
        "\n",
        "# The custom Chemometric Model objects for plotting\n",
        "from pyChemometrics.plotting_utils import manhattan_plot, interactive_manhattan\n",
        "from pyChemometrics.plotting_utils import _scatterplots as scatterplot\n",
        "from pyChemometrics.plotting_utils import _lineplots as Lineplot"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "hz4dU2VWSMMw"
      },
      "source": [
        "The next cell sets up the figure display mode. The *notebook* and *ipympl* modes allows interactive plotting. Another option is to select *inline*, to obtain static plots in a notebook cell."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "CnOt8lc6SMMw"
      },
      "outputs": [],
      "source": [
        "# Set the plot backend to support interactive plotting\n",
        "#%matplotlib notebook # run this line when running in jupyter notebook"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "ZZiIGYsCSMMw"
      },
      "outputs": [],
      "source": [
        "# Set the plot backend to support interactive plotting - Run this when running in Google Colab\n",
        "%matplotlib ipympl\n",
        "\n",
        "import plotly.io as pio\n",
        "pio.renderers.default = \"colab\"\n",
        "\n",
        "from google.colab import output\n",
        "output.enable_custom_widget_manager()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "vvoPw3M8SMMx"
      },
      "source": [
        "## Data Import"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zp53O9KGSMMx"
      },
      "source": [
        "We will now import the NMR data and the metadata (Y variables).\n",
        "\n",
        "X - NMR data matrix\n",
        "\n",
        "Y - Matrix with the 2 metadata outcomes\n",
        "\n",
        "ppm - Chemical shift axis for the NMR data in H $\\delta$ppm.\n",
        "\n",
        "#### Metadata\n",
        "Y1 - represents the sex (1: Male, 0: Female, in original Y data matrix)\n",
        "\n",
        "Y2 - represents the age (varies from 61 - 94, in original Y data matrix)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "DXkX2_iASMMx"
      },
      "outputs": [],
      "source": [
        "# Load the dataset\n",
        "dementia_nmr_dataset = pds.read_csv(\"./Data/Dementia U NMR_PQN_normalised_Data.csv\",delimiter=',')\n"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# Delete samples where outcome variable is unknown - Study Samples in standard NPC pipeline\n",
        "dementia_nmr_dataset = dementia_nmr_dataset[~dementia_nmr_dataset['Gender'].isnull()]\n",
        "\n",
        "# Create the X matrix\n",
        "X = dementia_nmr_dataset.iloc[:, 5::].values\n",
        "\n",
        "# Use pandas Categorical type to generate the dummy enconding of the Y1 vector (0 and 1)\n",
        "Y1 = pds.Categorical(dementia_nmr_dataset['Gender']).codes\n",
        "\n",
        "#Use any regular numpy array for continuous variables\n",
        "Y2 = np.array(dementia_nmr_dataset['Age'], dtype='int')\n",
        "\n",
        "ppm = np.array(dementia_nmr_dataset.columns[5::], dtype=np.float64)"
      ],
      "metadata": {
        "id": "rJ8ubvWhD7CK"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ukDA5q-RSMMx"
      },
      "source": [
        "**Note**: To apply the analyses exemplified in this notebook to any other dataset, just modify the cell above to import the data matrices and vectors X and Y from any other source file.\n",
        "\n",
        "The expected data types and formatting for **X** and **Y** are:\n",
        "\n",
        "   **X**: Any data matrix with n rows (observations/samples) and p columns (variables/features). The matrix should be provided as a [numpy.ndarray](https://numpy.org/doc/stable/reference/generated/numpy.ndarray.html) object, with 2 dimensions, and with shape = (n, p). We recommend using the *numpy* function [numpy.genfromtxt](https://numpy.org/devdocs/reference/generated/numpy.genfromtxt.html) or the *pandas* [pandas.read_csv](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.read_csv.html) function to read the data from a text file. When using the *pandas.read_csv* function, extract the data matrix as a *numpy.ndarray* from the pandas.DataFrame object using the `.values` attribute.\n",
        "```\n",
        "X_DataFrame = pds.read_csv(\"./data/X_spectra.csv\")\n",
        "X = X_DataFrame.values\n",
        "```\n",
        "   \n",
        "   **Y** vectors: Each **Y** vector should be a 1-dimensional [numpy.ndarray](https://numpy.org/doc/stable/reference/generated/numpy.ndarray.html) object, with a number and ordering of elements matching the rows in **X**. For continuous variables, any regular *numpy.ndarray* with a data type of `int` (integers only) or `float` can be used.\n",
        "   ```\n",
        "   Y_continuous = numpy.ndarray([23.4, 24, 0.3, -1.23], dtype='float')\n",
        "   ```\n",
        "To encode binary class labels, a *numpy.ndarray* of dtype `int`, with 0 and 1 as labels (e.g., 0 = Control, 1 = Case) must be used. The way in which classes are encoded will affect the model interpretation: the class labeled as 1 is used as the \"positive/case\" class by the *pyChemometrics* objects.\n",
        "   \n",
        "   In the example above, we used the *pandas* [Categorical](https://pandas.pydata.org/pandas-docs/stable/user_guide/categorical.html) datatype to handle the conversion of the original categorical variable (Male, Female) to the required (0, 1) labels. After converting a column to a `Categorical` datatype, the `.codes` attribute returns a vector with the same length of the original Y, but where each value is replaced by their integer (`int`) code. The correspondence between code and category can be inspected with the `categories` attribute. The order of the labels in `.codes` is the same as the order of the `categories` attribute (i.e. 0 is the first element in `categories`, 1 the second and so on).\n",
        "   ```\n",
        "   Y1 = pds.Categorical(Y.iloc[:, 1])\n",
        "   Y1.codes # The numerical label\n",
        "   Y1.categories # Original text or numerical description of the category\n",
        "   ```\n",
        "   [get_dummies](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.get_dummies.html) is another helpful function to perform dummy (0-1) encoding of variables.\n",
        "   \n",
        "For the linear regression models, we will also generate a pandas DataFrame contaning the metabolic and class variables. This is done to take advantage of the more explicit `formula` synthax from the *statsmodels* package."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "S-Z7Ld3nSMMx"
      },
      "outputs": [],
      "source": [
        "# Prepare the dataframe for use with model fitting\n",
        "dataset = pds.DataFrame(data = X, columns = ppm)\n",
        "\n",
        "# Add the Y vector at the beginning of the dataframe\n",
        "dataset.insert(0, 'Gender', Y1)\n",
        "dataset.insert(1, 'Age', Y2)\n",
        "\n",
        "# Display dataframe\n",
        "dataset"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0KQjC-HMSMMx"
      },
      "source": [
        "Plot the spectra in the dataset."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "ChDW0RnDSMMx"
      },
      "outputs": [],
      "source": [
        "# Plot the dataset\n",
        "plt.figure()\n",
        "plt.plot(ppm, X.T)\n",
        "plt.title(\"X matrix of spectra\")\n",
        "plt.xlabel(\"$\\delta$ppm\")\n",
        "plt.gca().invert_xaxis()\n",
        "plt.ylabel(\"Intensity\")\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "FMq9eNGXSMMy"
      },
      "source": [
        "## Linear regression modeling - effect of gender\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "IZ7GXDGgSMMy"
      },
      "source": [
        "We start by assessing the effect of gender on the intensity of each individual spectral signal using linear regression.\n",
        "\n",
        "This model is equivalent to a 2-sample (each of the gender) *t-test* on signal intensity or a 1-way ANOVA with Gender as a single discrete covariate/factor. Other modelling alternatives are possible, for example, logistic regression.\n",
        "\n",
        "The results from this model should be comparable with a PLS-DA model regressing the gender (Y) on the metabolic profile variables (X). Note how the $R^{2}Y$ in the PLS-DA model refers to how well the metabolic profile is able to predict gender, while for the linear models $r^{2}$ measures how much of \"metabolite\" variability is explained by gender.\n",
        "\n",
        "In next cell we will iteratively fit the same regression model, Metabolite(y) ~ Gender for each metabolic profile variable (a given point in the chemical shift *x-axis* in $\\delta$ppm) value.\n",
        "\n",
        "For each model, we obtain the calculated $r^{2}$ (variance explained),\n",
        "regression coefficients ($\\beta$) for each covariate (Gender only for this first model) and the *p-value* of a two-tailed Wald test ($\\beta$ significantly different from 0)."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "TcyFOt9kSMMy"
      },
      "outputs": [],
      "source": [
        "# Lists to store the information\n",
        "# p-value for the gender effect\n",
        "pval_genonly = list()\n",
        "# regression coefficient for the gender effect\n",
        "beta_genonly = list()\n",
        "# P-value for the F-test\n",
        "fpval_genonly = list()\n",
        "# r2 for the regression model\n",
        "r2_genonly = list()\n",
        "\n",
        "# Fit each column with a spectral variable\n",
        "for curr_variable in dataset.iloc[:, 2:]:\n",
        "    # Formula for current variable\n",
        "    fm = 'Q(curr_variable) ~ C(Gender)'\n",
        "    mod = smf.ols(formula = fm, data=dataset)\n",
        "    res = mod.fit()\n",
        "    pval_genonly.append(res.pvalues[1])\n",
        "    beta_genonly.append(res.params[1])\n",
        "    fpval_genonly.append(res.f_pvalue)\n",
        "    r2_genonly.append(res.rsquared)\n",
        "\n",
        "pval_genonly = np.array(pval_genonly)\n",
        "beta_genonly = np.array(beta_genonly)\n",
        "r2_genonly = np.array(r2_genonly)\n",
        "fpval_genonly = np.array(fpval_genonly)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zwrAMiksSMMy"
      },
      "source": [
        "## Linear modeling per variable - with adjustment for age"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "SLWYydXsSMMy"
      },
      "source": [
        "Adjustment for other covariates and potential confounders can be easily incorporated in linear regression.\n",
        "In this case, gender and age will be modelled together. The confounding variation caused by Age will be separated from the effect of gender.\n",
        "\n",
        "In next cell we will iteratively fit the same regression model, Metabolite(y) ~ Gender + Age for each metabolic profile variable (ppm) value.\n",
        "\n",
        "For each model, we obtain the calculated $r^{2}$ (variance explained),\n",
        "regression coefficients ($\\beta$) for each covariate (Gender and Age) and the *p-value* of a two-tailed Wald test ($\\beta$ significantly different from 0)."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "SspL8LqHSMMy"
      },
      "outputs": [],
      "source": [
        "# Generate lists to store the results for each varable.\n",
        "\n",
        "\n",
        "# beta_gen = the values for the regression coefficient associated with age\n",
        "# pval_gen = the p-value for the t-test\n",
        "# beta_age = the values for the age regression coefficients\n",
        "# pval_age = the p-values for the age regression coefficient\n",
        "\n",
        "# r2 = r-squared value for the entire regression model\n",
        "# fpval = the p-value for the F-test of the model's r-squared\n",
        "\n",
        "pval_gen = list()\n",
        "pval_age = list()\n",
        "beta_gen = list()\n",
        "beta_age = list()\n",
        "fpval = list()\n",
        "r2 = list()\n",
        "\n",
        "for curr_variable in dataset.iloc[:, 2:]:\n",
        "    # Fit each column\n",
        "    fm = 'Q(curr_variable) ~ C(Gender) + C(Age)'\n",
        "    mod = smf.ols(formula = fm, data=dataset)\n",
        "    res = mod.fit()\n",
        "    pval_gen.append(res.pvalues[1])\n",
        "    pval_age.append(res.pvalues[2])\n",
        "    beta_gen.append(res.params[1])\n",
        "    beta_age.append(res.params[2])\n",
        "    fpval.append(res.f_pvalue)\n",
        "    r2.append(res.rsquared)\n",
        "\n",
        "pval_gen = np.array(pval_gen)\n",
        "beta_gen = np.array(beta_gen)\n",
        "pval_age = np.array(pval_age)\n",
        "beta_age = np.array(beta_age)\n",
        "r2 = np.array(r2)\n",
        "fpval = np.array(fpval)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "V94ULhPGSMMy"
      },
      "outputs": [],
      "source": [
        "sum(pval_gen < 0.05)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6wA6XohxSMMy"
      },
      "source": [
        "## Multiple testing correction\n",
        "\n",
        "One of the problems associated with carrying out a large amount of sequential univariate tests is the potentially high number of false positives. In this instance, we have performed 17456 tests, and at an $\\alpha$ threshold of 0.05 the expected number of false positives due to chance alone is equal to $17456~\\textrm{(number of tests)} \\times 0.05m(\\alpha) = 872.8$ ~ $873$\n",
        "\n",
        "The set of *p-values* obtained by this multiple testing strategy should be corrected to restrict the expected large number of false positives.\n",
        "\n",
        "In this example we use the Benjamini-Yekutieli procedure to control the false discovery\n",
        "rate (FDR) of the sequential univariate testing procedure, the proportion of \"false discoveries\" (false positives) in the metabolic signature uncovered (the set of positive results/rejected null hypotheses from the 18910 univariate tests).\n",
        "\n",
        "Alternatively, Family-Wise error rate (FWER) correction methods, such as the Bonferroni correction can be used. These re-calibrate the significance threshold so that the FWER, the probability of having one or more false positives in the entirety of the analysis, is kept at a nominal $\\alpha$ level (typically 0.05). FWER correction methods tend to be more stringent than FDR correction methods, but the concept of False discovery rate control is particularly suited for metabolic profiling analyses, because FDR control methods aims at controlling the general \"quality\" (proportion of false positives) of the uncovered signatures.\n",
        "\n",
        "**Note**: When using FWER adjustment, we use the original *p-values* and simply adjust the $\\alpha$ significance level to meet the FWER requirements. The Benjamini-Hochberg and Benjamini-Yekutieli adjustment procedures return new quantities, a set of *q-values* (one per original *p-value*). This new *q-value* should be interpreted in the context of the entire signature, and it specifies the estimated false discovery rate if all features with a *q-value* smaller or equal than the cut-off are selected as hits. For example, setting a 0.05 cut-off implies that the obtained signature has an estimated false discovery rate of 5%."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "mVtQvAEESMMy"
      },
      "outputs": [],
      "source": [
        "# Adjusting the first analysis without age\n",
        "by_res_gen_only = multipletests(pval_genonly, alpha=0.05, method='fdr_by')\n",
        "p_genonly_byadj = by_res_gen_only[1]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "a4kkOT3RSMMy"
      },
      "outputs": [],
      "source": [
        "# Adjusting the analysis with gender and age\n",
        "by_res_gen = multipletests(pval_gen, alpha=0.05, method='fdr_by')\n",
        "p_gen_byadj = by_res_gen[1]\n",
        "\n",
        "by_res_f = multipletests(fpval, alpha=0.05, method='fdr_by')\n",
        "p_byadj_f = by_res_f[1]\n",
        "\n",
        "by_res_age = multipletests(pval_age, alpha=0.05, method='fdr_by')\n",
        "p_age_byadj = by_res_age[1]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "nUncf7HdSMMy"
      },
      "outputs": [],
      "source": [
        "#Number of significant features\n",
        "print(\"Number of features significantly associated with gender\")\n",
        "\n",
        "print(\"Without adjustment for age: {0}\".format(sum(by_res_gen_only[0])))\n",
        "\n",
        "print(\"With adjustment for age: {0}\".format(sum(by_res_gen[0])))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "-Nl4ZqlQSMMy"
      },
      "outputs": [],
      "source": [
        "results_dframe_genonly = pds.DataFrame(np.c_[ppm, pval_genonly, beta_genonly, p_genonly_byadj, r2_genonly, fpval_genonly],\n",
        "                               columns=['ppm', 'gender_p-value',\n",
        "                                        'gender_beta', 'gender_q-value', 'r2', 'f-test_pval'])\n",
        "results_dframe_genonly.to_csv('./Data/UnivariateAnalysis_Gender_NMR.csv', index=False)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "bS66pXxDSMMy"
      },
      "outputs": [],
      "source": [
        "results_dframe = pds.DataFrame(np.c_[ppm, pval_gen, beta_gen, p_gen_byadj, pval_age, beta_age, p_age_byadj,\n",
        "                                     r2, fpval],\n",
        "                               columns=['ppm', 'gender_p-value',\n",
        "                                        'gender_beta', 'gender_q-value', 'age_p-value', 'age_beta',\n",
        "                                        'age_q-value', 'r2', 'f-test_pval'])\n",
        "results_dframe.to_csv('./Data/UnivariateAnalysis_Gender_Age_NMR.csv', index=False)\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "g1nJZ7ocSMMz"
      },
      "source": [
        "A histogram of the *p-value* distribution shows a mostly uniform distribution, with a heavy tail at the lower end (towards 0)."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "r2rO9Uc8SMMz"
      },
      "outputs": [],
      "source": [
        "plt.figure()\n",
        "plt.hist(pval_gen, 50)\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YuIy94HhSMMz"
      },
      "source": [
        "Saving the results for comparison with the PLS-DA modelling results. The steps for this comparison are in the **Multivariate analysis with PLS-DA notebook.**"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "MVunnUk4SMMz"
      },
      "source": [
        "## Model resampling - bootstrapping"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "81dZ64GISMMz"
      },
      "source": [
        "In the multivariate modeling steps, cross-validation was used during model validation to obtain uncertainty estimates for the quality of the PLS model and its parameters.\n",
        "\n",
        "In univariate analysis using linear regression or any other statistical test, jacknifing, bootstrapping and other cross-validation and resampling techniques can also be used to obtain more reliable confidence intervals for model parameters and statistics.\n",
        "\n",
        "The confidence of the results can be further assessed by observing the bootstrap distribution for *p-values* or test statistics, such as the $\\beta$ coefficients."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jdMm_CYNSMMz"
      },
      "source": [
        "### Optional: model bootstrapping - Parallel implementation\n",
        "Non-parametric bootstrapping of cases. Run time is approximately an hour on a PC with 8 cores."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "qkBEtPxPSMMz"
      },
      "outputs": [],
      "source": [
        "# Define function that can be called by each worker:\n",
        "def bootstrap_model(variable, n_boot, dataset):\n",
        "    boot_stats = np.zeros((n_boot, 6))\n",
        "\n",
        "    for boot_iter in range(n_boot):\n",
        "        boot_sample = np.random.choice(dataset.shape[0], dataset.shape[0], replace=True)\n",
        "        fm = 'Q(dataset.columns[variable]) ~ C(Gender) + C(Age)'\n",
        "        mod = smf.ols(formula = fm, data=dataset.iloc[boot_sample, :])\n",
        "        res = mod.fit()\n",
        "        boot_stats[boot_iter, 0] = res.pvalues[1]\n",
        "        boot_stats[boot_iter, 1] = res.pvalues[2]\n",
        "        boot_stats[boot_iter, 2] = res.params[1]\n",
        "        boot_stats[boot_iter, 3] = res.params[2]\n",
        "        boot_stats[boot_iter, 4] = res.f_pvalue\n",
        "        boot_stats[boot_iter, 5] = res.rsquared\n",
        "    return boot_stats"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "hMPZDEEoSMMz"
      },
      "outputs": [],
      "source": [
        "# Skip this step if running in Google Colab and load pre-connputed results\n",
        "import joblib\n",
        "results = joblib.Parallel(n_jobs=8, verbose=5, pre_dispatch='2*n_jobs')(joblib.delayed(bootstrap_model)(i, 50, dataset) for i in range(2, dataset.shape[1]))\n",
        "# Save results.\n",
        "np.save('bootstrap_results_univariate_NMR.npy', results)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "nc12fod2SMMz"
      },
      "outputs": [],
      "source": [
        "#load pre-computed results.\n",
        "results = np.load('./Data/bootstrap_results_univariate_NMR.npy')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "aXQGKR43SMM2"
      },
      "outputs": [],
      "source": [
        "# mean and standard deviation for regression coefficients and wald test p-values\n",
        "# Gender\n",
        "mean_beta_gen = np.array([x[:, 2].mean() for x in results])\n",
        "std_beta_gen = np.array([x[:, 2].std() for x in results])\n",
        "mean_p_gen = np.array([x[:, 0].mean() for x in results])\n",
        "std_p_gen = np.array([x[:, 0].std() for x in results])\n",
        "# Age\n",
        "mean_beta_age = np.array([x[:, 3].mean() for x in results])\n",
        "std_beta_age = np.array([x[:, 3].std() for x in results])\n",
        "mean_p_age = np.array([x[:, 1].mean() for x in results])\n",
        "std_p_age = np.array([x[:, 1].std() for x in results])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9TFNnKtQSMM3"
      },
      "source": [
        "# Visualization of results\n",
        "We can now assess the regression coefficients obtained for each factor and the *p-value* corresponding to the two-tailed Wald t-test of $\\beta \\neq 0$"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "3GHjWTV9SMM3"
      },
      "outputs": [],
      "source": [
        "# Visualization of regression coefficients with Bootstrap confidence intervals\n",
        "Lineplot(mean=mean_beta_gen, error=std_beta_gen, xaxis=ppm)\n",
        "\n",
        "# Run the following line instead if bootstrapping wasn't done\n",
        "# Lineplot(mean=beta_gen, error=None, xaxis=ppm)\n",
        "plt.gca().invert_xaxis()\n",
        "plt.xlabel('$\\delta$ppm 1H')\n",
        "plt.ylabel('Beta coefficent')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Ll5LTyDOSMM3"
      },
      "source": [
        "A Manhattan plot combines the regression coefficient estimate with the magnitude of the *p-value*. The first (estimated $\\beta$) is encoded on the colorscale, while the base 10 logarithm of the *p-value* multiplied by the sign of the regression coefficient are plotted in the y axis.\n",
        "\n",
        "The next figures contain a Manhattan plot style visualization for the effect of the Age covariate..."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "tY40YCJbSMM3"
      },
      "outputs": [],
      "source": [
        "manhattan_plot(p_age_byadj, beta_age, instrument='nmr', xaxis=ppm)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "G5m1cA6iSMM3"
      },
      "source": [
        "... and for the effect of gender."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "gMLAgdb_SMM3"
      },
      "outputs": [],
      "source": [
        "manhattan_plot(p_gen_byadj, beta_gen, instrument = 'nmr', xaxis=ppm)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "VE0-ZlDTSMM3"
      },
      "source": [
        "The **interactive_manhattan** function generates an interactive version of Manhattan Plot, here exemplified for gender and age covariate."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "8gDFLnsASMM3"
      },
      "outputs": [],
      "source": [
        "iplot(interactive_manhattan(p_gen_byadj, beta_gen, sig=0.05, instrument = 'nmr', xaxis=ppm))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "F2lxCEWsSMM3"
      },
      "outputs": [],
      "source": [
        "iplot(interactive_manhattan(p_age_byadj, beta_age, sig=0.05, instrument = 'nmr', xaxis=ppm))"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python [conda env:npcproject]",
      "language": "python",
      "name": "conda-env-npcproject-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.18"
    },
    "colab": {
      "provenance": []
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}