shap_col1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING:tensorflow:From /home/kchen/.local/lib/python3.9/site-packages/tensorflow/python/compat/v2_compat.py:101: disable_resource_variables (from tensorflow.python.ops.variable_scope) is deprecated and will be removed in a future version.\n",
      "Instructions for updating:\n",
      "non-resource variables are not supported in the long term\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from sklearn.metrics import confusion_matrix, roc_auc_score, roc_curve, accuracy_score, auc, precision_recall_curve\n",
    "from sklearn.model_selection import train_test_split\n",
    "import tensorflow as tf\n",
    "import shap\n",
    "tf.compat.v1.disable_v2_behavior()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#read in the data\n",
    "pd.options.display.max_rows = 20\n",
    "pd.options.display.max_columns = 200\n",
    "data = pd.read_csv('clean_col.csv', index_col='CASEID')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#define anastomtoic leak as \"y\"\n",
    "y_names=['COL_ANASTOMOTIC']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#name the colectomy-specific variables\n",
    "cat = ['COL_INDICATION', 'COL_EMERGENT', 'COL_STEROID', 'COL_MECH_BOWEL_PREP', 'COL_ORAL_ANTIBIOTIC', 'COL_CHEMO','COL_APPROACH']\n",
    "cont = []\n",
    "features1 = cat + cont + y_names\n",
    "#narrow the dataset to only include colectomy-specific variables\n",
    "data = data[features1]\n",
    "col = {'COL_INDICATION':'Indication', 'COL_EMERGENT':'Emergent indication', 'COL_STEROID':'Steroid use', 'COL_MECH_BOWEL_PREP':'Mechanical bowel prep', 'COL_ORAL_ANTIBIOTIC':'Antibiotic bowel prep', 'COL_CHEMO':'Chemotherapy','COL_APPROACH':'Approach'}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#rename variables to make figure more readable\n",
    "data.rename(columns=col, inplace=True, errors='raise')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#define X and y\n",
    "y = data['COL_ANASTOMOTIC']\n",
    "X = data.drop(['COL_ANASTOMOTIC'], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#split data into train/validation sets\n",
    "X_train, X_valid, y_train, y_valid = train_test_split(X, y, train_size=0.8, test_size=0.2, random_state=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING:tensorflow:From /home/kchen/.local/lib/python3.9/site-packages/keras/layers/normalization/batch_normalization.py:520: _colocate_with (from tensorflow.python.framework.ops) is deprecated and will be removed in a future version.\n",
      "Instructions for updating:\n",
      "Colocations handled automatically by placer.\n"
     ]
    }
   ],
   "source": [
    "#create a model based on only the colectomy specific variables\n",
    "input_shape = [X.shape[1]]\n",
    "model4 = tf.keras.models.Sequential()\n",
    "model4.add(tf.keras.layers.Flatten(input_shape=input_shape))\n",
    "model4.add(tf.keras.layers.BatchNormalization())\n",
    "for _ in range(2):\n",
    "    model4.add(tf.keras.layers.Dense(1000))\n",
    "    model4.add(tf.keras.layers.BatchNormalization())\n",
    "    model4.add(tf.keras.layers.Dropout(0.8))\n",
    "    model4.add(tf.keras.layers.Activation(\"relu\"))\n",
    "model4.add(tf.keras.layers.Dense(1, activation=\"sigmoid\"))\n",
    "\n",
    "opt = tf.keras.optimizers.Adam(learning_rate=3e-4)\n",
    "\n",
    "metrics = [tf.keras.metrics.Recall(name='Sensitivity'), tf.keras.metrics.TrueNegatives(name='tn'), tf.keras.metrics.AUC(name='auc'), tf.keras.metrics.AUC(name='prc', curve='PR')]\n",
    "\n",
    "model4.compile(\n",
    "    optimizer=opt,\n",
    "    loss=tf.keras.losses.BinaryCrossentropy(from_logits=False),\n",
    "    metrics=metrics,)\n",
    "\n",
    "early_stopping = tf.keras.callbacks.EarlyStopping(\n",
    "    patience=10,\n",
    "    min_delta=1e-6,\n",
    "    restore_best_weights=True,)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Train on 155723 samples, validate on 38931 samples\n",
      "Epoch 1/200\n",
      "155723/155723 [==============================] - ETA: 0s - loss: 0.2527 - Sensitivity: 0.0417 - tn: 145505.0000 - auc: 0.5071 - prc: 0.0348"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "`Model.state_updates` will be removed in a future version. This property should not be used in TensorFlow 2.0, as `updates` are applied automatically.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "155723/155723 [==============================] - 4s 24us/sample - loss: 0.2527 - Sensitivity: 0.0417 - tn: 145505.0000 - auc: 0.5071 - prc: 0.0348 - val_loss: 0.1781 - val_Sensitivity: 0.0000e+00 - val_tn: 37630.0000 - val_auc: 0.6056 - val_prc: 0.0480\n",
      "Epoch 2/200\n",
      "155723/155723 [==============================] - 3s 18us/sample - loss: 0.1889 - Sensitivity: 0.0112 - tn: 149855.0000 - auc: 0.5238 - prc: 0.0378 - val_loss: 0.1684 - val_Sensitivity: 0.0000e+00 - val_tn: 37630.0000 - val_auc: 0.6164 - val_prc: 0.0520\n",
      "Epoch 3/200\n",
      "155723/155723 [==============================] - 3s 19us/sample - loss: 0.1733 - Sensitivity: 0.0025 - tn: 150194.0000 - auc: 0.5397 - prc: 0.0380 - val_loss: 0.1852 - val_Sensitivity: 0.0000e+00 - val_tn: 37630.0000 - val_auc: 0.6125 - val_prc: 0.0504\n",
      "Epoch 4/200\n",
      "155723/155723 [==============================] - 3s 19us/sample - loss: 0.1615 - Sensitivity: 0.0017 - tn: 150465.0000 - auc: 0.5566 - prc: 0.0406 - val_loss: 0.1926 - val_Sensitivity: 0.0000e+00 - val_tn: 37630.0000 - val_auc: 0.6067 - val_prc: 0.0494\n",
      "Epoch 5/200\n",
      "155723/155723 [==============================] - 3s 19us/sample - loss: 0.1557 - Sensitivity: 3.8775e-04 - tn: 150543.0000 - auc: 0.5585 - prc: 0.0408 - val_loss: 0.1857 - val_Sensitivity: 0.0000e+00 - val_tn: 37630.0000 - val_auc: 0.6039 - val_prc: 0.0497\n",
      "Epoch 6/200\n",
      "155723/155723 [==============================] - 3s 19us/sample - loss: 0.1513 - Sensitivity: 0.0000e+00 - tn: 150564.0000 - auc: 0.5719 - prc: 0.0421 - val_loss: 0.1846 - val_Sensitivity: 0.0000e+00 - val_tn: 37630.0000 - val_auc: 0.6061 - val_prc: 0.0500\n",
      "Epoch 7/200\n",
      "155723/155723 [==============================] - 3s 19us/sample - loss: 0.1493 - Sensitivity: 0.0000e+00 - tn: 150565.0000 - auc: 0.5765 - prc: 0.0431 - val_loss: 0.1812 - val_Sensitivity: 0.0000e+00 - val_tn: 37630.0000 - val_auc: 0.6066 - val_prc: 0.0495\n",
      "Epoch 8/200\n",
      "155723/155723 [==============================] - 3s 19us/sample - loss: 0.1478 - Sensitivity: 0.0000e+00 - tn: 150565.0000 - auc: 0.5824 - prc: 0.0452 - val_loss: 0.1787 - val_Sensitivity: 0.0000e+00 - val_tn: 37630.0000 - val_auc: 0.6053 - val_prc: 0.0494\n",
      "Epoch 9/200\n",
      "155723/155723 [==============================] - 3s 19us/sample - loss: 0.1474 - Sensitivity: 0.0000e+00 - tn: 150565.0000 - auc: 0.5811 - prc: 0.0446 - val_loss: 0.1756 - val_Sensitivity: 0.0000e+00 - val_tn: 37630.0000 - val_auc: 0.6031 - val_prc: 0.0499\n",
      "Epoch 10/200\n",
      "155723/155723 [==============================] - 3s 18us/sample - loss: 0.1468 - Sensitivity: 0.0000e+00 - tn: 150565.0000 - auc: 0.5844 - prc: 0.0457 - val_loss: 0.1724 - val_Sensitivity: 0.0000e+00 - val_tn: 37630.0000 - val_auc: 0.6053 - val_prc: 0.0496\n",
      "Epoch 11/200\n",
      "155723/155723 [==============================] - 3s 19us/sample - loss: 0.1465 - Sensitivity: 0.0000e+00 - tn: 150565.0000 - auc: 0.5858 - prc: 0.0456 - val_loss: 0.1758 - val_Sensitivity: 0.0000e+00 - val_tn: 37630.0000 - val_auc: 0.6049 - val_prc: 0.0495\n",
      "Epoch 12/200\n",
      "155723/155723 [==============================] - 3s 19us/sample - loss: 0.1463 - Sensitivity: 0.0000e+00 - tn: 150565.0000 - auc: 0.5857 - prc: 0.0460 - val_loss: 0.1725 - val_Sensitivity: 0.0000e+00 - val_tn: 37630.0000 - val_auc: 0.6034 - val_prc: 0.0495\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<keras.callbacks.History at 0x7f4cfb8f96a0>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model4.fit(X_train, y_train,\n",
    "    validation_data=(X_valid, y_valid),\n",
    "    batch_size=512,\n",
    "    epochs=200,\n",
    "    callbacks=[early_stopping])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "model4.save('keras_col.h5')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "model4 = tf.keras.models.load_model('keras_col.h5')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "import shap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "tf.config.list_physical_devices('GPU')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "keras is no longer supported, please use tf.keras instead.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING:tensorflow:From /home/kchen/.local/lib/python3.9/site-packages/shap/explainers/tf_utils.py:28: The name tf.keras.backend.get_session is deprecated. Please use tf.compat.v1.keras.backend.get_session instead.\n",
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmkAAAEaCAYAAABKL8AgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA7tklEQVR4nO3debgUxdn38W+JYEQWUXED5Ii4PGo0r94mkrjGLUZJiPsSDXGPJmqMW4wSXOK+RB81ajSighr3FVww8kSNS25NXGMU8aAiCCggiIpAv39UDTbDnJUznD6c3+e65jrTXdXVVd09Pfep6u4JWZYhIiIiIsWyTGtXQEREREQWpSBNREREpIAUpImIiIgUkII0ERERkQJSkCYiIiJSQArSRERERApo2daugEjegw8+mA0cOLC1qyEiIrIkhUoz1ZMmIiIiUkAK0kREREQKSEGaiIiISAEpSBMREREpIAVpIiIiIgWkIE1ERESkgBSkiYiIiBSQgjQRERGRAlKQJiIiIlJACtJERERECkhBmoiIiEgBKUgTERERKSAFaSIiIiIFpCBNREREpIAUpImIiIgUkII0ERERkQJSkCYiIiJSQArSRERERAooZFnW2nUQWSBcPFcHpIiIFFJ24rLVKjpUmqmeNBEREZECUpAmIiIiUkAK0kREREQKSEGaiIiISAEpSJMWZ2aDzWxsa9dDRESkLVOQViBm9jszy8zsZ61dFxEREWldCtIKwsyWAQ4HPgGOqPK6OlazfBEREVl8VXvghzTZLkAvYBDwkJlt7O6vAZhZBvwaGAysAzhwuLuPTeljgH8D/YHtgPHAie4+KqUPBbYBXgIOSn93NbM9gSFADVALDHX3e9MyvYHrgc2BTsArwPHu/mKpwma2B3BaWu8XwA3u/rtc+rHAycAKwB3A0e4+rwW2lYiIyFJPPWnFcQQwyt0fJgZER1ZI3wtYFXgdeMDMOuTSDwUuB1YEzgXuNbOaXPo2wESgD7CnmX0XGAGcCqxMDLZuM7PvpPzLAFcDfYHViYHdPaVeODPbFbgJGAqsAqwHjMqtry+wGjGo3ALYG9ivaZtERESk/VJPWgGY2ZrA7sRABuAG4EwzO9ndP0/zLsn1nJ0MTAO+A/wjpd/n7o+n9yPM7BfAAcSADWC8u1+S3s8xs8HA3aXeNuBhM7sXOAR43t3fA97L1fF04FhgXeAN4FfANe7+UMryKfB0rlmfA0NSz9lYM3sCMGJgKCIiIg1QT1oxHEq8Fq0U8AwHlgf2zeWpLb1x99nAFKB3pfTcdD59fFl6H+DdsnnvpPmY2SpmdrOZvWdmnwLvpzw9098a4K26m8TksqHNz4Cu9eQXERGRHAVprSzdMHAocZjyAzObROyp6sDCQ541uWU6E4OlDyql56bz6fPL0t+vsEw/vg7GzgPWAL7j7t1IwRtf/75YLbFXTURERKpAw52t7wfEAOjbwITc/E2BR8zsm2n61+kGgQnA+cA44Plc/kFmtgMwBtiHOLR4UD3rvQkYbWa3AKOBnYE9iDceAHQDZgPTzKwLcEHZ8lcBt5vZk8DjQGdgE3d/GhEREVlsCtJa35HE68leLJs/ycye5evetOuBe4i9XS8BPy4bTrwBOAG4n9gbtqe7lw9nLuDuz6TnsV1MvMh/PPBTd38uZRkCDAM+Bj5K00fkln/YzA4lXvP2V+Jw5vUsfF2aiIiINFPIsqy16yANSI/g2LquXqrUwzba3c9ZohWrgnDxXB2QIiJSSNmJVevbCpVm6po0ERERkQJSkCYiIiJSQBrulELRcKeIiBSVhjtFRERERHd3SrE8sP4oBg4c2NrVEBERaXXqSRMREREpIAVpIiIiIgWkIE1ERESkgBSkiYiIiBSQgjQRERGRAlKQJiIiIlJACtJERERECkhBmoiIiEgB6WehpFD0s1AiIrKkVfHnnhpLPwslIiIi0lYoSBMREREpIAVpIiIiIgWkIE1ERESkgBSkiYiIiBRQYYI0MxtqZqOrvI5ZZjagBcrZzszm1pNe9bY0lpkNNrOxrV0PERERaZoG7zk1szHAtsC+7n5Hbv53gOeA8e5eU60KtiR379LadRARERFpjMb2pP0HOLxs3uFpvrQjZtaxtesgIiLSHjT26W33AEeZWT93H2dmXYE9gXOBY0qZzGxZ4GRgMLAq8DpwnLt7Sg/E4O5XQF9gBnCBu1+Zighmdi5wWJr+k7v/Pi3bGRgOfBfoDIwFTnH3x1P6YOB04IpUhxWAO4Cj3X1eypMBW7v702l6W+AcYCNgPvCQuw9uaF2NFMzsMuBg4HPgSnc/P7ettgUuBDYAJgKXufu1Ke0B4Dl3PzdNvwfUuvs2afpqAHc/Ok0fDhwH9AHGpbo+1phKmtkwoGNq/4+BKcDZ7j4spQ8mbtdr0zpmABuZ2cbAJcBmqX0jgCHu/pWZ1QDvEvf1ycBqwBjgcHef3MjtJyIi0q41tiftC+KX8KFpen/g/4jBRd6ZxC/6HwArA38BHjGzHin9KGAo8AtgReD/Ac/nlt8GeA9YE/gRcJqZfS9X13uAdVPZtwF3m1nP3PJ9iQHBOsAWwN7AfpUaZGabAI8CNwBrEAOcYU1YV0O2AT5KZf8YOMHMDkjrXht4BPhTKn8wcJ6Z7Z2WHQ3smPKuD3QANjGz0nDtTilPKUA7BTgQ6AH8DrjHzPo3oa77ELfFSsCRwJ/M7Lu59BriPlkX2MLMViXu/3uAXsCAVKfflpV7cNoOfYhB4PAm1ElERKRda8rvIPwZeMzMfg8cAfyeGBQAC3rJjgV2c/dxafYNZnY8sBvxC/pXwB9KPVnA1PQqecvdr0nvnzOzfwMGPOPus1j4S/4iMzuFGIyNTPM+J/bmzAPGmtkTafkRFdpzFPBgqccoGQPQyHU1ZCKxlzADXjSz64jB2K3EIPel3LqfM7NriT2IdxIDsAvNbHlisPYoMRja1sxeBfoBf0vLHgec5e4vp+mRZvYkMTg9p5F1fc7dS+0dbWZ3p7r+I837CjjV3b8EMLOjgZdLPX/ABDM7D7gAOCtX7pnuPiktcxLwtpmt6e4fNrJeIiIi7VajgzR3f83MxgNnEIcyHyEGGyWrAF2AB9OwYklHoHd6XwO8Vc9qynvmPgO6AqSA5SLgh2ld81Navndrcmlos3z5CmqAf1VKaOS6GjI+BWgltcAe6X0f4nBg3jvEHjfc/Q0z+xjYmhik3UHchjsBqwMvuvv0tNzawFVmdkWurGWBD5pQ19oK05vlpieWArTcOr9nZtNz8wKxx6+uckvvewMK0kRERBrQ1F8UvY44PHiWu88zs3zaVGJQtKO7/7OO5WuJQ2ZNubar5ATi0NkOxOuzMjObSh0/StoIpbpUa119zSzkArUavg6c3icGgHn90vySJ4BdiHfWHknsSRtOHM7NP95jPPB7d7+zCXUrV1NhOh/kzS9LHw+MdvfdGlHuO2XraErwKCIi0m41NUi7jRhIvFiekAKZy4GLzewwd387XUP1PeDVNMR1FfE6s38Rr0VbCVi7nqAurxvwJfAx0CkNP67YxPrnXQs8b2YHEXuqlgG+4+5jWmhdawAnpZsHNiZeRH9CSrsNOMPMDiYOf25GDMR+kVt+NPC/wDh3n2xmU4g9mD8EfpLLdxkw1MzeBl4GvgFsDkx19zcbWdctzWx/4nbYlnhTyI715L8Z+I2ZHZLqP4cYhK3n7o/k8p1hZq8Rh6EvIAZ26kUTERFphCY9zNbdv3D30e4+rY4svwfuB+43s0+Bt4nXfpXWczVwHrE37lPgJeJ1Xo1xKTCdOFT2DjCbRYfpGi1dw/VDYmD0EfGGhYNacF1PEQO1ScBDwOXEgAZ3fzet+5fEQPAW4Iz8c+iIQVo3Uq9j6pF7kjh8/EyuHX8m3iV6IzAtteOMlK+x7kj1mUbcN8e4+zN1ZU7XmW0PDCJul2nAvcTewLzhxO3wPtCJr7eviIiINCBkWdZwLllqpUdwzHX3wxrK24Qya4jX3PVx9yYNb4aL5+qAFBGRJSo7sakDiy2u4uVUhflZKBERERH5moI0ERERkQLScKcUioY7RURkSdNwp4iIiIg0WquHjiJ5D6w/ioEDB7Z2NURERFqdetJERERECkhBmoiIiEgBKUgTERERKSAFaSIiIiIFpCBNREREpIAUpImIiIgUkII0ERERkQLSLw5IoegXB0RE2qcCPPW/NekXB0RERETaCgVpIiIiIgWkIE1ERESkgBSkiYiIiBSQgjQRERGRAmoTQZqZHWhmL+emx5jZ6Y3N3xLrbOKyQ81s9OKsv6WY2WAzG9va9RAREZGmqcr9rmb2O+AcYLC739TEZYcBc939sNI8dx8BjGhsGU3N3xLrFBEREWlJLd6TZmbLAIcDnwBHtHT50rrMrGNr10FERKQ9qEZP2i5AL2AQ8JCZbezur5USzSwDjgF+DmwAvE7scXvTzE4GDkz59kuLdAcOAk539/659axiZg8B2wHjgRPdfVRadnA+v5l1Bs4D9gCWB54GjnX39xq7zhScnAT8DFgTmAyc4u531bEdgpldBhwMfA5c6e7n57bDtsCFaRtMBC5z92tT2gPAc+5+bpp+D6h1923S9NUA7n50mj4cOA7oA4xL9XqsjnotJPUidgTmAz8GpgBnu/uw/LYErk3rmAFsZGYbA5cAm6X2jQCGuPtXZlYDvEsM1k8GVgPGAIe7++TG1EtERKS9q8Y1aUcAo9z9YeAV4MgKeQYDewKrAO8D/wvg7hcSv+xvcvcu6TWvjvUcClwOrAicC9ybgoNKLgO2TK++wFTgQTPr0IR1ngP8FNgb6AZsC7xV10YAtgE+AtYgBj8nmNkBAGa2NvAI8Cdg5bQ9zjOzvdOyo4EdU971gQ7AJmbWJaXvlPKUArRTiIFmD+B3wD1mlg9oG7IP8CiwEnF//cnMvptLryEGpusCW5jZqsD/AfcQA/IBqU6/LSv34LQd+hCDwOFNqJOIiEi71qI9aWa2JrA7MZABuAE408xOdvfPc1kvcvf30jLDaN6X933u/nh6P8LMfgEcQAzY8nVahtj7NdDdJ6R5xxOHY78NPNuIdgVi79++7v5Kmv1BetVlInCBu2fAi2Z2HTEYuxXYH3ip1FsFPGdm1wKHAXcSA7ALzWx5YrD2KDEY2tbMXgX6AX9Lyx4HnOXupZscRprZk8B+xMCyMZ5z99I+GG1md6e6/iPN+wo41d2/TNvjaODlUs8fMMHMzgMuAM7KlXumu09Ky5wEvG1ma7r7h42sl4iISLvV0sOdhxKDn4fS9HDikN6+wLBcvom5958BXZuxrtoK070r5OsJLEccfgPA3WeZ2WRiD0+DQVoqYwXq7zkrNz4FaPn67ZHe98nXJ3mH2OOGu79hZh8DWxODtDuIbdsJWB140d2np+XWBq4ysytyZS1L/QFkudoK05vlpieWArTcOr9nZtNz8wKxx6+uckvvewMK0kRERBrQYkFa6rE6lDj8+IGZlZI6EIfQhjWyqPmNzFdTYXpkhXxTgC9T+thU1y7AqsSh1sascwowmzjc93Yj69fXzEIuUKvh68DpfeCHZfn75eoD8ATx+r5tiduvFzHoXY001JmMB37v7nc2sl6V1FSYzgd55dtnPDDa3XdrRLnvlK2jKcGjiIhIu9WSPWk/IPYQfRuYkJu/KfCImX3T3V9tRDmTgC3NbBl3ry94GmRmOxAvSN8HMOLF/gtx9/lmdjNwtpm9AUwnXvD+JvBCY9bp7lm6WP/CdBH/68SgaaXc8Ge5NYCT0s0DGxMvoj8hpd0GnGFmBxOHPzcjBmK/yC0/mnit3jh3n2xmU4iB5Q+Bn+TyXQYMNbO3gZeBbwCbA1Pd/c066lZuSzPbn9hjty3xesEd68l/M/AbMzsk1X8OMQhbz90fyeU7w8xeI95YcAExsFMvmoiISCO05I0DRxKvE3vR3SflXo8ShxQr3UBQyfXEocWPzWy6mZUPoZXcQAx6ZgBDgD3dvXwIseTXgAP/BN4jBlA/yt0g0Jh1/o4YxNwHzCQGh/VdnP9UWs8k4vDv5cSAhlTPHwK/BD4GbgHOcPc7csuPJt6g8HhaJgOeJN6J+Uwpk7v/mTikfCMwLbXvjJSvse5I9ZlG3K7HuPszdWVO15ltT7yDtzYtdy+xNzBvOHE7vA90okIQLSIiIpWFLMsaziVLrUoP8m2BMmuI19z1cfcmDW+Gi+fqgBQRaYeyE6vyfP22IlSa2SZ+FkpERESkvVGQJiIiIlJAGu6UQtFwp4hI+6ThzkWpJ01ERESkgNp12CrF88D6oxg4cGBrV0NERKTVqSdNREREpIAUpImIiIgUkII0ERERkQJSkCYiIiJSQArSRERERApIQZqIiIhIASlIExERESkgBWkiIiIiBaSfhZJC0c9CiUiRtPOfKpIlRz8LJSIiItJWKEgTERERKSAFaSIiIiIFpCBNREREpIAUpLUSM7vGzK5czDLmmtl26f2BZvZyi1Su8rqqWr6IiIgsbInd3WlmY4ABwFdlSQPc/dUlUoklxMwyYGt3f7rK65kL7OjuY1q43GHAXHc/rCXLbQzd3SkiRaK7O2UJqXh355I++s5293OW8DoXMLOO7l4eJIqIiIgUTmH+RUg9bS8BawM7ApOBI4jR5R+BtYAngIPdfWZaZmXgQmBn4BvAk8Cv3P2jlF4L/AXYHtgCOMzMHgauAnYHZgJnpDwLeqTMbFCavw4wETjH3UektMHA6cAVwMnACsAdwNHuPi83JPiYmc0Hbq/UI1XeW5V6344Bfg5sALwODHb3N1N6V+BKYGCq95Cy8gYDp7t7/zTdETgJ+BmwZtqep7j7XWa2A3AusB4wN23XY919spmdDByYytgvFd8dOKis/M7AecAewPLA06mM93L780WgJu2fycAJ7n5/+bYQERGRRRXtmrSDgPOBFYG/ArcQA7VtiF/26wPHAphZAO4DMmBjoC8xeLm1rMzDgROArsD9wOVAP2Ig9E1gN6BDKbOZ7QTcABwPrEQMcq40s21yZfYFViMGcVsAewP7Abj7pinPzu7epYlDhoOBPYFVgPeB/82l/RFYF9gQ2AT4cb7eFZwD/DTVrRuwLfBWSvsS+CXQk7gN1iRuF9z9QmAEcFOqfxd3n1eh/MuALdOrLzAVeNDM8nX6GXAJMci7ErgpBXciIiLSgCXdk/Y7MzsxP8PdV8xN3uHuzwOY2XDgt8BF7v5JmvcQYCnv5um1o7t/mdJPBqaaWW93/yDl+7O7/yulzyH2Eu3q7pPTvNOAfXJ1OA643N2fStMvpLocDPw9zfscGJKCl7Fm9kSq14hmbpeSi3I9UcOA4en9Mqneu7n7pDTvFOAnlQpJAewxwL7u/kqa/UF6UXat3CQzu5DYm9goqT4/Awa6+4Q073jgE+DbwLMp61/d/R8p/TrgUmKgqRsQREREGrCkg7Q/NHBN2sTc+9l1zOua3q8NLAd8ZGa5LHxBHBotBWm1ubRVgE7A+Ny8/PtSudub2Qm5eR2Ap3LTk8t6lz7L1Wtx5NuaL7Mnsa21ufR36ymnJ3EY9q1KiWa2OXG4c1OgM3FIuUsT6lmqz4I6uPssM5sM9OHrIG1iLv2ztJ9aYjuJiIgs9QpzTVozjCcGMiu5+/x68uXTpgJziMNz76R5a1Uod5i7X7QYdWvpOxRL9a7h63rX1JN/CjGgXRd4u0L67cBdwN7u/qmZ7Q48mEuvb3uWyv8y1WEsgJl1AVYlDtOKiIjIYiraNWlN4cRhsyvSDQSYWc/cxe6LLhB7v24Fhqa8XYE/lGX7I/BrM9vazDqYWScz29zKuusaMIkYILWIXL3PNLPVzKwb8dq9uvJnwNXAhWa2sZkFM+ttZpukLN2AGcBMM1sLOLVC/fulYc1K5c8HbgbONrM103VmlwBvAi80v6UiIiJSsqSDtDPMbFbZa/fmFJQChR8Th+peNLOZwHPAdg0sehzwHnEo8DXgcWLP15ep3MeINxtcROzBmki8SL4pw4G/A84ys2lmdm0Tlmuo3u8SA6FXiT1flS7oz9fhDuLNFTOBMUD/lHYEcFiafw9wZ9my1xOHSz82s+llNwOU/JoYKP+TuD3XAH5Ux00GIiIi0kRL7GG2RWVm6xMDn17u/mFr16e908NsRaRI9DBbWUIK8TDbVmdm/YDVgeeJNxJcBvxdAZqIiIgUSVu+Jq25vgFcR7wm61XiBfYHtGqNRERERMq0u540d3+D+PBbERERkcJqjz1pIiIiIoXX7nrSpNgeWH8UAwcObO1qiIiItDr1pImIiIgUkII0ERERkQJSkCYiIiJSQArSRERERApIQZqIiIhIASlIExERESkgBWkiIiIiBaQgTURERKSAQpZlrV0HkQXCxXN1QIpIk2Un6tns0qaFSjPVkyYiIiJSQArSRERERApIQZqIiIhIASlIExERESkgBWlLMTN73cz2rSf9dDMbswSrJCIiIo2k22FamJn1Ay4Atga6ANMAB/Z19zlmNhg43d37V7su7r5RtdchIiIi1aGetJY3EpgIrA90BQYAj1LH7bXNYWYdzEz7TkREZCmm56S1IDNbGZgKbOTub1RIHwA8CXQCZqfZu7v7GDPbGLgE2Az4HBgBDHH3r8ysBngXOAz4DbAO0BdYDrgC+F5a5m7gt+7+eVpfLbHXbnia3g24CFgLGAOMBb7l7tvV0Z7y5Uv16OPuH5jZjqm8dYA5wL/dfceUtzNwFrAn0B14Afilu4+tbxvqOWki0hx6Tpq0cXpOWrW5+8fA68D1ZnawmW1oZiGX/ixwFDDO3buk1xgzWxX4P+AeoBex920n4LdlqzgA+D6xh24a8DAwiRiwbUkM1i6uVDczWyeVfy6wIjG4O3wxm3xzKqd7qvc5ubQ/Axukeq0OPA88ZGYdF3OdIiIi7YKCtJa3HbGX6njg38BHZnZGPlir4GDgZXe/1t3nuPsE4Lw0P+9Md5/k7nOAzYF1gRPc/bO0zOnAIXWsaz/gBXcf7u5z3f0x4L7mNjKZQ+xFW83dv3T3MQBmtgoxoDza3T9K9T0TWAP4zmKuU0REpF1Q/3ALc/epwGnAaWnIbx9ir9IE4C91LLY28D0zm56bF4AOZflqc+/7AFPc/bPcvHeAbwA9gclly/YuWx7i0GWvulvToB8T2/qqmU0BrnP3PxLbA/CKmeXzd0z1FhERkQYoSKsid58NDDOzXwHfSrPnV8g6Hhjt7rs1UGR+2feBnmbWOa0HoB/wBTClwrITgF3K5tU0sL6ZwAq56TXzie7+MrBv6rnbCnjMzF4BXktZ1nX3SnURERGRBihIa0Fm1gM4mXjR/3+BjNjbtDFwfso2CVjVzLq5+6dp3s3Ab8zsEOBW4jBiDbCeuz9Sx+peIF74f4mZ/YZ4ndnZwI3uXuni+9uBIWa2P3AncVh2EPHxIHV5EdjfzEYAywNn5NraCdgfeNjdp5rZNGIQOc/dJ5vZrcDVZna8u08wsxWB7YHH3X1WPesUERERdE1aS5sDrEq8QP8TYo/W6cCx7n5nyvMk8DjwrplNN7Nt3X0SMYAZRBySnAbcS+wZq8jd5wK7E4cx3yMGbc8DJ9aRfyywFzAEmA78Gri+gfacDswjPlJkDDHQy9sXeNPMZgEPAL939/9LaYcTA9UxZjYTeBXYmxi4ioiISAP0CA4pFD2CQ0SaQ4/gkDZOj+AQERERaSsUpImIiIgUkII0ERERkQJSkCYiIiJSQLrSUgrlgfVHMXDgwNauhoiISKtTT5qIiIhIASlIExERESkgBWkiIiIiBaQgTURERKSAFKSJiIiIFJCCNBEREZECUpAmIiIiUkAK0kREREQKKGRZ1tp1EFkgXDxXB6TIUiQ7Uc9MF2mEUGmmetJERERECkhBmoiIiEgBKUgTERERKSAFaSIiIiIFpCBNREREpIAUpNXBzMaY2emtXY+GmNlgMxvb2vUQERGRltWu7402MwNOB74HLAdMAkYCF7RmvepiZkOBrdx9x9aui4iIiFRXu+1JM7OdgKeB/wLfcvduwLbAx+lvu2RmwczadfAuIiJSBO32YbZm9jbwlLsfUkf6GOBFoAbYGZgMnODu9+fyDALOANYBJgLnuPuIlDaY2Et3FfAboDtwLXAecB2wE/AhcJi7P52WWRY4DRgM9ABeAo5z99fMbF/gFmJg/UWqwibANmk9VwAnAysAdwBHu/u8VO5awKXAVkAGPAj8xt1npvQMOB44CNgI2D6Vcy6wHjAXeAI41t0n57bPv4H+wHbAeOBEdx9lZj1S277r7v/Kba+/A4+7+9mVtjnoYbYiSxs9zFakUfQw2xIzW48YXNzaQNafAZcQA6wrgZvMrHMqYyfgBmJws1LKe6WZbZNbvi+wItCPGCD9ChgFXEQMwu4BbszlPwk4GPghsDrwFPC4mXVz978Sg6Yx7t4lvcbl1rMaMVjcAtgb2C/V8xvA34A3gLWBDYHewOVlbT0U2BfoAvwL+BL4JdAT+CawZh3LXJ7aeC5wr5nVuPs04E7gsFLGtM0HAH9BREREGtQugzRi4AEwoYF8f3X3f7j7fGLvV3dg3ZR2HHC5uz/l7vPd/QVgODHIKvkcONPd57j7y8DLwD/d/bnUyzUc6G9m3VP+nwMXuPub7v4lcBYwD9itgXp+Dgxx9y/dfSyx18tS2u5AcPch7v55CqDOAA40sw65Mi5293fcfV4q52l3/6e7z3X3ScCFwA5l673P3R9PeUYADhyQ0q4DDkhBIsSA7hF3b2ibi4iICO33xoEp6W8v4D/15JtYeuPun8X7DOiaZq0NbG9mJ+TydyD2fpVMTgFeyex8mWm6VOYMoA/wbm6d882sNs2vz+TS0GbyWVk91zKz6WXLZMTeulLQVJtPNLPNib1jmwKdiV2xXcrKqK0w3TvV/Wkz+xDYy8xuJ/Y0HtFAO0RERCRpl0Gau7+VHluxPzC6mcWMB4a5+0UtVzPeJ14DB4CZLZOm30+z5i+6SIPGA2+5+0YN5Csv+3bgLmBvd//UzHYnXsuWV1NhemRu+lpiD9osYo/gw42vtoiISPvWLoO05GjgQTP7CLjS3T80s9WAQ8j1ZtXjj8AwM3sO+AexF+2bxKFFb2adhgEnpwvsa4FTiPuoFNxMIvaKdXL3OY0s8yHgD2Z2GvC/xIBpTeDb7n5vPct1I/buzUw3HpxaIc8gM9sBGAPsQxxiPSiXfgvxRonfAzeW9faJiIhIPdrrNWm4++PEi/k3BF41s5nER3KsSgw6Glr+MeBw4k0AU4nDmJex6JBgU1wE3AY8BnwEfB/Y2d0/Tel3EnvVJpnZdDNbuxH1nJ3K2RB4kxh4PQF8q4FFjyBe+D+TeIPDnRXy3ACckMocAuzp7vnh2mnE3rhNU14RERFppHb7CA5ZPOkRHKPd/ZwG8g0lPopj58aUq0dwiCxd9AgOkUap+AgOfXqkatLw8eHohgEREZEma7fDnVJdZnYpMA540N11w4CIiEgTabhTCkXDnSJLFw13ijSKfnFAREREpK3QvzhSKA+sP4qBAwe2djVERERanXrSRERERApIQZqIiIhIASlIExERESkgBWkiIiIiBaQgTURERKSAFKSJiIiIFJCCNBEREZEC0i8OSKHoFwdEik+/IiDS4vSLAyIiIiJthYI0ERERkQJSkCYiIiJSQArSRERERApIQZqIiIhIASlIa4PMrLeZZWZWk6ZPM7MHq7i+qpYvIiIii9J91K3EzMYAo939nMUty93PXfwaRZXq1ZLli4iISOOoJ01ERESkgNST1srSkOW7wMHAb4E+wLPAz9x9YsqzOnAdsC3wEXBhWRlDga3cfcc03QUYCuwB9ATeB45096fMbL+0nrWBz4AHgBPc/TMzuxLYGhhgZqcCE9x9/QrlrwxcBuycqvAo8Gt3/ySl16b67gB8B6gFjnD3f7TENhMREWkP1JNWHPsC2wC9gBWAs3JpI4B5wFopz+AGyrqBGBztAHQDfgRMTGkzgAOAFYkB2dbA6QDu/kvgKeBsd+/i7uvXUf4IoAfwP+m1CnBLWZ5DgGOB7sDjwE0N1FlERERy1JNWHGe6+1QAM7sVOCy97wV8H+jv7jOAGWZ2JvBYpULMbFVgH2Bjd383zR5bSnf3UbnsY83samIvXqOY2ZrALsB67j4tzTsBeNPM1ij1/gHXuvvrKf164Hgz657aICIiIg1QkFYcE3PvPwO6pve909/xufR3qVtN+vtWpUQz2wkYAmwALAd0ACY3oZ59KtThnVxaqR3l7YHYJgVpIiIijaDhzuKbkP72zc2rqSd/bfq7bnmCmXUC7gNuB9Zy927AKSz8w67zG6jP+xXq0K8sTURERBaTetIKzt0/SI/FuNDMfg4sT+wJqyv/ZDO7C7jazAYTe+DWScmTiL1n09z9czPbEPhlWRGTgP71lP+hmT0GXGJmPyMGeJcAo3JDnSIiIrKY1JPWNhxADK7eJ17Yf3MD+Q8B/g38HzATuB9Y3d1nAb8gBnyzgKuAW8uWvQwwM5tuZq/XUf5PU7n/Bd4EptOE69pERESkYSHLstaug8gC4eK5OiBFCi47UYMwUkxDhw5l7NixDB8+vLWr0lSh0kz1pImIiEibceutt2JmdOnShTXWWINdd92Vp59+ulXqUltby/bbb0/nzp3ZYIMNGD16dIuWr3+HREREpF7h4rlVLb+xvbOXXnop559/Ptdccw277LILnTp14pFHHuH+++9nq622qmodK9l///0ZMGAAI0eOZOTIkey11168/fbb9OzZs0XKV0+aiIiIFN6MGTMYMmQIV111FXvssQcrrLACHTt2ZODAgVx00UUVl9l7771ZffXV6d69O9tssw2vv/71pdYjR45kww03pGvXrvTq1YuLL74YgKlTp7L77ruz4oorstJKK7H11lszf/6iDz546623eOmllzjzzDNZfvnl2XPPPfnmN7/J3Xff3WJtVpAmIiIihffss8/yxRdf8JOf/KTRy+y66668/fbbTJ48mc0224wDDzxwQdqhhx7Ktddey8yZM3nttdf4/ve/D8All1xC7969mTJlCh999BHnnnsuISx6ydjrr79Ov3796Nq164J5m2666UKB4OLScKcUygPrj2LgwIGtXQ0RESmYjz/+mFVWWYVll2186HLIIYcseD906FB69OjBjBkz6N69Ox07duSNN95g0003pUePHvTo0QOAjh07MnHiRMaPH0///v3ZeuutK5Y9a9YsunfvvtC87t27M2HChIr5m0M9aSIiIlJ4K6+8MlOnTmXu3MZdHzdv3jxOPfVU1llnHbp160ZNTQ0QhzMB7r77bkaOHEnfvn3ZdtttefbZZwE46aST6N+/PzvvvDP9+vXj/PPPr1h+ly5d+PTTTxea9+mnny7Us7a4FKSJiIhI4Q0YMIDllluO++67r1H5b731Vu6//35Gjx7NjBkzqK2tBaD06LEtttiC+++/n8mTJzNo0CD22WcfALp27coll1zCuHHjeOCBB7j00kt54oknFil/o402Yty4ccycOXPBvJdffpmNNtpo8RqaoyBNRERECq979+6cddZZHHPMMdx3333Mnj2br776ilGjRnHyyScvkn/mzJkst9xyrLzyysyePZvTTjttQdqcOXMYMWIEM2bMoGPHjnTr1o1llokh0UMPPcTYsWPJsozu3bvToUOHBWl56623Ht/61rc488wz+eKLL7j33nt55ZVX2HPPPVuszQrSREREpE34zW9+w6WXXso555xDz5496dOnD1deeSWDBg1aJO/BBx9M37596dWrFxtuuCFbbrnlQum33HILNTU1dOvWjWuuuYYRI0YA8Pbbb7PjjjvSpUsXBgwYwNFHH832229fsT6333477k6PHj049dRTueuuu1rs8RugXxyQgnnwwQcz3TggIiLtjH5xQERERKStUJAmIiIiUkAK0kREREQKSEGaiIiISAEpSBMREREpIAVpIiIiIgWkIE1ERESkgBSkiYiIiBSQgjQRERGRAlKQJiIiIlJACtJERERECkhBmoiIiEgBKUgTERERKaCQZVlr10FkgeWWW+61OXPmfNHa9Whpyy677Cpz586d2tr1aGlLY7uWxjaB2tXWLI3tWhrbBC3WrqlZlv1gkblZlumlV2Fem2++ubd2HdSu9t2upbFNalfbey2N7Voa21Ttdmm4U0RERKSAFKSJiIiIFJCCNCma61q7AlWidrUdS2ObQO1qa5bGdi2NbYIqtks3DoiIiIgUkHrSRERERApo2daugCydzGw94CZgZeBj4GB3f7ssTwfgCuAHQAac7+7XL05aG2/XGcB+wDzgK+A0d3+0rbcrt/z6wL+Aq939xOq2qPptMrN9gDOAkNJ3dPeP2nK7zGxV4EagD9AReBI41t3ntoF27QycC3wT+N/8MdbGzxn1tatVzhnVbFNu+SV6vkjrrGq7mnPOUE+aVMs1wFXuvh5wFXBthTwHAv2BdYEBwFAzq1nMtGqrZrteALZw902AQ4C/mtnyVWpHuWq2q3Riuxa4rzrVr6hqbTIzA4YCO7n7xsBWwIxqNaRMNffVacB/0jG4CbA5sEd1mrGIxW3XOOAw4KImLldt1WxXa50zqtmm1jpfQBXb1dxzhoI0aXHpv/HNgNvSrNuAzcysZ1nWfYE/u/t8d59C/EDuvZhpVVPtdrn7o+4+O+V7hfjf1spVas4CS2B/AZwKPAS8VZVGlFkCbfo1cLG7TwJw9xnuXvWHMC+BdmVAVzNbBlgO6ARMqFJzFmiJdrn7WHf/N1Cp16/NnjPqa1drnDOWwL6CJXy+gCXSrmadMxSkSTX0ASa4+zyA9PfDND9vLWB8bvq9XJ7mplVTtduVdzDwjrt/0AL1bkhV22VmmwK7AJe1eM3rVu19tSHQz8z+bmYvmdnpZhZauA2VVLtdZwPrAROBScCj7v5MSzagDi3Rrvq05XNGYy2pc0ZV29RK5wuo/r5q1jlDQZpIwZjZtsQvy/1buy6Ly8w6Em9PP6p08ltKdCAOB+4EbAvsChzUqjVqGXsTe2TWAHoB25jZXq1bJWnI0nLOWIrPF9DMc4aCNKmG94Fe6bqC0vUFa6b5ee8BfXPTa+XyNDetmqrdLsxsADAcGOTu/23R2tetmu1aA1gHGGlmtcDxwOFmVu3nJS2JY/Aud//S3WcC9wPfbtEWVFbtdv0KGJGGcmYQ27V9i7agspZoV33a8jmjXq1wzqhmm1rrfAFL5hhs8jlDd3dKi3P3yWb2b+J/dcPT33+l8fu8O4kfwHuI11EMArZezLSqqXa7zGwL4K/AXu7+UlUbk1PNdrn7e8AqpQLMbCjQpdp3ay2BY/BW4IdmdgvxPLoDcFe12lOyBNr1LvGutRfMrBOwI3BPtdpT0kLtqk9bPmfUqTXOGdVsU2udL9K6q30MNuucoSBNquUo4CYzGwJMI14vgZmNBIa4uwO3AN8BSrc4n+Xu76b3zU2rtmq262pgeeDaeCMQAAe5+6tVbE9JNdvVWqrZptsBA94A5gOPAjdUtzkLVLNdxwPXmNmrxOGZJ4E/V7c5CyxWu8xsK+J+6QYEM9sPONTjIyna7DmjgXa11jmjmm1qTdVsV7POGfrFAREREZEC0jVpIiIiIgWkIE1ERESkgBSkiYiIiBSQgjQRERGRAlKQJiIiIlJACtJEFlMIYZcQwlO56e1CCLWtWKUlJoQwLIRwfQuWVxNCyHLTPUMI40MIq9S3XMp7VAjhlpaqS1sQQtg6hDC9tevRHoUQftqUz3lLf1akftX6bDRjv58fQji7uetTkCayGEIIgfgbc79vIN8vQgivhRA+DSFMCyF4CGHfXHptCOGnFZZbZH6I3kpldSlL2y6EkIUQZqXXhyGEG0MIKy1eS1tHlmVTiA+BbGj7rgCcBQxdAtUqjCzLnsqybMXWrkddQghDQwijW7se7UG1tnUIYUwI4fSWLrfayj8brXgsXgAcE0Lo1ZyFFaSJLJ6dgU7Eh35WFELYnxhkHAp0J/7UyK+JD0tsju2BfsQHIlb6rb55WZZ1ybKsC7AVMAD4YzPXVQR/AX4eQuhWT56fAq9mWfbOEqrTQkIIHUIIOp+KyEKyLJsGjAKObM7yOqlIm5F6lU4PITyZeoleDSFsEkLYP4QwNoQwI4RwfQhh2dwya4UQ7gohTAohTAwhXBdC6JpLPzeEMC6V904I4fhcWk3qlToohPBGCGFmCOGxEMIauWoNAkZn9T8V+rvA37Msez6LPk//5T3WzE1xJPAI8cnX9X7wsywbBzwE/L/ytBDCsmmbDCqbPyyEcGN6v0MI4fnU+zclhHB7CGHVutaXttdWuentQghzy9Z5WuoJnB5CeCaEYJVLW9CGt4GpxJ8oqssg4PGyuhwXQngz7bf3QgjnhRA6pLSLQgj3leXfLuVdIU1vHEJ4NLW7tHzHlFY6Ng4NIbwBzAZWDSHsF0J4OfVyTgwhXFsqLy23egjhwXSsvpWWz0IINbk8h6de1xkhhH+FEHauq9EVtu+wEMItIYS/pO07IX0+vhVC+Gdq35MhhDVzy9SGEIaEEJ5OnwMPIWyRS6/3GAghdEz79L+p/HdCCHuF2FN8GrBd+Lpnt18d7dg2rWNG2mdH5tK2CyHMDSHsm8qeEUK4I/85rlBec84Vm4QQ/pbaOS4t3yGX/u20bWaFEJ4m/qOUX2fnEMLFIYR3QwifhBAeCSH0r6uOFeq8cgjh5hDPVZNCCDeFXA94KOtVzx2Dveva1iGEwam9p6TjcXII4ZIKx3HvXLmDQwhj0/sriT95dEYqs+Jvg4bYS/VECOGCdIx8HEI4IYTQN23TmSGEF0MI/5NbZrE+K+HrY/3P4etjfZHjJr2vd/uUtWWhYekW2u+PE89RTZdlmV56tYkXUEv8KY7/AToSf1/tHeA6YAXiD91OBg5M+b8BjCUOgy0P9ABGAn/JlflTYs9WAL4PfA7sktJqgIwY5KxC/KmPZ4A/55Z/Hji2rJ7bAbW56b2BL4BziL/XtmIdbftpQ/OBnsCXwB7EwCsDNi9b99zcdH/gv/k2l5V/IXBfbroLMAvYOk1vBWxB/Am51YG/A7fl8g8Drs9NZ8BW9dTnD2mb9SP+7NChxACsR36bV6jng8A59RwbHwE/Kpu3J7B22rf/L+U5MqVtCMwBeuby3wTckN6vCnxMDII7Ab0AB4aUHRtPpO3SKbVnV2Aj4j/A/Yk/AXNebh1PAHenY2lVYEwqpyalH048ZjdNZfww7Y/+dbS7fPsOIx7Du6Xlj0rLPwD0BjoDf2PhY7gW+BDYPLXjVGAK0K2Rx8AFqZ2bpG3dG9gkpQ0l/hNT3+d67VTnwWkdWwKfAHvn2pgRf0KnC7Aa8TzwuxY8V3RPx8cZwHJpuXHASbn0j9O26ZS2xyQW/pyPIJ4rVkt5zgTeBDpW+qxUqPMjxOO8R3o9DDxcz7mgJm2X3nVt67RNvwKuIp4D1wHeAk6rVEZumbG56THA6Q3sw6FpPYfx9edgHjC6bB88nltmcT8rw4jHzY9SGXukOvSt47NR1/YZWzZvwX5qif2e8mxOHPnoVN92rLhtm7qAXnq11iudpE7KTf8wfWjzX7R3AJel93sB75SVsTkxyOlQxzruAi5M70snsC1y6ccA/8pNvwUMLitju/yHOM3bnfhD1R+lk9eTwMZlbfsMmF72ms/CJ+aTiV8upRP/S8C1ZevO0rLTiD+YfQ0VAsOU/3+IwcqqafoQ4K169sHuwOTc9IITWpquM0gjfoHPBLYpK/PVUhupO0gbAVxdT73mANs1cPxcDNyRm34e+HV63zVt/++l6ROBv5UtvyfphJ47NrZpYJ2/BF5I73unZfrl0ndg4S+e14CDy8p4kDq+JKkcpOW/2Dun8vfOzTuahY/hWuDs3HQA3gMOaOgYSHlnAbvVkXcoDQdppwHPlM07D3i07JjOf84vAu6tp8xamnauOAB4n/RTiWnekcB/0/sD0zbJp/+B9Dkn/hOXAWvl0pcBZpA+D9QTpBH/UcyAdXPz1k/z1si1qTlB2pdA59y8w0if8fIycss0J0h7vWze5Ar7YFoLflaGkTvW07wpwI/r+GzUtX3qC9IWe7+neeumfKvWtx0rvfQD69LWTMy9n028/mpK2bzSMMjawFph0Tt8MmKPwIQQwrHE3ovexC+c5YkXqte1zs9y5UMMhOq7ViquMMseIv63RQhhA+IPIz8UQlg7S59iYi/P8PxyIXcXUQghpLoOz7LsqzT7BuD8EMKJWZbNTPPmZY28mDzLsv+EEF4i9iheCvwcuDG3zs2Bc4k9O52J26hLhaIaY5W07IMhdwcn8b/s3pUXWaAbMeCsyyL7IcRrAU8g9totS/wv97lclhuBXxBv/NgH+CDLsmdS2trA98qOnUDsJcirLVvnTsAQYANij0wH4pcVxN44iCf9kvFl5a0NXBVCuCI3b1ngAxpvwfGaZdnseNgs8rkpHyqszS2ThRDeI+2TBo6BnsSeqbeaUL9yfVh0374D/Dg3Xf45L/8cVtKUc0UfYHzus1iqQ5/0vneF9Hyd105/X0nbu6Rjroz6lPLky3wnlzaR5pucZdns3HQtDX/emqO8jrOp57hrgc9KpXU25rhoipba7934+p/nJtE1abI0G0/8j3HFstc3siybEEL4HnGo5khglRTYPEj8EmqsfxGHzhoty7I3iYFBX+KwRmN9nzgscEjpuhVi13oXYk9Ac90IDE7XUWwJ3JxLu53YW7delmXdqHyjQt4s4pd2yZq591OJJ9Edy/bHClmWnd9AuRsTt3VdFtoPIYQ+xOGVc4g9Ed2JQz75fXs7sF4IYTPif9Q35tLGE//rztezexZvxsibn1tnJ+C+VO5aaXudklvnhPR3rdzy+fel9R5Stt4uWZb9op62t4Sa0pv0z8BafB0Y1ncMTCF++a5bR7nz65if935+/Um/NH9JeR/oGxb+ps3XYUKF9Jrc+1IAsW7ZvuucZdltjVx/eZn9ytJmUvdnC+re1quGEDqX1bu0b0v/2DWn3GZroc9KU1VqR/k2hYXb31L7fWNiT+OcplZaQZoszR4COoV4UXPXEPUKIfwkpXcjDj1OAbIQwm7E6ySa4j5iN3ydQgiHhBD2DulZX+ki3aOAN7Is+6QJ6zqSeD3QBsC30mtjYnBxRBPrnXc7Mfi7gnjNyIRcWjdi1/3MEMJaxGsz6vMi8LMQQqd0ge8JpYT03+jlwMUhhHUBQghdQnzOXPkXwwIpeOxJvL6lLvex8I0FXYjntynAVyGELYGD8gtkWTYduJcYyG1JvCat5GbA0r77RghhmXSh8Q/qqUMnYo/AtCzLPg8hbEgcwimt7wPi0NH56XjsCZQ/2uAyYGiIF/qHEMLyIYStUu9rNR0SQtgsxAvKTyL2mD2c0uo8BtI+vRq4MMQbLUKIF7JvkrJMIvZmd6pn3bcBm4cQDg7xxpJvE4/1G1q0hfV7mLjvTkvH7vrEoKFUh4eIx9RJId4osRnxekoAsiybTOyBvzqkRy2EEFYMIfwklD0mp5Isyz4EHgMuScv1AC4BRmVZVuotehHYP31mehKvn8ura1svA1yQjqV+xKH8m9J6Pyb9YxDiHcrfJPbWl5fb6BsgGqklPitNVWn7/JsYxO6ePuM/AbbJpbfUft+JeI5qMgVpstRKXfzfJ/awvEn8onmCGNwAPEr8Mn6B2MuzF/FLuykeBeaGELarJ8804rDaf0IInxGvhZpOvLanUUK8m24QcHGWZZPyL2Jv4P8LDdwlWZcsy2YQ270r8XEXeUcQr2GZSbym7s4Givsl8YT+CfGan2Fl6b8H7gfuDyF8Sry4+yjqPxcdAgxL9azLLcCm6UuILMv+k1vXdGJgUalH40Ziux/NfRmStuv2xG1eS9yH91J2Z1delmWziPv5whDCLGLPXfnQ+QHEAOgD4k0ope35ZSrjz8SbOW5M63yP+GXcsZ62t4TriEH6NGBf4jVmpe3d0DHwO+K+vi/lGcPXX+p3EnuCJoV4B97aZcuSZdm7xOuVfkm8SPsW4Iwsy+5oobY1KLV1Z2Kg/xFfnxsuTenTiTdj7EvcRlcAfyor5nDiTTpjQggzidda7k0c5mqMnxK333+J56vpwMG59NOJ/1ROJG7j28uWr2tbjyceb+8Szz2PEI+xkp8Rz0UzUnvLg+PLiP+wTA8hvN7IttSrJT4rzbDI9sniI3uOIx7/nwA/IN6sUKrndBZzv4cQViQe39c0p9Jh4aFWEWmq1LtyWpZl26Tp7YhBRU0rVqtNSr1v72ZZFtJ0T+JdlVZ2PVGlZY8iXvh/UH35iiSEsAsxkFw+a6WTcYjXPZ5efj2ktH0hhMHEfdvSPWFLXBE+K80RQjiPeD1ks3oCdeOAyGLKsuwR4n+n0sJSYNa3kXmvoZn/rS4pIYRvEa+NeZV40fE5wF/b0peOyJKwtHxWsiz77eIsr+FOkZZXS9t+wn9rmk68GWJp1YM4ZDgLeBp4hTjcIiIL02cFDXeKiIiIFJJ60kREREQKSEGaiIiISAEpSBMREREpIAVpIiIiIgWkIE1ERESkgBSkiYiIiBTQ/wcJlE9POCj/YAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 576x309.6 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#generate the feature importance plot using SHAP deep explainer, sample values set to maximum to create robust results\n",
    "#explainer requires numpy array rather than pandas dataframe -> .values\n",
    "explainer = shap.DeepExplainer(model4, X_train.sample(n=5000))\n",
    "shap_values = explainer.shap_values(X_valid.sample(n=5000).values)\n",
    "shap.summary_plot(shap_values, X_valid.sample(n=5000))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "shap_values[0].shape\n",
    "vals = np.abs(shap_values[0]).mean(0)\n",
    "feature_names = X_train.columns\n",
    "feature_importance = pd.DataFrame(list(zip(feature_names, vals)),\n",
    "                                  columns=['col_name','feature_importance_vals'])\n",
    "feature_importance.sort_values(by=['feature_importance_vals'],\n",
    "                               ascending=False, inplace=True)\n",
    "feature_importance.head(10)\n",
    "feature_importance.to_excel('shap_col.xlsx')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LogisticRegression\n",
    "lr = LogisticRegression()\n",
    "lr.fit(X_train, y_train)\n",
    "coef = lr.coef_\n",
    "ors = np.exp(coef)\n",
    "ors = ors.squeeze()\n",
    "or_df = pd.DataFrame(list(zip(X_train.columns, ors)),columns=['col_name','odds ratio'])\n",
    "or_df['diff'] = abs(1 - or_df['odds ratio'])\n",
    "or_df.sort_values(by='diff', ascending=False).to_excel('ors_hep.xlsx')"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "916dbcbb3f70747c44a77c7bcd40155683ae19c65e1c03b4aa3499c5328201f1"
  },
  "kernelspec": {
   "display_name": "Python 3.9.5 64-bit",
   "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"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}