bTB-Diagnostic_2020_v4_crossVal+tuning_VetOnly.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8454912-7ef8-48d5-bc40-ab4a2d5ed557",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import geopandas as gp\n",
    "import geoplot\n",
    "import seaborn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f12f0587-9d1e-4a2a-8aa7-57062f389cc4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import time\n",
    "import shapely\n",
    "import rtree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2b85c35f-6acf-4ed9-9ad5-8101ca434a8a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sklearn\n",
    "from sklearn.ensemble import HistGradientBoostingClassifier as GBT\n",
    "from sklearn.metrics import roc_curve, auc, roc_auc_score, make_scorer\n",
    "from sklearn.inspection import permutation_importance\n",
    "from sklearn.model_selection import RandomizedSearchCV, train_test_split, cross_val_score, GridSearchCV\n",
    "from scipy.stats import randint,uniform,loguniform"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d842e20-252c-4ddf-8602-dab12137bfe5",
   "metadata": {},
   "outputs": [],
   "source": [
    "from joblib import dump, load"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "216ed604-493f-47c0-933e-69c7be17d13d",
   "metadata": {},
   "source": [
    "# Load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "66c256bb-b5f4-4bb9-8a77-b70e76263187",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Load data\n",
    "data = pd.read_csv('/Data/TB_Diagnostics/inputVars_VetOnly.csv',parse_dates=['dateOfTest'],dtype=float)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "956b91a3-a22d-46b8-aa0c-868cc802dadc",
   "metadata": {},
   "outputs": [],
   "source": [
    "min(data.dateOfTest) , max(data.dateOfTest)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "777af021-1f6a-4270-846c-9d7a90603a4e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get target feature (confirmed breakdowns) as binary class\n",
    "data_y = data.confirmedBreakdown.to_numpy().astype(bool)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19ea9381-9e2e-42c5-8a89-c6480c5b0fb5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get observed features\n",
    "data_X = data.drop(columns=['confirmedBreakdown'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19559dca-d5e5-4e73-907f-c0496474067d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert dates to float\n",
    "data_X.dateOfTest = data_X.dateOfTest.astype(int).astype(float)\n",
    "# Add Random features\n",
    "data_X['rand'] = np.random.random_sample(len(data_X))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5b09bee9-c8d9-41f2-9b00-dd2a06b28933",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Detect categorical features (<= 3 categories and explicit named features)\n",
    "named_cat_features = ['vetPractice','batchBovine','batchAvian']\n",
    "cat_features = []\n",
    "for c in data_X.columns:\n",
    "    catf = len(data_X[c].unique())<=3\n",
    "    if c in named_cat_features:\n",
    "        catf = True\n",
    "    cat_features.append(catf)\n",
    "\n",
    "# NB: this is fine for boolean features (inc. missing values)\n",
    "#      but needs a proper encoding for true categorical features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "247e72c9-11bc-47c4-af7e-4acd727e40f3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convery all to float matrix\n",
    "data_X = data_X.to_numpy()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b230243-181a-4b99-8718-904fd94f1409",
   "metadata": {},
   "source": [
    "# Training and testing sets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0332c21f-1377-45f6-8deb-6bb33117bbde",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Hold a final test set (random)\n",
    "X_train, X_test, y_train, y_test = train_test_split(data_X, data_y, test_size=0.20)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34be244f-ecdc-405f-ae24-722d23be7fdd",
   "metadata": {},
   "source": [
    "# Model scoring functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "291c4fc2-90ac-4e40-bc81-dbe5dde7f058",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Function: sensitivity(prediction,target)\n",
    "# returns sensitivity of prediction vs. target\n",
    "# Se = TP / (TP + FN)\n",
    "def sensitivity(p,t):\n",
    "    TP = (p&t).sum()\n",
    "    FN = (~p&t).sum()\n",
    "    return TP / (TP + FN)\n",
    "\n",
    "## Function: specificity(prediction,target)\n",
    "# returns specificity of prediction vs. target\n",
    "# Sp = TN / (TN + FP)\n",
    "def specificity(p,t):\n",
    "    TN = (~p&~t).sum()\n",
    "    FP = (p&~t).sum()\n",
    "    return TN / (TN + FP)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "06670e4a-081b-4b54-b4d7-b950a3899292",
   "metadata": {},
   "source": [
    "### SICCT Test performance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a4e9c32-3bdc-45a2-9ef7-3e35f89686e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "sicct = X_test[:,1].astype(bool)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b957aa88-d5c6-49a9-954e-f3830c88956d",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Sensitivity\n",
    "Se_sicct = sensitivity(sicct,y_test)\n",
    "Se_sicct"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5afe6254-69ff-447e-916d-acfe449583c3",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Specificity\n",
    "Sp_sicct = specificity(sicct,y_test)\n",
    "Sp_sicct"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "050ee423-09dc-49ab-b423-12ab4023ffec",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Accuracy\n",
    "(sicct==y_test).sum() / len(y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d8e8dc4c-0587-4922-b53e-068d866aead4",
   "metadata": {},
   "source": [
    "### Custom model scoring function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b680b7d-78f9-4d5c-900b-b9c4426c0332",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set specificity threshold to level for SICCT-only prediction\n",
    "specificity_threshold = Sp_sicct"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ffd35bd-fbee-47e5-812d-db40934a6032",
   "metadata": {},
   "outputs": [],
   "source": [
    "# define a custom score function:\n",
    "#    score keeps specificity above SICCT, maximises sensitivity\n",
    "def sensspec_score(t,p): #input: true (t) and predicted (p) classes\n",
    "    if specificity(p,t) < specificity_threshold:\n",
    "        return 0\n",
    "    else:\n",
    "        return sensitivity(p,t)\n",
    "\n",
    "custom_score = make_scorer(sensspec_score)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed335a56-8bdc-4f75-be84-464f18452c2c",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Hyperparameter tuning / cross-validation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d33da3a-cdf3-4552-8518-c9228bcaf606",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Create model\n",
    "gbt = GBT(categorical_features=cat_features, class_weight='balanced')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d5aaff0-11de-4a39-a9dc-3f2fe925c646",
   "metadata": {},
   "outputs": [],
   "source": [
    "# define parameter spaces\n",
    "param_grid = {'learning_rate':[1.0,0.1,0.01,0.001,0.0005,0.0001,0.00005,0.00001],\n",
    "              'max_leaf_nodes':[2,5,10,20,30,50,100,500,1000]}\n",
    "\n",
    "#param_dists = {'learning_rate':loguniform(0.00001,1.0),\n",
    "#               'max_leaf_nodes':randint(2,10000)}\n",
    "\n",
    "param_dists = {'learning_rate':loguniform(0.01,1.0),\n",
    "               'max_leaf_nodes':randint(2,2000)}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "acbf1509-6d32-446b-849a-aef78fbe5b72",
   "metadata": {},
   "outputs": [],
   "source": [
    "# perform grid search\n",
    "#model = GridSearchCV(gbt, param_grid, n_jobs=-1, cv=5)#, scoring='recall') #5-fold cross-validation #n_jobs -1 for all procs\n",
    "#start_time = time.time()\n",
    "#model.fit(X_train, y_train)\n",
    "#print(\"%0.2f seconds\" % (time.time() - start_time))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "002e9e98-19ad-4b73-911e-aee0aa8b002f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# perform random search\n",
    "model = RandomizedSearchCV(gbt, param_dists, n_jobs=-1, cv=10, verbose=1, n_iter=100, scoring='roc_auc')\n",
    "start_time = time.time()\n",
    "model.fit(X_train, y_train)\n",
    "print(\"%0.2f seconds\" % (time.time() - start_time))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c82261ac-4794-4946-9a3c-9250ca2abd3c",
   "metadata": {},
   "outputs": [],
   "source": [
    "model.best_params_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "43f701c4-5c64-404d-becb-589bc570b120",
   "metadata": {},
   "outputs": [],
   "source": [
    "model.score(X_test,y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "020db8a3-b880-44b0-b216-9d729c24e81d",
   "metadata": {},
   "outputs": [],
   "source": [
    "tuning_results = pd.DataFrame(model.cv_results_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10ae9194-678c-4e45-8cc3-e5bb0890b47b",
   "metadata": {},
   "outputs": [],
   "source": [
    "seaborn.relplot(tuning_results,x='param_learning_rate',y='mean_test_score',hue='param_max_leaf_nodes')\n",
    "plt.xscale('log')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe93d0e5-084e-47bb-9d34-413613301ad4",
   "metadata": {},
   "outputs": [],
   "source": [
    "seaborn.relplot(tuning_results,x='param_max_leaf_nodes',y='mean_test_score',hue='param_learning_rate')\n",
    "#plt.xscale('log')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b2268492-5484-4e06-a3a6-38caa03eb236",
   "metadata": {},
   "outputs": [],
   "source": [
    "seaborn.relplot(tuning_results,x='param_learning_rate',y='param_max_leaf_nodes', hue='mean_fit_time',size='mean_test_score',sizes=(1,200))\n",
    "#plt.axhline(model.best_params_['max_leaf_nodes'],ls='--', label='Best fit')\n",
    "#plt.axvline(model.best_params_['learning_rate'],ls='--')\n",
    "plt.scatter(model.best_params_['learning_rate'],model.best_params_['max_leaf_nodes'], marker='+', c='b', s=300)\n",
    "plt.xscale('log')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07c6a8d8-02ff-442c-9392-b1bc522b923f",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Evaluate performance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "53f99b2b-20ea-42b7-b052-8879a4e193b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Model score on testing set: (score is metric set at training time)\n",
    "model.score(X_test,y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "534591d2-c5b7-4a15-b2a1-1b66b3b8c088",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Get test predictions for more detailed evaluation:\n",
    "y_test_result = model.predict(X_test)\n",
    "y_score = model.decision_function(X_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7cb3cb9d-4a82-4092-97b1-111bcc15a410",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Sensitivity\n",
    "Se = sensitivity(y_test_result,y_test)\n",
    "Se"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2c803c1-0b0a-4d14-8661-ac1e0346d836",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Specificity\n",
    "Sp = specificity(y_test_result,y_test)\n",
    "Sp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4787cebd-043c-4a6e-8c43-4a30d8e9dfd0",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Accuracy\n",
    "(y_test_result==y_test).sum() / len(y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a178d34b-4fd9-46ef-bba1-3ea392f64a6f",
   "metadata": {},
   "source": [
    "---\n",
    "# ROC Curves"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "253e1f2f-dd84-4635-8dc2-1c47c24ff60f",
   "metadata": {},
   "outputs": [],
   "source": [
    "fpr, tpr, _ = roc_curve(y_test,y_score)\n",
    "roc_auc = auc(fpr,tpr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41a1dafb-75b4-42f6-bfac-1ce9fb521cfe",
   "metadata": {},
   "outputs": [],
   "source": [
    "roc_auc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f1ac1678-e96c-4203-b014-fe9532219038",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "lw = 2\n",
    "plt.plot(\n",
    "    fpr,\n",
    "    tpr,\n",
    "    lw=lw,\n",
    "    label=\"ROC curve, model (area = %0.2f)\" % roc_auc,\n",
    ")\n",
    "plt.plot(1-Sp_sicct,Se_sicct,'+', label=\"SICCT only\", ms='15')\n",
    "plt.plot([0, 1], [0, 1], lw=lw, linestyle=\"--\")\n",
    "plt.xlim([0.0, 1.0])\n",
    "plt.ylim([0.0, 1.0])\n",
    "plt.xlabel(\"(1 - Specificity)\")\n",
    "plt.ylabel(\"Sensitivity\")\n",
    "plt.title(\"Receiver operating characteristic\")\n",
    "plt.legend(loc=\"lower right\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ce0029f-bb51-4ee0-a579-4e57192c6266",
   "metadata": {},
   "source": [
    "---\n",
    "# Decision threshold choice"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b4ad4c6-eeb2-4b9b-b0a7-f4be6d2b0ed1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# function to apply decision threshold\n",
    "def predict_with_threshold(X, model, decision_threshold):\n",
    "    return model.predict_proba(X)[:,1]>=decision_threshold"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ded871a5-0f67-44a1-af96-269dbca425a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# try different thresholds\n",
    "thresholds = np.linspace(0.0,1.0,101)\n",
    "sens = np.zeros(len(thresholds)) #sensitivity at threshold\n",
    "spec = np.zeros(len(thresholds)) #specificity at threshold\n",
    "for x in range(len(thresholds)):\n",
    "    y_th = predict_with_threshold(X_test,model,thresholds[x])\n",
    "    sens[x] = sensitivity(y_th,y_test)\n",
    "    spec[x] = specificity(y_th,y_test)\n",
    "\n",
    "best_sens = max(sens[spec >= Sp_sicct]) #sensitivity s.t. specificity >= SICCT\n",
    "best_thresh = min(thresholds[spec >= Sp_sicct]) #threshold with max sensitivity s.t. specificity >= SICCT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6932b32c-2505-4edf-887d-529006b3454f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot thresholds\n",
    "plt.plot(thresholds,sens,label='Model sensitivity')\n",
    "plt.plot(thresholds,spec,label='Model specificity')\n",
    "best_sens_label = 'Best sensitivity = '+str(round(best_sens*100,1))+'%'\n",
    "sicct_sens_label = 'SICCT sensitivity = '+str(round(Se_sicct*100,1))+'%'\n",
    "sicct_spec_label = 'SICCT specificity = '+str(round(Sp_sicct*100,1))+'%'\n",
    "best_thresh_label = 'Best threshold = '+str(round(best_thresh,3))\n",
    "plt.axvline(best_thresh,c='k',ls='-.',label=best_thresh_label)\n",
    "plt.axhline(best_sens,c='k',ls='--',label=best_sens_label)\n",
    "plt.axhline(Se_sicct,c='tab:blue',ls=':',label=sicct_sens_label)\n",
    "plt.axhline(Sp_sicct,c='tab:orange',ls=':',label=sicct_spec_label)\n",
    "plt.xlabel('Decision Threshold')\n",
    "plt.legend(bbox_to_anchor=(1.0, 0.7))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8c12558-1b39-4a81-a6bb-896ed2c142d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Increase in sensitivity\n",
    "str(round((best_sens-Se_sicct)/Se_sicct * 100,1))+'% increase in sensitivity over SICCT alone.'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7db487d9-8f91-484c-9e62-c2f4eba7bee6",
   "metadata": {},
   "source": [
    "---\n",
    "# Save model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ea33123-351d-40d9-a3eb-c84b69ffc628",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save training / testing datasets to disk\n",
    "dump((X_train, X_test, y_train, y_test), '/Data/TB_Diagnostics/final_data_split_VetOnly.data')\n",
    "# Save model to disk\n",
    "dump(model, '/Data/TB_Diagnostics/final_model_VetOnly.model')"
   ]
  }
 ],
 "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.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}