3-DataAnalysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import shap\n",
    "from copy import deepcopy\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.metrics import balanced_accuracy_score, roc_auc_score, make_scorer, confusion_matrix, plot_confusion_matrix\n",
    "from sklearn.metrics import classification_report, auc, roc_curve, precision_recall_curve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load the model and get the data for the app\n",
    "\n",
    "fpath = \"./models/xgb_classifier_train_test_without_specialty_5y.pkl\"\n",
    "with open(fpath,\"rb\") as open_file:\n",
    "    vars = pickle.load(open_file)\n",
    "# x_train and y_train should just be x and y\n",
    "classify_xgb,X_train,y_train,X_train_ids,train_ids = vars\n",
    "\n",
    "y_train_pred = classify_xgb.predict(X_train)\n",
    "y_train_pred_prob = classify_xgb.predict_proba(X_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Figure 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "explainer = shap.TreeExplainer(classify_xgb)\n",
    "# explainer = shap.TreeExplainer(classify_xgb, feature_perturbation='tree_path_dependent', model_output=\"raw\")\n",
    "shap_values = explainer(X_train)\n",
    "shap_ival = explainer.shap_interaction_values(X_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "feature_name_dict = {\n",
    "    \"tenure\":\"Tenure\",\n",
    "    \"age_group\": \"Age group\",\n",
    "    \"EWA_avg_risk_avg\":\"Exp. weighted panel complexity\",\n",
    "    \"EWA_avg_note_quality_manual_value\": \"Exp. weighted note quality\",\n",
    "    \"EWA_avg_note_quality_contribution_value\":\"Exp. weighted note quality contribution\",\n",
    "    \"note_quality_manual_value\": \"Note quality\",\n",
    "    \"panel_cnt\": \"Panel count\",\n",
    "    \"risk_avg\": \"Panel complexity\",\n",
    "    \"EWA_avg_teamwork_on_inbox_value\":\"Exp. weighted teamwork on inbox - value\",\n",
    "    \"r_slope_panel_cnt\":'Roll. slope panel count',\n",
    "    \"teamwork_on_inbox_value\": \"Teamwork on inbox - value\",\n",
    "    \"gender\": \"Gender\",\n",
    "    \"calendar_month\": \"Calendar month\",\n",
    "    \"covid_wave\": \"Covid wave\",\n",
    "    \"patient_volume\": \"Patient volume\",\n",
    "    \"physician_demand\": \"Physician demand\",\n",
    "    'EWA_avg_order_time_8':'Exp. weighted order time',\n",
    "    'EWA_avg_wow_time_8':'Exp. weighted work outside of work time',\n",
    "    'EWA_avg_physician_demand':'Exp. weighted physician demand',\n",
    "    'EWA_avg_ib_time_8':'Exp. weighted inbox time',\n",
    "    'EWA_avg_note_time_8':'Exp. weighted note time',\n",
    "    'EWA_avg_ehr_time_8':'Exp. weighted EHR time',\n",
    "    'r_slope_wow_time_8':'Roll. slope work outside of work time',\n",
    "    'EWA_avg_patient_volume': 'Exp. weighted patient volume',\n",
    "    'physician_demand': 'physician demand',\n",
    "    \n",
    "    \n",
    "    \n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams[\"figure.figsize\"] = (20,20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# all features summarized - figure 2\n",
    "dc_shap_obj = deepcopy(shap_values)\n",
    "dc_shap_obj.feature_names = [feature_name_dict[x] if x in feature_name_dict else x for x in dc_shap_obj.feature_names]\n",
    "shap.plots.beeswarm(dc_shap_obj, max_display=10, plot_size=(20,20), order=shap_values.abs.mean(0))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dc_shap_obj = deepcopy(shap_values)\n",
    "dc_shap_obj.feature_names = [feature_name_dict[x] if x in feature_name_dict else x for x in dc_shap_obj.feature_names]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Top 4 features explored - figure 3\n",
    "plt.rcParams[\"figure.figsize\"] = (20,20)\n",
    "ax = plt.gca()\n",
    "shap.plots.scatter(dc_shap_obj[:,feature_name_dict['tenure']], color = dc_shap_obj, x_jitter=0.5, ax=ax, show=False) # male is 0\n",
    "plt.show()\n",
    "ax = plt.gca()\n",
    "shap.plots.scatter(dc_shap_obj[:,feature_name_dict['EWA_avg_risk_avg']], color = dc_shap_obj, ax=ax, show=False) # male is 0\n",
    "plt.show()\n",
    "ax = plt.gca()\n",
    "shap.plots.scatter(dc_shap_obj[:,feature_name_dict['age_group']], color = dc_shap_obj, x_jitter=0.25, ax=ax, show=False) # male is 0\n",
    "plt.xticks(ticks=[0,1,2,3,4], labels=['25-34', '35-44', '45-54', '55-64', '65+'])\n",
    "plt.show()\n",
    "ax = plt.gca()\n",
    "shap.plots.scatter(dc_shap_obj[:,feature_name_dict['EWA_avg_physician_demand']], color = dc_shap_obj, ax=ax, show=False) # male is 0\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train.columns = [feature_name_dict[x] if x in feature_name_dict else x for x in X_train.columns]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Interactions with top feature. Potentially new figure?\n",
    "# fig 4, need to justify what interactions we show\n",
    "# is there a good way to rank them in the shap documentation?\n",
    "# here i've just done some key EHR use metrics because that makes sense\n",
    "\n",
    "# shap.plots.scatter(shap_values[:,'tenure']) # male is 0\n",
    "ax = plt.gca()\n",
    "shap.dependence_plot((feature_name_dict[\"tenure\"], feature_name_dict[\"EWA_avg_ehr_time_8\"]),shap_ival, X_train,x_jitter = 0.5, ax=ax)\n",
    "plt.show()\n",
    "ax = plt.gca()\n",
    "shap.dependence_plot((feature_name_dict[\"tenure\"], feature_name_dict[\"EWA_avg_ib_time_8\"]),shap_ival, X_train,x_jitter = 0.5, ax=ax)\n",
    "plt.show()\n",
    "ax = plt.gca()\n",
    "shap.dependence_plot((feature_name_dict[\"tenure\"], feature_name_dict[\"EWA_avg_order_time_8\"]),shap_ival, X_train,x_jitter = 0.5, ax=ax)\n",
    "plt.show()\n",
    "ax = plt.gca()\n",
    "shap.dependence_plot((feature_name_dict[\"tenure\"], feature_name_dict[\"EWA_avg_note_time_8\"]),shap_ival, X_train,x_jitter = 0.5, ax=ax)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# tricky bit, going to remove some special words from the columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "'Teamwork on inbox - value'.replace(' - value', '')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train.columns = [feature_name_dict[x].replace(' - value', '') if x in feature_name_dict else x.replace(' - value', '') for x in X_train.columns]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compile_physician_data(X,y,y_pred,y_prob,X_ids,ids):\n",
    "    P = X.copy()\n",
    "    P['id'] = X_ids\n",
    "    P['prob'] = y_prob[:,1]\n",
    "    P['pred'] = y_pred\n",
    "    P['depart'] = y\n",
    "    P['phys_depart'] = P.id.isin(ids[ids['depart']]['id'])\n",
    "    # P['month_sync'] = P.groupby('id')['study_day'].transform(lambda x: round((x-max(x))/30))\n",
    "    P['month_sync'] = P.groupby('id').cumcount()\n",
    "    P['prob_rm'] = P.groupby('id')['prob'].rolling(3).mean().to_list()\n",
    "    return(P)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "P_train = compile_physician_data(X_train,y_train,y_train_pred,y_train_pred_prob,X_train_ids,train_ids)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shap_values_2 = explainer.shap_values(X_train)\n",
    "shap_ival2 = explainer.shap_interaction_values(X_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shap_ival_flat = []\n",
    "for i in range(len(shap_ival2)):\n",
    "    shap_ival_flat.append(shap_ival2[i][np.triu_indices(75, k=0)])\n",
    "shap_ival_flat = np.array(shap_ival_flat)\n",
    "triu_cols = X_train.columns[np.triu_indices(75, k=0)[0]] + '<>' + X_train.columns[np.triu_indices(75, k=0)[1]]\n",
    "shap_ival_flat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train_ids.index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(X_train_ids)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.logical_not(y_train_pred)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train_ids[np.logical_not(y_train_pred)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## subset the shap data set to only physicians for whom we have predictions of both stay and leave\n",
    "\n",
    "# create a id'd shap matrix, include prediction for later grouping\n",
    "shap_values_pd = pd.DataFrame(shap_ival_flat, columns = triu_cols)\n",
    "shap_values_pd['id'] = X_train_ids.values\n",
    "shap_values_pd['pred'] = y_train_pred\n",
    "\n",
    "# find ids with both depart and non-depart\n",
    "depart_ids = pd.unique(X_train_ids[y_train_pred])\n",
    "stay_ids = pd.unique(X_train_ids[np.logical_not(y_train_pred)])\n",
    "\n",
    "select_ids = set(depart_ids).intersection(set(stay_ids))\n",
    "select_ids = pd.Series(list(select_ids))\n",
    "\n",
    "# # do a check - all looks good\n",
    "# select_rows = X_train_ids.isin(select_ids)\n",
    "# phys_check = pd.DataFrame({\n",
    "#     'phys': X_train_ids[select_rows],\n",
    "#     'pred': y_train_pred[select_rows]\n",
    "# })\n",
    "\n",
    "# subset the shap matrix\n",
    "shap_values_pd_sub = shap_values_pd.loc[shap_values_pd['id'].isin(select_ids)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "depart_ids"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## find the mean per physician of shap scores for quit and nonquit months\n",
    "\n",
    "shap_paired = shap_values_pd_sub.groupby(['id','pred']).mean()\n",
    "shap_diff = shap_paired.groupby('id').diff()\n",
    "shap_diff = shap_diff.groupby('id').nth(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shap_paired"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sort columns by magnitude and plot\n",
    "shap_diff_np = shap_diff.to_numpy()\n",
    "#feat_vals = np.mean(np.abs(shap_diff_np),axis=0)\n",
    "feat_vals = np.mean(shap_diff_np,axis=0)\n",
    "feat_sort = np.argsort(feat_vals)\n",
    "\n",
    "# sort the overall means\n",
    "feat_vals_sort = feat_vals[feat_sort]\n",
    "\n",
    "# sort the individual physicians\n",
    "shap_diff_np_sort = shap_diff_np[:,feat_sort]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# individual contribution version\n",
    "\n",
    "num_feats = 15\n",
    "y_vec = np.arange(num_feats)\n",
    "x_vec = feat_vals_sort\n",
    "plt.barh(y_vec,x_vec[-num_feats:],color='#4286DE',lw=2)\n",
    "plt.yticks(y_vec,triu_cols[feat_sort][-num_feats:])\n",
    "plt.rcParams['figure.figsize'] = [20, 20]\n",
    "plt.xlabel('SHAP value change')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.10.7"
  },
  "vscode": {
   "interpreter": {
    "hash": "aebeff7af16a3e74dfb8d2ca55a8624561b514dde17729b916609a121bafd8b7"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}