doc/source/examples/cfproto_housing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Counterfactuals guided by prototypes on Boston housing dataset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook goes through an example of [prototypical counterfactuals](../methods/CFProto.ipynb) using [k-d trees](https://en.wikipedia.org/wiki/K-d_tree) to build the prototypes. Please check out [this notebook](./cfproto_mnist.ipynb) for a more in-depth application of the method on MNIST using (auto-)encoders and trust scores.\n",
    "\n",
    "In this example, we will train a simple neural net to predict whether house prices in the Boston area are above the median value or not. We can then find a counterfactual to see which variables need to be changed to increase or decrease a house price above or below the median value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TF version:  2.2.0\n",
      "Eager execution enabled:  False\n"
     ]
    }
   ],
   "source": [
    "import tensorflow as tf\n",
    "tf.get_logger().setLevel(40) # suppress deprecation messages\n",
    "tf.compat.v1.disable_v2_behavior() # disable TF2 behaviour as alibi code still relies on TF1 constructs \n",
    "from tensorflow.keras.layers import Dense, Input\n",
    "from tensorflow.keras.models import Model, load_model\n",
    "from tensorflow.keras.utils import to_categorical\n",
    "import matplotlib\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import os\n",
    "from sklearn.datasets import load_boston\n",
    "from alibi.explainers import CounterfactualProto\n",
    "\n",
    "print('TF version: ', tf.__version__)\n",
    "print('Eager execution enabled: ', tf.executing_eagerly()) # False"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load and prepare Boston housing dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "boston = load_boston()\n",
    "data = boston.data\n",
    "target = boston.target\n",
    "feature_names = boston.feature_names"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Transform into classification task: target becomes whether house price is above the overall median or not"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = np.zeros((target.shape[0],))\n",
    "y[np.where(target > np.median(target))[0]] = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remove categorical feature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = np.delete(data, 3, 1)\n",
    "feature_names = np.delete(feature_names, 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Explanation of remaining features:\n",
    "\n",
    "- CRIM: per capita crime rate by town\n",
    "- ZN: proportion of residential land zoned for lots over 25,000 sq.ft.\n",
    "- INDUS: proportion of non-retail business acres per town\n",
    "- RM: average number of rooms per dwelling\n",
    "- AGE: proportion of owner-occupied units built prior to 1940\n",
    "- DIS: weighted distances to five Boston employment centres\n",
    "- RAD: index of accessibility to radial highways\n",
    "- TAX: full-value property-tax rate per USD10,000\n",
    "- PTRATIO: pupil-teacher ratio by town\n",
    "- B: 1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town\n",
    "- LSTAT: % lower status of the population"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Standardize data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu = data.mean(axis=0)\n",
    "sigma = data.std(axis=0)\n",
    "data = (data - mu) / sigma"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define train and test set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "idx = 475\n",
    "x_train,y_train = data[:idx,:], y[:idx]\n",
    "x_test, y_test = data[idx:,:], y[idx:]\n",
    "y_train = to_categorical(y_train)\n",
    "y_test = to_categorical(y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Train model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(42)\n",
    "tf.random.set_seed(42)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def nn_model():\n",
    "    x_in = Input(shape=(12,))\n",
    "    x = Dense(40, activation='relu')(x_in)\n",
    "    x = Dense(40, activation='relu')(x)\n",
    "    x_out = Dense(2, activation='softmax')(x)\n",
    "    nn = Model(inputs=x_in, outputs=x_out)\n",
    "    nn.compile(loss='categorical_crossentropy', optimizer='sgd', metrics=['accuracy'])\n",
    "    return nn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Model: \"model\"\n",
      "_________________________________________________________________\n",
      "Layer (type)                 Output Shape              Param #   \n",
      "=================================================================\n",
      "input_1 (InputLayer)         [(None, 12)]              0         \n",
      "_________________________________________________________________\n",
      "dense (Dense)                (None, 40)                520       \n",
      "_________________________________________________________________\n",
      "dense_1 (Dense)              (None, 40)                1640      \n",
      "_________________________________________________________________\n",
      "dense_2 (Dense)              (None, 2)                 82        \n",
      "=================================================================\n",
      "Total params: 2,242\n",
      "Trainable params: 2,242\n",
      "Non-trainable params: 0\n",
      "_________________________________________________________________\n"
     ]
    }
   ],
   "source": [
    "nn = nn_model()\n",
    "nn.summary()\n",
    "nn.fit(x_train, y_train, batch_size=64, epochs=500, verbose=0)\n",
    "nn.save('nn_boston.h5', save_format='h5')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test accuracy:  0.83870965\n"
     ]
    }
   ],
   "source": [
    "nn = load_model('nn_boston.h5')\n",
    "score = nn.evaluate(x_test, y_test, verbose=0)\n",
    "print('Test accuracy: ', score[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate counterfactual guided by the nearest class prototype"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Original instance:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = x_test[1].reshape((1,) + x_test[1].shape)\n",
    "shape = X.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run counterfactual:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# define model\n",
    "nn = load_model('nn_boston.h5')\n",
    "\n",
    "# initialize explainer, fit and generate counterfactual\n",
    "cf = CounterfactualProto(nn, shape, use_kdtree=True, theta=10., max_iterations=1000,\n",
    "                         feature_range=(x_train.min(axis=0), x_train.max(axis=0)), \n",
    "                         c_init=1., c_steps=10)\n",
    "\n",
    "cf.fit(x_train)\n",
    "explanation = cf.explain(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The prediction flipped from 0 (value below the median) to 1 (above the median):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Original prediction: 0\n",
      "Counterfactual prediction: 1\n"
     ]
    }
   ],
   "source": [
    "print(f'Original prediction: {explanation.orig_class}')\n",
    "print('Counterfactual prediction: {}'.format(explanation.cf['class']))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a look at the counterfactual. To make the results more interpretable, we will first undo the pre-processing step and then check where the counterfactual differs from the original instance:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "AGE: -6.526239960695747\n",
      "LSTAT: -4.799340540220259\n"
     ]
    }
   ],
   "source": [
    "orig = X * sigma + mu\n",
    "counterfactual = explanation.cf['X'] * sigma + mu\n",
    "delta = counterfactual - orig\n",
    "for i, f in enumerate(feature_names):\n",
    "    if np.abs(delta[0][i]) > 1e-4:\n",
    "        print('{}: {}'.format(f, delta[0][i]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So in order to increase the house price, the proportion of owner-occupied units built prior to 1940 should decrease by ~11-12%. This is not surprising since the proportion for the observation is very high at 93.6%. Furthermore, the % of the population with \"lower status\" should decrease by ~5%."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "% owner-occupied units built prior to 1940: 93.6\n",
      "% lower status of the population: 18.68\n"
     ]
    }
   ],
   "source": [
    "print('% owner-occupied units built prior to 1940: {}'.format(orig[0][5]))\n",
    "print('% lower status of the population: {}'.format(orig[0][11]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clean up:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.remove('nn_boston.h5')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}