select_outcome_model_for_diabetes_experiment.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Make Figure 9 and Table 2\n",
    "Select from logistic regression, decision tree, and random forest for outcome model for diabetes experiment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "from sklearn.metrics import roc_auc_score, brier_score_loss, log_loss\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.tree import DecisionTreeClassifier, DecisionTreeRegressor\n",
    "from sklearn.linear_model import LogisticRegression, Ridge\n",
    "from sklearn.ensemble import RandomForestClassifier, RandomForestRegressor\n",
    "from sklearn.inspection import PartialDependenceDisplay\n",
    "\n",
    "from helpers import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "path_to_data = # TODO where diabetes data is located\n",
    "with open(path_to_data) as f:\n",
    "    df = pickle.load(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Filter out providers with fewer than 4 samples. Similar to load_diabetes_data function in real_data_loader.py\n",
    "prv_counts = df['prv'].value_counts().sort_index() # Sorting to make deterministic\n",
    "prv_counter = 0\n",
    "prv_dict = dict()\n",
    "orig_prvs = []\n",
    "for prv in prv_counts.keys():\n",
    "    if prv_counts[prv] >= 4:\n",
    "        prv_dict[prv] = prv_counter\n",
    "        prv_counter += 1\n",
    "        orig_prvs.append(prv)\n",
    "df = df[df['prv'].isin(set(list(prv_dict.keys())))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Scale and split data\n",
    "X = df[['egfr','creatinine','heart_disease','treatment_date_sec']].values\n",
    "d_orig = df['prv'].values\n",
    "d = np.empty(d_orig.shape)\n",
    "for i in range(len(d_orig)):\n",
    "    d[i] = prv_dict[d_orig[i]]\n",
    "d = d.astype(int)\n",
    "t = df['y'].values\n",
    "train_fold_idxs, valid_fold_idxs, test_idxs = split_real_data_into_folds(d)\n",
    "train_valid_idxs = np.concatenate((train_fold_idxs[0], valid_fold_idxs[0]))\n",
    "train_valid_X = X[train_valid_idxs]\n",
    "scaler = StandardScaler()\n",
    "scaler.fit(train_valid_X)\n",
    "X = scaler.transform(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Table 2: AUCs on fold 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for fold_idx in range(4):\n",
    "    print('Fold ' + str(fold_idx))\n",
    "    X_train = X[train_fold_idxs[fold_idx]]\n",
    "    t_train = t[train_fold_idxs[fold_idx]]\n",
    "    X_valid = X[valid_fold_idxs[fold_idx]]\n",
    "    t_valid = t[valid_fold_idxs[fold_idx]]\n",
    "    X_test = X[test_idxs]\n",
    "    t_test = t[test_idxs]\n",
    "    \n",
    "    logreg_model = train_logreg(X_train, t_train, X_valid, t_valid)\n",
    "    logreg_train_pred = logreg_model.predict_proba(X_train)\n",
    "    logreg_train_auc = roc_auc_score(t_train, logreg_train_pred[:,1])\n",
    "    logreg_valid_pred = logreg_model.predict_proba(X_valid)\n",
    "    logreg_valid_auc = roc_auc_score(t_valid, logreg_valid_pred[:,1])\n",
    "    logreg_test_pred = logreg_model.predict_proba(X_test)\n",
    "    logreg_test_auc = roc_auc_score(t_test, logreg_test_pred[:,1])\n",
    "    \n",
    "    if X_train.shape[0] > 10000:\n",
    "        min_samples_leaf_options = [100, 500, 1000, 5000, 10000]\n",
    "    else:\n",
    "        min_samples_leaf_options = [10, 25, 100]\n",
    "    dectree_model = train_decision_tree(X_train, t_train, X_valid, t_valid, min_samples_leaf_options)\n",
    "    dectree_train_pred = dectree_model.predict_proba(X_train)\n",
    "    dectree_train_auc = roc_auc_score(t_train, dectree_train_pred[:,1])\n",
    "    dectree_valid_pred = dectree_model.predict_proba(X_valid)\n",
    "    dectree_valid_auc = roc_auc_score(t_valid, dectree_valid_pred[:,1])\n",
    "    dectree_test_pred = dectree_model.predict_proba(X_test)\n",
    "    dectree_test_auc = roc_auc_score(t_test, dectree_test_pred[:,1])\n",
    "    \n",
    "    randforest_model = train_random_forest(X_train, t_train, X_valid, t_valid)\n",
    "    randforest_train_pred = randforest_model.predict_proba(X_train)\n",
    "    randforest_train_auc = roc_auc_score(t_train, randforest_train_pred[:,1])\n",
    "    randforest_valid_pred = randforest_model.predict_proba(X_valid)\n",
    "    randforest_valid_auc = roc_auc_score(t_valid, randforest_valid_pred[:,1])\n",
    "    randforest_test_pred = randforest_model.predict_proba(X_test)\n",
    "    randforest_test_auc = roc_auc_score(t_test, randforest_test_pred[:,1])\n",
    "    \n",
    "    print('Logistic regression AUC: Train {0:.4f}'.format(logreg_train_auc) \n",
    "          + ', Valid {0:.4f}'.format(logreg_valid_auc) + ', Test {0:.4f}'.format(logreg_test_auc))\n",
    "    print('Decision tree AUC: Train {0:.4f}'.format(dectree_train_auc) + ', Valid {0:.4f}'.format(dectree_valid_auc)\n",
    "          + ', Test {0:.4f}'.format(dectree_test_auc))\n",
    "    print('Random forest AUC: Train {0:.4f}'.format(randforest_train_auc) \n",
    "          + ', Valid {0:.4f}'.format(randforest_valid_auc) + ', Test {0:.4f}'.format(randforest_test_auc))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Figure 9: partial dependence plot of selected random forest model with features at original scale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "unscaled_X = df[['egfr','creatinine','heart_disease','treatment_date_sec']].values\n",
    "unscaled_X_train = unscaled_X[train_fold_idxs[3]]\n",
    "t_train = t[train_fold_idxs[3]]\n",
    "unscaled_X_valid = unscaled_X[valid_fold_idxs[3]]\n",
    "t_valid = t[valid_fold_idxs[3]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(nrows=2, ncols=2, sharey=True, figsize=(6,6))\n",
    "feature_names=['eGFR', 'Creatinine', 'Heart disease', 'Treatment date']\n",
    "plt.rcParams.update({'font.size': 10})\n",
    "PartialDependenceDisplay.from_estimator(randforest_model, unscaled_X_train, [0,1,2,3], feature_names=feature_names, \n",
    "                                        ax=ax)\n",
    "plt.tight_layout()\n",
    "plt.savefig('random_forest_pdp_correct_scale.pdf')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python3 (sard_cji_clone)",
   "language": "python",
   "name": "sard_cji_clone"
  },
  "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}