tutorials/HBR_FCON/HBR_NormativeModel_FCONdata_Tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "4GrH5Af4381L"
   },
   "source": [
    "# [Predictive Clinical Neuroscience Toolkit](https://github.com/amarquand/PCNtoolkit) \n",
    "# Hierarchical Bayesian Regression Normative Modelling and Transfer onto unseen site.\n",
    "\n",
    "This notebook will go through basic data preparation (training and testing set, [see Saige's tutorial](https://github.com/predictive-clinical-neuroscience/PCNtoolkit-demo/blob/main/tutorials/BLR_protocol/BLR_normativemodel_protocol.ipynb) on Normative Modelling for more detail), the actual training of the models, and will finally describe how to transfer the trained models onto unseen sites. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "fkYKxCSo381O"
   },
   "source": [
    "### Created by [Saige Rutherford](https://twitter.com/being_saige) \n",
    "### adapted/edited by Andre Marquand and Pierre Berthet"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ptugU56E381P"
   },
   "source": [
    "<div>\n",
    "<img src=\"https://github.com/predictive-clinical-neuroscience/PCNtoolkit-demo/blob/main/data/NormModelSetup.png?raw=1\" width=\"500\"/>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "yjRPC3Qq381P"
   },
   "source": [
    "## Step 0: Install necessary libraries & grab data files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "! pip install pcntoolkit==0.28"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "XnQzHXgT381Q"
   },
   "source": [
    "For this tutorial we will use data from the [Functional Connectom Project FCON1000](http://fcon_1000.projects.nitrc.org/) to create a multi-site dataset. \n",
    "\n",
    "The dataset contains some cortical measures (eg thickness), processed by Freesurfer 6.0, and some covariates (eg age, site, gender).\n",
    "\n",
    "First we import the required package, and create a working directory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "CBjOrufm381R"
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import pandas as pd\n",
    "import pcntoolkit as ptk\n",
    "import numpy as np\n",
    "import pickle\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "YklC6BFp381S"
   },
   "outputs": [],
   "source": [
    "processing_dir = \"<path-to-my>/PCNtoolkit-demo/HBR_demo/\"    # replace with desired working directory\n",
    "if not os.path.isdir(processing_dir):\n",
    "    os.makedirs(processing_dir)\n",
    "os.chdir(processing_dir)\n",
    "processing_dir = os.getcwd()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "zcNvTEhT381S"
   },
   "source": [
    "#### Overview\n",
    "Here we get the FCON dataset, remove the ICBM site for later transfer, assign some site id to the different scanner sites and print an overview of the left hemisphere mean raw cortical thickness as a function of age, color coded by the various sites:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "id": "XibIZLaT381T",
    "outputId": "916e8fd8-5b42-4252-b15a-cbfd6f2917e8",
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "fcon = pd.read_csv('https://raw.githubusercontent.com/predictive-clinical-neuroscience/PCNtoolkit-demo/main/data/fcon1000.csv')\n",
    "\n",
    "# extract the ICBM site for transfer\n",
    "icbm = fcon.loc[fcon['site'] == 'ICBM']\n",
    "icbm['sitenum'] = 0\n",
    "\n",
    "# remove from the training set (also Pittsburgh because it only has 3 samples)\n",
    "fcon = fcon.loc[fcon['site'] != 'ICBM']\n",
    "fcon = fcon.loc[fcon['site'] != 'Pittsburgh']\n",
    "\n",
    "sites = fcon['site'].unique()\n",
    "fcon['sitenum'] = 0\n",
    "\n",
    "f, ax = plt.subplots(figsize=(12, 12))\n",
    "\n",
    "for i,s in enumerate(sites):\n",
    "    idx = fcon['site'] == s\n",
    "    fcon['sitenum'].loc[idx] = i\n",
    "    \n",
    "    print('site',s, sum(idx))\n",
    "    ax.scatter(fcon['age'].loc[idx], fcon['lh_MeanThickness_thickness'].loc[idx])\n",
    "    \n",
    "ax.legend(sites)\n",
    "ax.set_ylabel('LH mean cortical thickness [mm]')\n",
    "ax.set_xlabel('age')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "R41KrcwG381U"
   },
   "source": [
    "## Step 1: Prepare training and testing sets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "FNF8I4_S381U"
   },
   "source": [
    "Then we randomly split half of the samples (participants) to be either in the training or in the testing samples. We do this for the remaing FCON dataset and for the ICBM data. The transfer function will also require a training and a test sample. \n",
    "\n",
    "The numbers of samples per sites used for training and for testing are then displayed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "txJK1fHt381V",
    "outputId": "5aba80e6-d783-444c-bb20-3578ef3def0c"
   },
   "outputs": [],
   "source": [
    "tr = np.random.uniform(size=fcon.shape[0]) > 0.5\n",
    "te = ~tr\n",
    "    \n",
    "fcon_tr = fcon.loc[tr]\n",
    "fcon_te = fcon.loc[te]\n",
    "    \n",
    "tr = np.random.uniform(size=icbm.shape[0]) > 0.5\n",
    "te = ~tr\n",
    "    \n",
    "icbm_tr = icbm.loc[tr]\n",
    "icbm_te = icbm.loc[te]\n",
    "    \n",
    "print('sample size check')\n",
    "for i,s in enumerate(sites):\n",
    "    idx = fcon_tr['site'] == s\n",
    "    idxte = fcon_te['site'] == s\n",
    "    print(i,s, sum(idx), sum(idxte))\n",
    "    \n",
    "fcon_tr.to_csv(processing_dir + '/fcon1000_tr.csv')\n",
    "fcon_te.to_csv(processing_dir + '/fcon1000_te.csv')\n",
    "icbm_tr.to_csv(processing_dir + '/fcon1000_icbm_tr.csv')\n",
    "icbm_te.to_csv(processing_dir + '/fcon1000_icbm_te.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "QlTN6JAW381V"
   },
   "source": [
    "Otherwise you can just load these pre defined subsets:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "uXBxuLK7381V"
   },
   "outputs": [],
   "source": [
    "# Optional\n",
    "#fcon_tr = pd.read_csv('https://raw.githubusercontent.com/predictive-clinical-neuroscience/PCNtoolkit-demo/main/data/fcon1000_tr.csv')\n",
    "#fcon_te = pd.read_csv('https://raw.githubusercontent.com/predictive-clinical-neuroscience/PCNtoolkit-demo/main/data/fcon1000_te.csv')\n",
    "#icbm_tr = pd.read_csv('https://raw.githubusercontent.com/predictive-clinical-neuroscience/PCNtoolkit-demo/main/data/fcon1000_icbm_tr.csv')\n",
    "#icbm_te = pd.read_csv('https://raw.githubusercontent.com/predictive-clinical-neuroscience/PCNtoolkit-demo/main/data/fcon1000_icbm_te.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "qwp5XfII381W"
   },
   "source": [
    "## Step 2: Configure HBR inputs: covariates, measures and batch effects"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "NnIR5Sqq381W"
   },
   "source": [
    "We will here only use the mean cortical thickness for the Right and Left hemisphere: two idps.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "4_EqQFOe381W"
   },
   "outputs": [],
   "source": [
    "idps = ['rh_MeanThickness_thickness','lh_MeanThickness_thickness']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "wLu-CoUl381W"
   },
   "source": [
    "As input to the model, we need covariates (used to describe predictable source of variability (fixed effects), here 'age'), measures (here cortical thickness on two idps), and batch effects (random source of variability, here 'scanner site' and 'sex').\n",
    "\n",
    "`X` corresponds to the covariate(s) \n",
    "\n",
    "`Y` to the measure(s)\n",
    "\n",
    "`batch_effects` to the random effects\n",
    "\n",
    "We need these values both for the training (`_train`) and for the testing set (`_test`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "W2xf-Zc_381X"
   },
   "outputs": [],
   "source": [
    "X_train = (fcon_tr['age']/100).to_numpy(dtype=float)\n",
    "Y_train = fcon_tr[idps].to_numpy(dtype=float)\n",
    "\n",
    "# configure batch effects for site and sex\n",
    "#batch_effects_train = fcon_tr[['sitenum','sex']].to_numpy(dtype=int)\n",
    "\n",
    "# or only site\n",
    "batch_effects_train = fcon_tr[['sitenum']].to_numpy(dtype=int)\n",
    "    \n",
    "with open('X_train.pkl', 'wb') as file:\n",
    "    pickle.dump(pd.DataFrame(X_train), file)\n",
    "with open('Y_train.pkl', 'wb') as file:\n",
    "    pickle.dump(pd.DataFrame(Y_train), file) \n",
    "with open('trbefile.pkl', 'wb') as file:\n",
    "    pickle.dump(pd.DataFrame(batch_effects_train), file) \n",
    "\n",
    "\n",
    "X_test = (fcon_te['age']/100).to_numpy(dtype=float)\n",
    "Y_test = fcon_te[idps].to_numpy(dtype=float)\n",
    "#batch_effects_test = fcon_te[['sitenum','sex']].to_numpy(dtype=int)\n",
    "batch_effects_test = fcon_te[['sitenum']].to_numpy(dtype=int)\n",
    "    \n",
    "with open('X_test.pkl', 'wb') as file:\n",
    "    pickle.dump(pd.DataFrame(X_test), file)\n",
    "with open('Y_test.pkl', 'wb') as file:\n",
    "    pickle.dump(pd.DataFrame(Y_test), file) \n",
    "with open('tsbefile.pkl', 'wb') as file:\n",
    "    pickle.dump(pd.DataFrame(batch_effects_test), file) \n",
    "\n",
    "# a simple function to quickly load pickle files    \n",
    "def ldpkl(filename: str): \n",
    "    with open(filename, 'rb') as f:\n",
    "        return pickle.load(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "batch_effects_test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "6WYLCEmV381X"
   },
   "source": [
    "## Step 3: Files and Folders grooming"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "I18buNGh381X"
   },
   "outputs": [],
   "source": [
    "respfile = os.path.join(processing_dir, 'Y_train.pkl')       # measurements  (eg cortical thickness) of the training samples (columns: the various features/ROIs, rows: observations or subjects)\n",
    "covfile = os.path.join(processing_dir, 'X_train.pkl')        # covariates (eg age) the training samples (columns: covariates, rows: observations or subjects)\n",
    "\n",
    "testrespfile_path = os.path.join(processing_dir, 'Y_test.pkl')       # measurements  for the testing samples\n",
    "testcovfile_path = os.path.join(processing_dir, 'X_test.pkl')        # covariate file for the testing samples\n",
    "\n",
    "trbefile = os.path.join(processing_dir, 'trbefile.pkl')      # training batch effects file (eg scanner_id, gender)  (columns: the various batch effects, rows: observations or subjects)\n",
    "tsbefile = os.path.join(processing_dir, 'tsbefile.pkl')      # testing batch effects file\n",
    "\n",
    "output_path = os.path.join(processing_dir, 'Models/')    #  output path, where the models will be written\n",
    "log_dir = os.path.join(processing_dir, 'log/')           #\n",
    "if not os.path.isdir(output_path):\n",
    "    os.mkdir(output_path)\n",
    "if not os.path.isdir(log_dir):\n",
    "    os.mkdir(log_dir)\n",
    "\n",
    "outputsuffix = '_estimate'      # a string to name the output files, of use only to you, so adapt it for your needs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "wykFbsV_381Y"
   },
   "source": [
    "## Step 4: Estimating the models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "HvcGW4vK381Y"
   },
   "source": [
    "Now we have everything ready to estimate the normative models. The `estimate` function only needs the training and testing sets, each divided in three datasets: covariates, measures and batch effects. We obviously specify `alg=hbr` to use the hierarchical bayesian regression method, well suited for the multi sites datasets. The remaining arguments are basic data management: where the models, logs, and output files will be written and how they will be named."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "fEYt38W9381Y",
    "outputId": "66664835-e06a-4ed9-c41c-6225dfb6d7d0"
   },
   "outputs": [],
   "source": [
    "ptk.normative.estimate(covfile=covfile, \n",
    "                       respfile=respfile,\n",
    "                       tsbefile=tsbefile, \n",
    "                       trbefile=trbefile, \n",
    "                       inscaler='standardize',\n",
    "                       outscaler='standardize',\n",
    "                       linear_mu='True',\n",
    "                       random_intercept_mu='True',\n",
    "                       centered_intercept_mu='True',\n",
    "                       alg='hbr', \n",
    "                       log_path=log_dir, \n",
    "                       binary=True,\n",
    "                       output_path=output_path, \n",
    "                       testcov= testcovfile_path,\n",
    "                       testresp = testrespfile_path,\n",
    "                       outputsuffix=outputsuffix, \n",
    "                       savemodel=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "cu4Umis3381Z"
   },
   "source": [
    "Here some analyses can be done, there are also some error metrics that could be of interest. This is covered in step 6 and in [Saige's tutorial](https://github.com/predictive-clinical-neuroscience/PCNtoolkit-demo/blob/main/tutorials/BLR_protocol/BLR_normativemodel_protocol.ipynb) on Normative Modelling."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "FyTrHWij381Z"
   },
   "source": [
    "## Step 5: Transfering the models to unseen sites"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ZRmYyQry381Z"
   },
   "source": [
    "Similarly to what was done before for the FCON data, we also need to prepare the ICBM specific data, in order to run the transfer function: training and testing set of covariates, measures and batch effects:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "ujwfzth4381Z"
   },
   "outputs": [],
   "source": [
    "X_adapt = (icbm_tr['age']/100).to_numpy(dtype=float)\n",
    "Y_adapt = icbm_tr[idps].to_numpy(dtype=float)\n",
    "#batch_effects_adapt = icbm_tr[['sitenum','sex']].to_numpy(dtype=int)\n",
    "batch_effects_adapt = icbm_tr[['sitenum']].to_numpy(dtype=int)\n",
    "    \n",
    "with open('X_adaptation.pkl', 'wb') as file:\n",
    "    pickle.dump(pd.DataFrame(X_adapt), file)\n",
    "with open('Y_adaptation.pkl', 'wb') as file:\n",
    "    pickle.dump(pd.DataFrame(Y_adapt), file) \n",
    "with open('adbefile.pkl', 'wb') as file:\n",
    "    pickle.dump(pd.DataFrame(batch_effects_adapt), file) \n",
    "\n",
    "# Test data (new dataset)\n",
    "X_test_txfr = (icbm_te['age']/100).to_numpy(dtype=float)\n",
    "Y_test_txfr = icbm_te[idps].to_numpy(dtype=float)\n",
    "#batch_effects_test_txfr = icbm_te[['sitenum','sex']].to_numpy(dtype=int)\n",
    "batch_effects_test_txfr = icbm_te[['sitenum']].to_numpy(dtype=int)\n",
    "    \n",
    "with open('X_test_txfr.pkl', 'wb') as file:\n",
    "    pickle.dump(pd.DataFrame(X_test_txfr), file)\n",
    "with open('Y_test_txfr.pkl', 'wb') as file:\n",
    "    pickle.dump(pd.DataFrame(Y_test_txfr), file) \n",
    "with open('txbefile.pkl', 'wb') as file:\n",
    "    pickle.dump(pd.DataFrame(batch_effects_test_txfr), file) \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "wqzmwQ0B381Z"
   },
   "outputs": [],
   "source": [
    "respfile = os.path.join(processing_dir, 'Y_adaptation.pkl')\n",
    "covfile = os.path.join(processing_dir, 'X_adaptation.pkl')\n",
    "testrespfile_path = os.path.join(processing_dir, 'Y_test_txfr.pkl')\n",
    "testcovfile_path = os.path.join(processing_dir, 'X_test_txfr.pkl')\n",
    "trbefile = os.path.join(processing_dir, 'adbefile.pkl')\n",
    "tsbefile = os.path.join(processing_dir, 'txbefile.pkl')\n",
    "\n",
    "log_dir = os.path.join(processing_dir, 'log_transfer/')\n",
    "output_path = os.path.join(processing_dir, 'Transfer/')\n",
    "model_path = os.path.join(processing_dir, 'Models/')  # path to the previously trained models\n",
    "outputsuffix = '_transfer'  # suffix added to the output files from the transfer function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "VlP1kyjj381a"
   },
   "source": [
    "Here, the difference is that the transfer function needs a model path, which points to the models we just trained, and new site data (training and testing). That is basically the only difference."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "hwpmdprx381a"
   },
   "outputs": [],
   "source": [
    "yhat, s2, z_scores = ptk.normative.transfer(covfile=covfile, \n",
    "                                            respfile=respfile,\n",
    "                                            tsbefile=tsbefile, \n",
    "                                            trbefile=trbefile,\n",
    "                                            inscaler='standardize',\n",
    "                                            outscaler='standardize',\n",
    "                                            linear_mu='True',\n",
    "                                            random_intercept_mu='True',\n",
    "                                            centered_intercept_mu='True',\n",
    "                                            model_path = model_path,\n",
    "                                            alg='hbr', \n",
    "                                            log_path=log_dir, \n",
    "                                            binary=True,\n",
    "                                            output_path=output_path, \n",
    "                                            testcov= testcovfile_path,\n",
    "                                            testresp = testrespfile_path,\n",
    "                                            outputsuffix=outputsuffix, \n",
    "                                            savemodel=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "output_path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "EV = pd.read_pickle('EXPV_estimate.pkl')\n",
    "print(EV)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "QBtsSlDG381a"
   },
   "source": [
    "And that is it, you now have models that benefited from prior knowledge about different scanner sites to learn on unseen sites."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "INPI9Tnt381a"
   },
   "source": [
    "## Step 6: Interpreting model performance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "eMCWLa3b381a"
   },
   "source": [
    "Output evaluation metrics definitions: \n",
    "* yhat - predictive mean\n",
    "* ys2 - predictive variance\n",
    "* nm - normative model\n",
    "* Z - deviance scores\n",
    "* Rho - Pearson correlation between true and predicted responses\n",
    "* pRho - parametric p-value for this correlation\n",
    "* RMSE - root mean squared error between true/predicted responses\n",
    "* SMSE - standardised mean squared error\n",
    "* EV - explained variance\n",
    "* MSLL - mean standardized log loss\n",
    "    * See page 23 in http://www.gaussianprocess.org/gpml/chapters/RW2.pdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "ux4ZjZtj381b"
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "name": "HBR_NormativeModel_FCONdata_Tutorial.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "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.10.9"
  },
  "vscode": {
   "interpreter": {
    "hash": "40d3a090f54c6569ab1632332b64b2c03c39dcf918b08424e98f38b5ae0af88f"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}