4-ShapQuantification.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "127c8cde-d554-4ca5-aed2-051301fb5a58",
   "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\n",
    "from statsmodels.stats.diagnostic import linear_reset\n",
    "import statsmodels.api as sm\n",
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b882ed2d-7b6a-4492-adfb-6078f7c6cbbb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# load the model and get the data for the app\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": "code",
   "execution_count": null,
   "id": "af892163-37c7-4e80-aa29-6e43ea4b52d7",
   "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",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e26f6239-922f-4151-9c1f-c5261039a58f",
   "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)\n",
    "# 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",
    "dc_shap_ival_obj = deepcopy(shap_ival)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d0bd138-1c09-4288-812f-e1dc68082b65",
   "metadata": {},
   "source": [
    "### Tenure Tests\n",
    "Here we fit a linear model to two variables: tenure and EWA_avg_ehr_time_8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "236a733c-22d8-4c5f-953c-3c9a439f6414",
   "metadata": {},
   "outputs": [],
   "source": [
    "class shap_quant:\n",
    "    def __init__(self, shap_obj, shap_ival_obj, df, feature_name_dict, main_feat, interaction_feat):\n",
    "        self.shap_obj = shap_obj\n",
    "        self.df = df\n",
    "        self.feature_name_dict = feature_name_dict\n",
    "        self.main_feat = main_feat\n",
    "        self.interaction_feat = interaction_feat\n",
    "        self.SHO_main_feat = feature_name_dict[main_feat] if main_feat in feature_name_dict else main_feat\n",
    "        self.SHO_interaction_feat = feature_name_dict[interaction_feat] if interaction_feat in feature_name_dict else interaction_feat\n",
    "        self.shap_ival_obj = shap_ival_obj\n",
    "        self.fit_ols()\n",
    "    def fit_ols(self):\n",
    "        self.shap_value = self.shap_obj[:,self.SHO_main_feat].values\n",
    "        self.feature_value = self.df[self.main_feat].to_numpy()\n",
    "        self.interaction_value = self.df[self.interaction_feat].to_numpy()\n",
    "        model_ols = sm.OLS(self.feature_value,self.shap_value,missing=\"drop\")\n",
    "        model_ols_fit = model_ols.fit()\n",
    "        self.model_ols_fit = model_ols_fit\n",
    "        return model_ols_fit\n",
    "    def is_linear(self) :# is it linear?\n",
    "        res = linear_reset(self.model_ols_fit)\n",
    "        print('Linear test: ', res)\n",
    "        #return res\n",
    "    '''\n",
    "    # here is the plot of tenure shap value vs. value of the interaction feature\n",
    "    # plt.scatter(tenure_sv[tenure_fv==5.0],tenure_iv[tenure_fv==5.0])\n",
    "    # When sv=shap value, fv=feature value, iv=interaction value\n",
    "    # so when the fv -s a value, x, we find the correlation coeeficient between the two to see if there is a positive or negative relationship\n",
    "    # works by ranking the values of x and y, (lowest value = 1, highest value = n)\n",
    "    # push this into the formula and you get the coeficcients.\n",
    "    # so these all seem to have a positive relationship\n",
    "    # It assesses how well the relationship between two variables can be described using a monotonic function.\n",
    "    # so it looks at both the linear relationship between the variables but also how well it can be described by some monotonic function\n",
    "    # if a var goes up either the other also goes up or the other goes down\n",
    "    # here we can say as the relationship \n",
    "    '''\n",
    "    def sv_iv_correlation(self): # check if the shap value and interaction value between the main feat and interaction feat\n",
    "        vals = [x for x in np.unique(self.feature_value).tolist() if isinstance(x, (int,float)) and not np.isnan(x)]\n",
    "        for v in vals:\n",
    "            print(f'for value {v}', stats.spearmanr(self.shap_value[self.feature_value==v],self.interaction_value[self.feature_value==v]))\n",
    "    '''\n",
    "    # look at the correlation between the interaction shap value (i.e. the y\n",
    "    # coordinate on plot) and the interaction feature value (i.e. the color of the\n",
    "    # point) for a given primary feature value (i.e. given tenure)\n",
    "    '''\n",
    "    def isv_ifv_correlation(self):\n",
    "        idx1 = X_train.columns.get_loc(self.main_feat)\n",
    "        idx2 = X_train.columns.get_loc(self.interaction_feat)\n",
    "        isv = self.shap_ival_obj[:,idx1,idx2] #isv interaction shap value\n",
    "        fv = self.df[self.main_feat]\n",
    "        ifv = self.df[self.interaction_feat] # ifv interaction feature value\n",
    "        vals = [x for x in np.unique(self.feature_value).tolist() if isinstance(x, (int,float)) and not np.isnan(x)]\n",
    "        for v in vals:\n",
    "            print(f'for value {v}', stats.spearmanr(isv[fv==v],ifv[fv==v]))\n",
    "        plt.scatter(fv,isv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6da17b95-5a2d-426a-b1c8-91eea9a23ead",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "sq = shap_quant(dc_shap_obj,dc_shap_ival_obj, X_train, feature_name_dict, 'tenure', 'EWA_avg_ehr_time_8')\n",
    "sq.is_linear()\n",
    "sq.sv_iv_correlation()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d1361660-3e24-44f2-b4f8-349d68490926",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "sq = shap_quant(dc_shap_obj,dc_shap_ival_obj, X_train, feature_name_dict, 'EWA_avg_risk_avg', 'EWA_avg_physician_demand')\n",
    "sq.is_linear()\n",
    "sq.sv_iv_correlation()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "00b5e598-07b7-42d6-a1ac-30ba908de109",
   "metadata": {},
   "outputs": [],
   "source": [
    "sq = shap_quant(dc_shap_obj,dc_shap_ival_obj, X_train, feature_name_dict, 'age_group', 'EWA_avg_ehr_time_8')\n",
    "sq.is_linear()\n",
    "sq.sv_iv_correlation()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e403d38b-8324-420f-b6fd-ce09343824e8",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "sq = shap_quant(dc_shap_obj,dc_shap_ival_obj, X_train, feature_name_dict, 'EWA_avg_physician_demand', 'panel_cnt')\n",
    "sq.is_linear()\n",
    "sq.sv_iv_correlation()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05e9ea47-ead2-412c-aa65-7834b019221a",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "sq = shap_quant(dc_shap_obj,dc_shap_ival_obj, X_train, feature_name_dict, 'EWA_avg_physician_demand', 'panel_cnt')\n",
    "sq.is_linear()\n",
    "sq.sv_iv_correlation()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0416aa71-620b-440f-bdb3-70e834dfad80",
   "metadata": {},
   "source": [
    "### lets look at the interaction shap and interaction feature values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "265362da-7ae8-41a1-b564-3e69bd51907f",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "sq = shap_quant(dc_shap_obj,dc_shap_ival_obj, X_train, feature_name_dict, 'tenure', 'EWA_avg_ehr_time_8')\n",
    "sq.is_linear()\n",
    "sq.isv_ifv_correlation()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61ef8df6-ef08-4889-9ac4-82f741100946",
   "metadata": {},
   "outputs": [],
   "source": [
    "sq = shap_quant(dc_shap_obj,dc_shap_ival_obj, X_train, feature_name_dict, 'tenure', 'EWA_avg_ib_time_8')\n",
    "sq.is_linear()\n",
    "sq.isv_ifv_correlation()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0f822961-7dea-4714-9045-a9669005d4c4",
   "metadata": {},
   "outputs": [],
   "source": [
    "sq = shap_quant(dc_shap_obj,dc_shap_ival_obj, X_train, feature_name_dict, 'tenure', 'EWA_avg_order_time_8')\n",
    "sq.is_linear()\n",
    "sq.isv_ifv_correlation()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e5c98ff-f227-4c29-bc78-0d274a2815bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "sq = shap_quant(dc_shap_obj,dc_shap_ival_obj, X_train, feature_name_dict, 'tenure', 'EWA_avg_note_time_8')\n",
    "sq.is_linear()\n",
    "sq.isv_ifv_correlation()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "66dfeaf9-6c0d-4e6f-beb7-e59b08cbdb9a",
   "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}