evaluation/Classification_Analysis_sota.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calibrated Explanations for Binary Classification\n",
    "## Stability and Robustness\n",
    "\n",
    "Author: Tuwe Löfström (tuwe.lofstrom@ju.se)  \n",
    "Copyright 2023 Tuwe Löfström  \n",
    "License: BSD 3 clause\n",
    "Sources:\n",
    "1. [\"Calibrated Explanations: with Uncertainty Information and Counterfactuals\"](https://arxiv.org/abt/2305.02305) by [Helena Löfström](https://github.com/Moffran), [Tuwe Löfström](https://github.com/tuvelofstrom), Ulf Johansson, and Cecilia Sönströd.\n",
    "\n",
    "### 1. Import packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle\n",
    "import numpy as np\n",
    "from scipy import stats as st"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2 Import results from the pickled result file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('results_sota.pkl', 'rb') as f:\n",
    "    results = pickle.load(f)\n",
    "data_characteristics = {'colic': 60, \n",
    "                        'creditA': 43, \n",
    "                        'diabetes': 9, \n",
    "                        'german': 28, \n",
    "                        'haberman': 4, \n",
    "                        'haberman': 4,\n",
    "                        'heartC': 23,\n",
    "                        'heartH': 21,\n",
    "                        'heartS': 14,\n",
    "                        'hepati': 20,\n",
    "                        'iono': 34,\n",
    "                        'je4042': 9,\n",
    "                        'je4243': 9, \n",
    "                        'kc1': 22,\n",
    "                        'kc2': 22,\n",
    "                        'kc3': 40,\n",
    "                        'liver': 7,\n",
    "                        'pc1req': 9,\n",
    "                        'pc4': 38,\n",
    "                        'sonar': 61,\n",
    "                        'spect': 23,\n",
    "                        'spectf': 45,\n",
    "                        'transfusion': 5,\n",
    "                        'ttt': 28,\n",
    "                        'vote': 17,\n",
    "                        'wbc': 10,}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3 Stability and Robustness\n",
    "In order to verify some of the claims of CE, two experiments have been run. The two experiments evaluate the stability and the robustness of the method. Stability is evaluated through experiments where the same model, calibration set and test set have been explained 30 times per data set. The only source of variation was the random seed. Robustness is evaluated through experiments where training and calibration sets have been randomly resampled before a new model was trained and explained. The experiment where run 30 times per data set with the same test set for all runs. A test set with 20 stratified instances (making sure both classes being equally represented) were used. Robustness is measured in this way to avoid inferring perturbed instances which are not from the same distribution as the test instances being explained. The expectation is that a stable and robust explanation method should result in low variance in the feature weights. \n",
    "\n",
    "Both random forests and xGBoost are used and both factual and counterfactual explanations are evaluated. The probability estimate of each of the models was computed on the same test set, as comparison to the robustness results. Furthermore, the probability estimates from each model calibrated using VA was also computed on the same test set for comparison.\n",
    "\n",
    "The two state-of-the-art (sota) techniques LIME (version 0.2.0.1 and using `LimeTabularExplainer` class) and SHAP (version 0.44.0 `Explainer` class) were also evaluated in the same way and using the same instances and models. These two techniques were selected as the two obvious sota techniques based on their accessibility (through e.g., `pip` installation) and their large user base. No obvious sota technique for counterfactuals were identified, as all the proposed algorithms, like e.g., LORE or MACE, seem to lack in either accessibility, user base, or both.\n",
    "\n",
    "Everything was run on 25 datasets. See the `Classification_Experiment_sota.py` for details on the experiment.\n",
    "\n",
    "The tabulated results are the mean variance of the stability and robustness measured over the 30 runs and 20 instances. The variance is measured per instance and computed over the 30 runs on the feature importance weight of the most influential feature, defined as the feature most often having highest absolute feature importance weight. The average variance is computed over the 20 instances. The most influential feature is used since it is the feature that is most likely to be used in a decision but also the feature with the greatest expected variation (as a consequence of the weights having the highest absolute values). \n",
    "\n",
    "The results are printed as a latex table for inclusion in the paper."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Stability & xGB & xGB & xGB & xGB & xGB & xGB \\\\\n",
      "Dataset & CE & CCE & L_C & S_C & L_U & S_U \\\\\n",
      "\\hline\n",
      "haberman & 4e-34 & 5e-34 & 1e-63 & 0e+00 & 1e-04 & 3e-33 & 1e-33 & 9e-34 & 9e-64 & 0e+00 & 4e-05 & 1e-33 \\\\\n",
      "heartC & 8e-34 & 0e+00 & 2e-10 & 1e-32 & 7e-05 & 1e-04 & 4e-33 & 3e-33 & 2e-10 & 5e-33 & 1e-05 & 5e-06 \\\\\n",
      "heartH & 6e-37 & 6e-37 & 2e-10 & 0e+00 & 5e-05 & 9e-05 & 5e-33 & 5e-33 & 2e-10 & 2e-32 & 8e-06 & 8e-06 \\\\\n",
      "heartS & 2e-34 & 2e-33 & 7e-11 & 6e-33 & 6e-05 & 5e-05 & 4e-33 & 5e-33 & 7e-11 & 1e-32 & 1e-05 & 3e-06 \\\\\n",
      "hepati & 6e-34 & 0e+00 & 7e-05 & 5e-33 & 3e-05 & 6e-05 & 4e-33 & 4e-33 & 8e-05 & 0e+00 & 1e-05 & 2e-06 \\\\\n",
      "je4243 & 7e-34 & 6e-34 & 5e-06 & 2e-34 & 9e-05 & 2e-33 & 2e-33 & 3e-33 & 1e-05 & 2e-34 & 1e-05 & 8e-34 \\\\\n",
      "kc2 & 7e-34 & 1e-33 & 4e-64 & 0e+00 & 7e-05 & 2e-04 & 1e-33 & 2e-33 & 1e-63 & 0e+00 & 5e-06 & 2e-05 \\\\\n",
      "kc3 & 4e-35 & 4e-35 & 4e-06 & 0e+00 & 3e-05 & 9e-05 & 4e-34 & 6e-34 & 4e-06 & 0e+00 & 1e-05 & 5e-06 \\\\\n",
      "pc1req & 2e-34 & 8e-35 & 1e-11 & 3e-33 & 8e-05 & 5e-33 & 2e-33 & 1e-33 & 2e-11 & 2e-33 & 4e-05 & 2e-33 \\\\\n",
      "spect & 0e+00 & 0e+00 & 3e-12 & 1e-32 & 1e-05 & 3e-05 & 6e-34 & 6e-34 & 4e-12 & 1e-32 & 3e-06 & 7e-06 \\\\\n",
      "transfusion & 9e-34 & 5e-34 & 7e-64 & 0e+00 & 7e-05 & 3e-33 & 2e-33 & 4e-33 & 7e-64 & 0e+00 & 3e-05 & 1e-33 \\\\\n",
      "vote & 0e+00 & 0e+00 & 3e-11 & 4e-33 & 6e-05 & 8e-05 & 8e-33 & 8e-33 & 3e-11 & 0e+00 & 1e-05 & 2e-05 \\\\\n",
      "wbc & 2e-33 & 2e-34 & 3e-05 & 0e+00 & 1e-04 & 2e-33 & 1e-33 & 5e-33 & 3e-05 & 0e+00 & 4e-05 & 8e-34 \\\\\n",
      "\\hline\n",
      "Average & 5e-34 & 4e-34 & 8e-06 & 3e-33 & 6e-05 & 6e-05 & 3e-33 & 3e-33 & 1e-05 & 4e-33 & 2e-05 & 5e-06 \\\\\n"
     ]
    }
   ],
   "source": [
    "stab_rank = {}\n",
    "stab_val = {}\n",
    "average_results = {}\n",
    "for a in ['xGB','RF']:\n",
    "    average_results[a+'_stab_ce'] = {}\n",
    "    average_results[a+'_stab_cce'] = {}\n",
    "    average_results[a+'_stab_lime'] = {}\n",
    "    average_results[a+'_stab_lime_va'] = {}\n",
    "    average_results[a+'_stab_shap'] = {}\n",
    "    average_results[a+'_stab_shap_va'] = {}\n",
    "\n",
    "n = results['test_size']\n",
    "r = results['num_rep']\n",
    "print('Stability & xGB & xGB & xGB & xGB & xGB & xGB \\\\\\\\')\n",
    "print('Dataset & CE & CCE & L_C & S_C & L_U & S_U \\\\\\\\\\n\\\\hline')\n",
    "for d in np.sort([k for k in results.keys()]):\n",
    "    if d in ['test_size', 'num_rep']:\n",
    "        continue\n",
    "    print(d, end='')\n",
    "    algorithms = results[d].keys()\n",
    "    for a in algorithms:\n",
    "        stability = results[d][a]['stability']\n",
    "        \n",
    "        for key in ['ce', 'cce']:    \n",
    "            ranks = []\n",
    "            for j in range(n):\n",
    "                rank = []\n",
    "                for i in range(r):\n",
    "                    rank.append(np.argsort(np.abs(stability[key][i][j]['predict']))[-1:][0])\n",
    "                ranks.append(rank)\n",
    "            stab_rank[key] = st.mode(ranks, axis=1)[0] # Find most important feature per instance\n",
    "            value = []\n",
    "            for j in range(n):\n",
    "                value.append([np.mean([stability[key][i][j]['predict'][stab_rank[key][j]] for i in range(r)]), np.var([stability[key][i][j]['predict'][stab_rank[key][j]] for i in range(r)])])\n",
    "            stab_val[key] = value \n",
    "\n",
    "        stability = results[d][a]['stability']\n",
    "                    \n",
    "        for key in ['lime', 'lime_va', 'shap', 'shap_va']:    \n",
    "            ranks = []\n",
    "            for j in range(n):\n",
    "                rank = []\n",
    "                for i in range(r):\n",
    "                    rank.append(np.argsort(np.abs(stability[key][i][j]))[-1:][0])\n",
    "                ranks.append(rank)\n",
    "            stab_rank[key] = st.mode(ranks, axis=1)[0] # Find most important feature per instance\n",
    "            value = []\n",
    "            for j in range(n):\n",
    "                value.append([np.mean([stability[key][i][j][stab_rank[key][j]] for i in range(r)]), np.var([stability[key][i][j][stab_rank[key][j]] for i in range(r)])])\n",
    "            stab_val[key] = value \n",
    "        \n",
    "        for key in ['ce', 'cce', 'lime', 'lime_va', 'shap', 'shap_va']:\n",
    "            average_results[a+'_stab_'+key][d] = np.mean([t[1] for t in stab_val[key]])\n",
    "        # print(f'{np.mean([t[1] if t[1] > 1e-20 else 0 for t in stab_val[\"ce\"]]):.1e} & {np.var([t[1] if t[1] > 1e-20 else 0 for t in stab_val[\"cce\"]]):.1e} & ',end='')\n",
    "        \n",
    "        for key in ['ce', 'cce', 'lime_va', 'shap_va', 'lime', 'shap']:\n",
    "            print(f' & {average_results[a+\"_stab_\"+key][d]:.0e}',end='')\n",
    "    print(' \\\\\\\\')\n",
    "print('\\\\hline\\nAverage', end='')\n",
    "for a in algorithms:\n",
    "    for key in ['ce', 'cce', 'lime_va', 'shap_va', 'lime', 'shap']:\n",
    "        print(f' & {np.mean([v for v in average_results[a+\"_stab_\"+key].values()]):.0e}',end='')\n",
    "print(' \\\\\\\\')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As can be seen above, the stability is practically 0 for both factual CE (CE) and counterfactual CE (CCE), illustrating that the method is stable by definition. Explanations extracted using SHAP (S_C) from calibrated models are also practically 0. LIME on calibrated models (L_C) and both LIME (L_U) and SHAP (S_U) on uncalibrated models are clearly less stable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Robustness & xGB & xGB & xGB & xGB & xGB & xGB & xGB & xGB \\\\\n",
      "Dataset & CE & CCE & L_C & S_C & C & L_U & S_U & U \\\\\n",
      "\\hline\n",
      "haberman & 0.011 & 0.009 & 0.000 & 0.000 & 0.010 & 0.008 & 0.011 & 0.033 & 0.011 & 0.009 & 0.000 & 0.000 & 0.009 & 0.004 & 0.004 & 0.024 \\\\\n",
      "heartC & 0.014 & 0.014 & 0.001 & 0.001 & 0.006 & 0.006 & 0.004 & 0.025 & 0.012 & 0.011 & 0.001 & 0.001 & 0.005 & 0.001 & 0.001 & 0.020 \\\\\n",
      "heartH & 0.016 & 0.015 & 0.171 & 0.002 & 0.004 & 0.006 & 0.004 & 0.022 & 0.011 & 0.011 & 0.186 & 0.001 & 0.003 & 0.001 & 0.001 & 0.019 \\\\\n",
      "heartS & 0.020 & 0.020 & 0.001 & 0.001 & 0.003 & 0.005 & 0.003 & 0.021 & 0.017 & 0.015 & 0.002 & 0.001 & 0.003 & 0.002 & 0.001 & 0.020 \\\\\n",
      "hepati & 0.021 & 0.018 & 0.000 & 0.001 & 0.003 & 0.011 & 0.006 & 0.021 & 0.014 & 0.014 & 0.000 & 0.001 & 0.003 & 0.003 & 0.001 & 0.020 \\\\\n",
      "je4243 & 0.011 & 0.009 & 0.003 & 0.000 & 0.003 & 0.005 & 0.008 & 0.021 & 0.010 & 0.007 & 0.002 & 0.000 & 0.003 & 0.001 & 0.001 & 0.020 \\\\\n",
      "kc2 & 0.019 & 0.010 & 0.000 & 0.000 & 0.003 & 0.011 & 0.009 & 0.022 & 0.007 & 0.005 & 0.000 & 0.000 & 0.004 & 0.000 & 0.001 & 0.021 \\\\\n",
      "kc3 & 0.006 & 0.005 & 0.013 & 0.000 & 0.004 & 0.003 & 0.001 & 0.020 & 0.005 & 0.005 & 0.012 & 0.000 & 0.005 & 0.001 & 0.000 & 0.019 \\\\\n",
      "pc1req & 0.013 & 0.012 & 0.009 & 0.004 & 0.005 & 0.017 & 0.011 & 0.020 & 0.014 & 0.013 & 0.009 & 0.004 & 0.005 & 0.004 & 0.004 & 0.020 \\\\\n",
      "spect & 0.007 & 0.007 & 0.003 & 0.001 & 0.005 & 0.003 & 0.002 & 0.019 & 0.006 & 0.006 & 0.004 & 0.002 & 0.005 & 0.001 & 0.000 & 0.018 \\\\\n",
      "transfusion & 0.008 & 0.006 & 0.000 & 0.000 & 0.005 & 0.004 & 0.008 & 0.019 & 0.008 & 0.006 & 0.000 & 0.000 & 0.005 & 0.003 & 0.003 & 0.019 \\\\\n",
      "vote & 0.014 & 0.014 & 0.001 & 0.001 & 0.005 & 0.002 & 0.002 & 0.019 & 0.011 & 0.011 & 0.001 & 0.001 & 0.005 & 0.001 & 0.001 & 0.019 \\\\\n",
      "wbc & 0.020 & 0.012 & 0.001 & 0.000 & 0.004 & 0.008 & 0.003 & 0.018 & 0.018 & 0.009 & 0.001 & 0.000 & 0.004 & 0.003 & 0.001 & 0.018 \\\\\n",
      "\\hline\n",
      "Average & 0.014 & 0.012 & 0.016 & 0.001 & 0.005 & 0.007 & 0.005 & 0.022 & 0.011 & 0.010 & 0.017 & 0.001 & 0.004 & 0.002 & 0.001 & 0.020 \\\\\n"
     ]
    }
   ],
   "source": [
    "rob_rank = {}\n",
    "rob_val = {}\n",
    "rob_proba = []\n",
    "rob_proba_va = []\n",
    "# average_results = {}\n",
    "for a in ['xGB','RF']:\n",
    "    average_results[a+'_rob_ce'] = {}\n",
    "    average_results[a+'_rob_cce'] = {}\n",
    "    average_results[a+'_rob_lime'] = {}\n",
    "    average_results[a+'_rob_lime_va'] = {}\n",
    "    average_results[a+'_rob_shap'] = {}\n",
    "    average_results[a+'_rob_shap_va'] = {}\n",
    "    average_results[a+'_rob_proba'] = {}\n",
    "    average_results[a+'_rob_proba_va'] = {}\n",
    "\n",
    "n = results['test_size']\n",
    "r = results['num_rep']\n",
    "\n",
    "print('Robustness & xGB & xGB & xGB & xGB & xGB & xGB & xGB & xGB \\\\\\\\')\n",
    "print('Dataset & CE & CCE & L_C & S_C & C & L_U & S_U & U \\\\\\\\\\n\\\\hline')\n",
    "for d in np.sort([k for k in results.keys()]):\n",
    "    if d in ['test_size', 'num_rep']:\n",
    "        continue\n",
    "    print(d, end='')\n",
    "    algorithms = results[d].keys()\n",
    "    for a in algorithms:\n",
    "        robustness = results[d][a]['robustness']\n",
    "        \n",
    "        for key in ['ce', 'cce']:               \n",
    "            ranks = []\n",
    "            values = []\n",
    "            for j in range(n):\n",
    "                rank = []\n",
    "                value = []\n",
    "                for i in range(r):\n",
    "                    rank.append(np.argsort(np.abs(robustness[key][i][j]['predict']))[-1:][0])\n",
    "                ranks.append(rank)\n",
    "                values.append(value)\n",
    "            rob_rank[key] = st.mode(ranks, axis=1)[0] # Find most important feature per instance\n",
    "            value = []\n",
    "            for j in range(n):\n",
    "                value.append([np.mean([robustness[key][i][j]['predict'][rob_rank[key][j]] for i in range(r)]), np.var([robustness[key][i][j]['predict'][rob_rank[key][j]] for i in range(r)])])\n",
    "            rob_val[key] = value\n",
    "\n",
    "        robustness = results[d][a]['robustness']\n",
    "            \n",
    "        for key in ['lime', 'lime_va', 'shap', 'shap_va']:    \n",
    "            ranks = []\n",
    "            for j in range(n):\n",
    "                rank = []\n",
    "                for i in range(r):\n",
    "                    rank.append(np.argsort(np.abs(robustness[key][i][j]))[-1:][0])\n",
    "                ranks.append(rank)\n",
    "            rob_rank[key] = st.mode(ranks, axis=1)[0] # Find most important feature per instance\n",
    "            value = []\n",
    "            for j in range(n):\n",
    "                value.append([np.mean([robustness[key][i][j][rob_rank[key][j]] for i in range(r)]), np.var([robustness[key][i][j][rob_rank[key][j]] for i in range(r)])])\n",
    "            rob_val[key] = value \n",
    "        \n",
    "        for inst in range(n):\n",
    "            rob_proba.append(np.var([robustness['proba'][j][inst] for j in range(r)]))\n",
    "            rob_proba_va.append(np.var([robustness['proba_va'][j][inst] for j in range(r)]))\n",
    "        \n",
    "        for key in ['ce', 'cce', 'lime', 'lime_va', 'shap', 'shap_va']:\n",
    "            average_results[a+'_rob_'+key][d] = np.mean([t[1] for t in rob_val[key]])\n",
    "        average_results[a+'_rob_proba'][d] = np.mean(rob_proba)\n",
    "        average_results[a+'_rob_proba_va'][d] = np.mean(rob_proba_va)\n",
    "        # print(f'{np.mean([t[1] if t[1] > 1e-20 else 0 for t in rob_val[\"ce\"]]):.1e} & {np.mean([t[1] if t[1] > 1e-20 else 0 for t in rob_val[\"cce\"]]):.1e} & ',end='')\n",
    "        \n",
    "        for key in ['ce', 'cce', 'lime_va', 'shap_va', 'proba_va', 'lime', 'shap', 'proba']:\n",
    "            print(f' & {average_results[a+\"_rob_\"+key][d]:.3f}',end='')\n",
    "    print(' \\\\\\\\')\n",
    "print('\\\\hline\\nAverage', end='')\n",
    "for a in algorithms:\n",
    "    for key in ['ce', 'cce', 'lime_va', 'shap_va', 'proba_va', 'lime', 'shap', 'proba']:\n",
    "        print(f' & {np.mean([v for v in average_results[a+\"_rob_\"+key].values()]):.3f}',end='')\n",
    "print(' \\\\\\\\')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, the picture is somewhat different, as CE and CCE have higher variability in their feature weights compared to LIME and SHAP. The results for CE and CCE are clearly lower than the uncalibrated model (U) from which they have been extracted and slightly more than the VA calibrated model (C). While these results could be interpreted as indicating low robustness, we argue that the experiment shows that the method updates its feature weights in accordance with how much the underlying model is changing its predictions. A similar pattern could be seen when comparing LIME and SHAP explanations from calibrated vs uncalibrated models, with low variability for the explanations of calibrated models having themselves lower variability and vice versa."
   ]
  }
 ],
 "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.11.5"
  },
  "vscode": {
   "interpreter": {
    "hash": "168dd73c7a7b76a0355e35f33a90e68c167b1dbb1e524891be00dd5c7b8524eb"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}