{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0983ddd4",
   "metadata": {},
   "source": [
    "# Training and Evaluating a Meta-Predictor"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2871419",
   "metadata": {},
   "source": [
    "## Load packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45ab908c",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.preprocessing import StandardScaler"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d356bdcc",
   "metadata": {},
   "source": [
    "## Read in dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "92a99be0",
   "metadata": {},
   "outputs": [],
   "source": [
    "clinvarSet1 = pd.read_excel(\"data/arsaData/ClinvarSet1Processed.xlsx\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4e885912",
   "metadata": {},
   "outputs": [],
   "source": [
    "clinvarSet2 = pd.read_excel(\"data/arsaData/ClinvarSet2Processed.xlsx\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a4b4614",
   "metadata": {},
   "source": [
    "## Examine datasets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "01e9d42c",
   "metadata": {},
   "outputs": [],
   "source": [
    "clinvarSet1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e8f11c13",
   "metadata": {},
   "source": [
    "### Look at row with index 0 (first row in dataframe)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a1e586cf",
   "metadata": {},
   "source": [
    "Note:\n",
    "\n",
    "    - `dataframe.loc[x]` retreives the row in the dataframe with index (bold, first column) `x`\n",
    "    - `dataframe.iloc[x]` retreives the `x`th row (0-indexed) in the dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03e82f3f",
   "metadata": {},
   "outputs": [],
   "source": [
    "clinvarSet1.loc[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9bbd13c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "NUMERICAL_FEATURES = ['CADD_raw', 'DANN_score', 'Eigen_PC_raw',\n",
    "                      'Eigen_raw', 'FATHMM_score', 'GERP___RS',\n",
    "                      'GM12878_fitCons_score', 'GenoCanyon_score',\n",
    "                      'H1_hESC_fitCons_score', 'HUVEC_fitCons_score',\n",
    "                      'LRT_score', 'M_CAP_score', 'MetaLR_score',\n",
    "                      'MetaSVM_score', 'MutationAssessor_score',\n",
    "                      'MutationTaster_score', 'PROVEAN_score',\n",
    "                      'Polyphen2_HDIV_score', 'Polyphen2_HVAR_score',\n",
    "                      'REVEL_score', 'SIFT_score', 'SiPhy_29way_logOdds',\n",
    "                      'VEST3_score', 'fathmm_MKL_coding_score',\n",
    "                      'integrated_fitCons_score', 'phastCons100way_vertebrate',\n",
    "                      'phastCons20way_mammalian', 'phyloP100way_vertebrate',\n",
    "                      'phyloP20way_mammalian',\"gnomAD_exomes_AF\",\"P_rec_\",\"RVIS_EVS\"]\n",
    "\n",
    "CATEGORICAL_FEATURES = ['Ancestral_allele', 'Ensembl_geneid',\n",
    "                        'Pathway_ConsensusPathDB_',\n",
    "                        'Pathway_KEGG__full', 'Function_description',\n",
    "                        'MIM_phenotype_id', 'Trait_association_GWAS_',\n",
    "                        'GO_biological_process', 'GO_cellular_component',\n",
    "                        'GO_molecular_function']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a79f5fe",
   "metadata": {},
   "outputs": [],
   "source": [
    "clinvarSet1.loc[0,NUMERICAL_FEATURES]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a285240",
   "metadata": {},
   "outputs": [],
   "source": [
    "clinvarSet1.loc[0,CATEGORICAL_FEATURES]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0faa73c",
   "metadata": {},
   "outputs": [],
   "source": [
    "clinvarSet1.loc[0,EXTRA_NUMERICAL_FEATURES]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e999e9e",
   "metadata": {},
   "source": [
    "## Train-Test Split\n",
    "> We want to reserve some of our data to simulate our method's performance on unseen data. The rest of the data we will use for model training and selection"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "94d20739",
   "metadata": {},
   "source": [
    "randomly sample 80\\% of the instances to use as training data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b7ad6ea6",
   "metadata": {},
   "outputs": [],
   "source": [
    "trainIndices = np.random.choice(np.arange(0,len(clinvarSet1)),\n",
    "                                size=int(np.round(len(clinvarSet1) * .8)),\n",
    "                                replace=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4354823",
   "metadata": {},
   "source": [
    "Use the remaining 20\\% as the test data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "88108509",
   "metadata": {},
   "outputs": [],
   "source": [
    "testIndices = np.array(list(set(range(0,len(clinvarSet1))) - set(trainIndices)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "86e1e22e",
   "metadata": {},
   "outputs": [],
   "source": [
    "trainIndices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8621fee",
   "metadata": {},
   "outputs": [],
   "source": [
    "testIndices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44894cc1",
   "metadata": {},
   "outputs": [],
   "source": [
    "clinvarSet1Train = clinvarSet1.iloc[trainIndices]\n",
    "clinvarSet1Test = clinvarSet1.iloc[testIndices]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ec6d5cc",
   "metadata": {},
   "source": [
    "## Pre-process Features"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9bb6d0de",
   "metadata": {},
   "source": [
    "### Z-Score Standardize Numerical Features\n",
    "> We will transform each numerical feature such that each feature's distribution appears to have been drawn from a standard normal distribution N(0,1). As to not use the test data in our model development, we will calculate these statistics on the training data.\n",
    "\n",
    "**Resources**\n",
    "- [Wikipedia: Standard Score](https://en.wikipedia.org/wiki/Standard_score#Standardizing_in_mathematical_statistics)\n",
    "- [Sklearn: Standard Scaler](https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.StandardScaler.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62bcf64e",
   "metadata": {},
   "outputs": [],
   "source": [
    "numericalScalers = {col: StandardScaler() for col in NUMERICAL_FEATURES}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4dad5dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "for f in NUMERICAL_FEATURES:\n",
    "    numericalScalers[f].fit(clinvarSet1Train[f].values.reshape((-1,1)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2fb0b17",
   "metadata": {},
   "source": [
    "Before Scaling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8b116d53",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(clinvarSet1Train[\"CADD_raw\"].values)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d29a7cb",
   "metadata": {},
   "source": [
    "After Scaling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e9c2c28a",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(numericalScalers[\"CADD_raw\"].transform(clinvarSet1Train[\"CADD_raw\"].values.reshape((-1,1))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9feacc6",
   "metadata": {},
   "outputs": [],
   "source": [
    "for f in NUMERICAL_FEATURES:\n",
    "    clinvarSet1Train.loc[:,f] = numericalScalers[f].transform(clinvarSet1Train[f].values.reshape((-1,1)))\n",
    "    clinvarSet1Test.loc[:,f] = numericalScalers[f].transform(clinvarSet1Test[f].values.reshape((-1,1)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e1444d72",
   "metadata": {},
   "outputs": [],
   "source": [
    "clinvarSet1Train[NUMERICAL_FEATURES].mean(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "66dad5d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "clinvarSet1Test[NUMERICAL_FEATURES].mean(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b4dda946",
   "metadata": {},
   "source": [
    "### Missing Numerical Features\n",
    "> As you can see below, some of our numerical features are missing. To train a classifier, we need to first replace those missing numerical features with a value. As we have already standard scaled each feature, we can simply replace each NaN value with the columns mean, which is 0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88f41a8c",
   "metadata": {},
   "source": [
    "Names of columns with missing values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ee29026",
   "metadata": {},
   "outputs": [],
   "source": [
    "clinvarSet1Train.columns[pd.isna(clinvarSet1Train).sum(0)/len(clinvarSet1Train) > 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6289e167",
   "metadata": {},
   "outputs": [],
   "source": [
    "def replaceNan(col):\n",
    "    mask = np.isnan(col.values)\n",
    "    retval = np.zeros(len(col.values))\n",
    "    retval[~mask] = col.values[~mask]\n",
    "    return retval"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52bff3f3",
   "metadata": {},
   "outputs": [],
   "source": [
    "clinvarSet1Train.loc[:,NUMERICAL_FEATURES] = clinvarSet1Train[NUMERICAL_FEATURES].apply(replaceNan,axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "73e182ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "clinvarSet1Test.loc[:,NUMERICAL_FEATURES] = clinvarSet1Test[NUMERICAL_FEATURES].apply(replaceNan,axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f12f5ccf",
   "metadata": {},
   "source": [
    "## Pre-process Non-Numerical Features\n",
    "\n",
    "> Non-numerical data can be more complicated. Lets investigate each column"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db6922aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "clinvarSet1Train.loc[:,CATEGORICAL_FEATURES]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "73a4cc93",
   "metadata": {},
   "outputs": [],
   "source": [
    "len(CATEGORICAL_FEATURES)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "93e654a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "for column,values in clinvarSet1Train.loc[:,CATEGORICAL_FEATURES].iteritems():\n",
    "    print(column)\n",
    "    print(np.array(list(set(values))[:3]).reshape((-1,1)))\n",
    "    print(\"~~~~~~\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "803856f9",
   "metadata": {},
   "source": [
    "Many of these columns are lists of lists, where list values are separated by `;`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb692c2c",
   "metadata": {},
   "outputs": [],
   "source": [
    "class NonNumericalProcessor:\n",
    "    \"\"\"\n",
    "    Encode a categorical feature where each value is a single category or multiple categories separated by `;`\n",
    "    \"\"\"\n",
    "    def fit(self,values):\n",
    "        # Find the set of unique values\n",
    "        uniqueitems = set()\n",
    "        for v in values:\n",
    "            uniqueitems.update(v.upper().split(\";\"))\n",
    "        self.valueMap = {v:valueNum for valueNum,v in enumerate(uniqueitems)}\n",
    "\n",
    "    def transform(self,values):\n",
    "        encoded = np.zeros((len(values), len(self.valueMap)))\n",
    "        for i,v in enumerate(values):\n",
    "            for vi in v.upper().split(\";\"):\n",
    "                encoded[i,self.valueMap[vi]] = 1\n",
    "        return encoded\n",
    "    \n",
    "    def fit_transform(self,values):\n",
    "        self.fit(values)\n",
    "        return self.transform(values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2c3648b",
   "metadata": {},
   "outputs": [],
   "source": [
    "nnp = NonNumericalProcessor()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8eea0dd6",
   "metadata": {},
   "outputs": [],
   "source": [
    "nnp.fit_transform([\";\".join(list(\"ABC\")),\n",
    "                   \";\".join(list(\"CBa\")),\n",
    "                   \";\".join(list(\"AC\"))])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a25b4f9e",
   "metadata": {},
   "outputs": [],
   "source": [
    "nonNumericalProcessors = {col: NonNumericalProcessor() for col in CATEGORICAL_FEATURES}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7187ccd1",
   "metadata": {},
   "outputs": [],
   "source": [
    "for col in CATEGORICAL_FEATURES:\n",
    "    nonNumericalProcessors[col].fit(clinvarSet1Train[col].values)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "590a3bf5",
   "metadata": {},
   "source": [
    "Get the number of categories in each column"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a3a5f6d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "[len(p.valueMap) for p in nonNumericalProcessors.values()]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8fff1283",
   "metadata": {},
   "source": [
    "### Get Processed Matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f5d81b25",
   "metadata": {},
   "outputs": [],
   "source": [
    "XDev = clinvarSet1Train.loc[:,NUMERICAL_FEATURES].values\n",
    "yDev = clinvarSet1Train.loc[:,\"Y\"].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cac4668c",
   "metadata": {},
   "outputs": [],
   "source": [
    "XTest = clinvarSet1Test.loc[:,NUMERICAL_FEATURES].values\n",
    "yTest = clinvarSet1Test.loc[:,\"Y\"].values"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "76936671",
   "metadata": {},
   "source": [
    "## Train Classifier\n",
    "> First, we will train a logistic regression classifier and evaluate using 5-fold cross-validation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec20477b",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import KFold"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0317b201",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LogisticRegression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fac016a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import roc_auc_score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "39ba440a",
   "metadata": {},
   "outputs": [],
   "source": [
    "for k,(trainInd, valInd) in enumerate(KFold(shuffle=True,random_state=16).split(XDev)):\n",
    "    XTrain = XDev[trainInd]\n",
    "    yTrain = yDev[trainInd]\n",
    "    XVal = XDev[valInd]\n",
    "    yVal = yDev[valInd]\n",
    "    lr = LogisticRegression()\n",
    "    lr.fit(XTrain,yTrain)\n",
    "    print(f\"fold {k} val AUC: {roc_auc_score(yVal, lr.predict_proba(XVal)[:,1]):.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fda9e5b2",
   "metadata": {},
   "source": [
    "We are doing very well. Let us go ahead and use all the development data to train our final classifier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b6c5269e",
   "metadata": {},
   "outputs": [],
   "source": [
    "clf = LogisticRegression()\n",
    "clf.fit(XDev,yDev)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30ae82db",
   "metadata": {},
   "source": [
    "## Get final Performance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "47c84c48",
   "metadata": {},
   "outputs": [],
   "source": [
    "testPreds = clf.predict_proba(XTest)[:,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e5484262",
   "metadata": {},
   "outputs": [],
   "source": [
    "roc_auc_score(yTest,testPreds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4eacd9b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.calibration import calibration_curve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf183e2d",
   "metadata": {},
   "outputs": [],
   "source": [
    "acc,conf = calibration_curve(yTest,testPreds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce9cba0f",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(conf,acc)\n",
    "plt.plot([0,1],[0,1])\n",
    "plt.xlabel(\"Confidence\")\n",
    "plt.ylabel(\"Accuracy\")\n",
    "plt.title(\"Calibration Curve of Our Final Classifier\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0f3309f",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.9.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}