#!/usr/bin/env python
# coding: utf-8

# # <center>Model Interpretability on Random Forest using LIME</center>

# ## Table of Contents
# 
# 1. [Problem Statement](#section1)<br><br>
# 2. [Importing Packages](#section2)<br><br>
# 3. [Loading Data](#section3)
#   - 3.1 [Description of the Dataset](#section301)<br><br>
# 4. [Data Preprocessing](#section4)<br><br>
# 5. [Data train/test split](#section5)<br><br>
# 6. [Random Forest Model](#section6)
#   - 6.1 [Random Forest in scikit-learn](#section601)<br>  <br> 
#   - 6.2 [Using the Model for Prediction](#section602)<br><br>
# 7. [Model Evaluation](#section7)
#    - 7.1 [Accuracy Score](#section701)<br><br>
# 8. [Model Interpretability using LIME](#section8)  
#   - 8.1 [Setup LIME Algorithm](#section801)<br><br>
#   - 8.2 [Explore Key Features in Instance-by-Instance Predictions](#section802)<br>

# <a id=section1></a>
# ## 1. Problem Statement
# 
# - We have often found that **Machine Learning (ML)** algorithms capable of capturing **structural non-linearities** in training data - models that are sometimes referred to as **'black box' (e.g. Random Forests, Deep Neural Networks, etc.)** - perform far **better at prediction** than their **linear counterparts (e.g. Generalized Linear Models)**. 
# 
# 
# - They are, however, much **harder to interpret** - in fact, quite often it is **not possible to gain any insight into why a particular prediction has been produced**, when given an **instance of input data (i.e. the model features)**. 
# 
# 
# - Consequently, it has **not been possible to use 'black box' ML algorithms** in situations where clients have sought **cause-and-effect explanations for model predictions**, with end-results being that sub-optimal predictive models have been used in their place, as their explanatory power has been more valuable, in relative terms.
# 
# 
# - The **problem with model explainability** is that it’s **very hard to define a model’s decision boundary in human understandable manner**. 
# 
# 
# - **LIME** is a **python library** which tries to **solve for model interpretability by producing locally faithful explanations**. 
# 
# <br>
# <img src="../../images/lime.jpg" width="700" height="550"/> <br><br>
# 
# 
# - We will use **LIME** to **interpret** our **RandomForest model**.

# ---

# <a id=section2></a>
# ## 2. Importing Packages

# In[ ]:


# Install LIME using the following command.

get_ipython().system('pip install lime')


# In[1]:


import numpy as np
np.set_printoptions(precision=4)                    # To display values only upto four decimal places. 

import pandas as pd
pd.set_option('mode.chained_assignment', None)      # To suppress pandas warnings.
pd.set_option('display.max_colwidth', -1)           # To display all the data in the columns.
pd.options.display.max_columns = 40                 # To display all the columns. (Set the value to a high number)

import matplotlib.pyplot as plt
plt.style.use('seaborn-whitegrid')                  # To apply seaborn whitegrid style to the plots.
plt.rc('figure', figsize=(10, 8))                   # Set the default figure size of plots.
get_ipython().run_line_magic('matplotlib', 'inline')

import warnings
warnings.filterwarnings('ignore')                   # To suppress all the warnings in the notebook.

from sklearn.preprocessing import LabelEncoder
from sklearn.preprocessing import OneHotEncoder
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import accuracy_score


# ---

# <a id=section3></a>
# ## 3. Loading Data

# In[2]:


df = pd.read_csv('../../data/mushrooms.csv')
df.head()


# <a id=section301></a>
# ### 3.1 Description of the Dataset

# - This **dataset includes descriptions** of hypothetical samples corresponding to **23 species** of **gilled mushrooms** in the **Agaricus and Lepiota Family Mushroom** drawn from **The Audubon Society Field Guid**e to **North American Mushrooms (1981)**. 
# 
# 
# - **Each species** is **identified as definitely edible**, **definitely poisonous**, or **of unknown edibility** and **not recommended**. This **latter class was combined with** the **poisonous one**. 
# 
# 
# - The **Guide clearly states** that there is no **simple rule for determining** the **edibility of a mushroom**; no rule like **"leaflets three, let it be''** for **Poisonous Oak and Ivy**.

# In[3]:


df.columns


# | **Column Name**                  | **Description**                                                                          |
# | ---------------------------------|:----------------------------------------------------------------------------------------:| 
# | class                            | classes: edible=e, poisonous=p.                                                          |
# | cap-shape                        | bell=b,conical=c, convex=x, flat=f, knobbed=k, sunken=s.                                 |
# | cap-surface                      | fibrous=f, grooves=g, scaly=y, smooth=s.                                                 |
# | cap-color                        | brown=n, buff=b, cinnamon=c, gray=g, green=r, pink=p, purple=u, red=e, white=w, yellow=y.|
# | bruises                          | bruises=t, no=f.                                                                         |
# | odor                             | almond=a, anise=l, creosote=c, fishy=y, foul=f, musty=m ,none=n, pungent=p, spicy=s.     |
# | gill-attachment                  | attached=a, descending=d, free=f, notched=n.                                             |
# | gill-spacing                     | close=c, crowded=w, distant=d.                                                           |
# | gill-size                        | broad=b, narrow=n.                                                                       |
# | gill-color                       | black=k, brown=n ,buff=b, chocolate=h, gray=g, green=r, orange=o, pink=p, purple=u, red=e, white=w, yellow=y.                                                                                                            |
# | stalk-shape                      | enlarging=e, tapering=t.                                                                 |
# | stalk-root                       | bulbous=b, club=c, cup=u, equal=e, rhizomorphs=z, rooted=r, missing=?.                   |
# | stalk-surface-above-ring         | fibrous=f, scaly=y, silky=k, smooth=s.                                                   |
# | stalk-surface-below-ring         | fibrous=f, scaly=y, silky=k, smooth=s.                                                   |
# | stalk-color-above-ring           | brown=n, buff=b, cinnamon=c, gray=g, orange=o, pink=p, red=e, white=w, yellow=y.         |
# | stalk-color-below-ring           | brown=n, buff=b, cinnamon=c, gray=g, orange=o, pink=p, red=e, white=w, yellow=y.         |
# | veil-type                        | partial=p ,universal=u.                                                                  |
# | veil-color                       | brown=n, orange=o, white=w, yellow=y.                                                    |
# | ring-number                      | none=n, one=o, two=t.                                                                    |
# | ring-type                        | cobwebby=c, evanescent=e, flaring=f, large=l, none=n, pendant=p, sheathing=s, zone=z.    |
# | spore-print-color                | black=k, brown=n, buff=b, chocolate=h, green=r, orange=o, purple=u, white=w, yellow=y.   |
# | population                       | abundant=a, clustered=c, numerous=n, scattered=s, several=v, solitary=y.                 |
# | habitat                          | grasses=g, leaves=l, meadows=m, paths=p, urban=u, waste=w, woods=d.                      |

# In[4]:


df.info()


# In[5]:


df.describe()


# ---

# <a id=section4></a>
# ## 4 Data Preprocessing

# In[36]:


df.head()


# In[37]:


# Creating labels array from the class column.

labels = df.iloc[:, 0].values
labels


# In[38]:


# Creating a LabelEncoder object le and fitting labels array into it.

le = LabelEncoder()
le.fit(labels)


# In[39]:


# Transforming the labels array to have numerical values.

labels = le.transform(labels)
labels


# In[40]:


# Storing the different classes found by LabelEncoder in labels array into class_names.

class_names = le.classes_
class_names


# In[41]:


# Dropping the class column from the df dataframe.

df.drop(['class'], axis=1, inplace=True)


# In[42]:


df.head()


# In[43]:


# Creating a range form 0 upto the number of categorical features. Since all the features in df are categorical using len(). 

categorical_features = range(len(df.columns))
categorical_features


# In[44]:


# Creating an array of feature names.

feature_names = df.columns.values
feature_names


# In[45]:


# We expand the characters into words, using the dataset description provided in the beginning.

categorical_names = '''bell=b,conical=c,convex=x,flat=f,knobbed=k,sunken=s
fibrous=f,grooves=g,scaly=y,smooth=s
brown=n,buff=b,cinnamon=c,gray=g,green=r,pink=p,purple=u,red=e,white=w,yellow=y
bruises=t,no=f
almond=a,anise=l,creosote=c,fishy=y,foul=f,musty=m,none=n,pungent=p,spicy=s
attached=a,descending=d,free=f,notched=n
close=c,crowded=w,distant=d
broad=b,narrow=n
black=k,brown=n,buff=b,chocolate=h,gray=g,green=r,orange=o,pink=p,purple=u,red=e,white=w,yellow=y
enlarging=e,tapering=t
bulbous=b,club=c,cup=u,equal=e,rhizomorphs=z,rooted=r,missing=?
fibrous=f,scaly=y,silky=k,smooth=s
fibrous=f,scaly=y,silky=k,smooth=s
brown=n,buff=b,cinnamon=c,gray=g,orange=o,pink=p,red=e,white=w,yellow=y
brown=n,buff=b,cinnamon=c,gray=g,orange=o,pink=p,red=e,white=w,yellow=y
partial=p,universal=u
brown=n,orange=o,white=w,yellow=y
none=n,one=o,two=t
cobwebby=c,evanescent=e,flaring=f,large=l,none=n,pendant=p,sheathing=s,zone=z
black=k,brown=n,buff=b,chocolate=h,green=r,orange=o,purple=u,white=w,yellow=y
abundant=a,clustered=c,numerous=n,scattered=s,several=v,solitary=y
grasses=g,leaves=l,meadows=m,paths=p,urban=u,waste=w,woods=d'''.split('\n')

categorical_names[0]


# In[46]:


for j, names in enumerate(categorical_names):
    values = names.split(',')
    values = dict([(x.split('=')[1], x.split('=')[0]) for x in values])
    df.iloc[:, j] = df.iloc[:, j].map(values)


# In[47]:


df.head()


# In[48]:


# LabelEncoding all the features. Capturing the different class values for each feature in the categorical_names dictionary.

categorical_names = {}

for feature in categorical_features:
    le = LabelEncoder()
    le.fit(df.iloc[:, feature])
    df.iloc[:, feature] = le.transform(df.iloc[:, feature])
    categorical_names[feature] = le.classes_


# In[49]:


categorical_names[0]


# ---

# <a id=section5></a>
# ## 5. Data train/test split
# 
# - Now that the entire **data** is of **numeric datatype**, lets begin our modelling process.
# 
# 
# - Firstly, **splitting** the complete **dataset** into **training** and **testing** datasets.

# In[50]:


df.head()


# In[51]:


X = df.iloc[:, :]
X.head()


# In[52]:


y = labels[:]
y[:10]


# In[53]:


# Using scikit-learn's train_test_split function to split the dataset into train and test sets.
# 80% of the data will be in the train set and 20% in the test set, as specified by test_size=0.2

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)


# In[54]:


# Checking the shapes of all the training and test sets for the dependent and independent features.

print(X_train.shape)
print(y_train.shape)
print(X_test.shape)
print(y_test.shape)


# In[55]:


# Finally, we use a One-hot encoder, so that the classifier does not take our categorical features as continuous features. 
# We will use this encoder only for the classifier, not for the explainer - 
# and the reason is that the explainer must make sure that a categorical feature only has one value.

ohe = OneHotEncoder(categorical_features=categorical_features)
ohe.fit(df)


# In[56]:


X_train_encoded = ohe.transform(X_train)


# In[57]:


X_test_encoded = ohe.transform(X_test)


# In[58]:


print(X_train_encoded.shape)
print(X_test_encoded.shape)


# ---

# <a id=section6></a>
# ## 6.  Random Forest Model
# 
# <a id=section601></a>
# ### 6.1  Random Forest with Scikit-Learn

# In[74]:


# Creating a Random Forest Classifier.

classifier_rf = RandomForestClassifier(n_estimators=500, random_state=0, oob_score=True, n_jobs=-1)


# In[75]:


# Fitting the model on the dataset.

classifier_rf.fit(X_train_encoded, y_train)


# In[76]:


classifier_rf.oob_score_


# - From the **OOB score** we can see how our model's gonna perform against the **test set or new** samples.

# ---

# <a id=section602></a>
# ### 6.2 Using the Model for Prediction

# In[124]:


# Making predictions on the training set.

y_pred_train = classifier_rf.predict(X_train_encoded)
y_pred_train[:10]


# In[125]:


# Making predictions on test set.

y_pred_test = classifier_rf.predict(X_test_encoded)
y_pred_test[:10]


# ---

# <a id=section7></a>
# ## 7. Model Evaluation
# 
# **Error** is the deviation of the values predicted by the model with the true values.

# <a id=section701></a>
# ### 7.1 Accuracy Score

# In[126]:


# Accuracy score on the training set.

print('Accuracy score for train data is:', accuracy_score(y_train, y_pred_train))


# In[127]:


# Accuracy score on the test set.

print('Accuracy score for test data is:', accuracy_score(y_test, y_pred_test))


# - We get an **accuracy** of **100%** on our train set and an **accuracy** of **100%** on our test set.
# 
# 
# - We can notice that the **accuracy** obtained on the **test set (1.0)** is similar to the one obtained using the **oob_score_ (1.0)**, so we can use the **oob_score_** as a **validation** before testing our model on the **test set**. 

# ---

# <a id=section8></a>
# ## 8. Model Interpretability using LIME
# 
# 
# - **LIME** stands for **Local Interpretable Model-Agnostic Explanations** is a technique to **explain the predictions of any machine learning classifier**, and **evaluate its usefulness** in various **tasks** related to **trust**. 

# In[132]:


# Our predict function first transforms the data into the one-hot representation. 
# Then it calculates the prediction probability for each class of target variable.

predict_fn = lambda x: classifier_rf.predict_proba(ohe.transform(x))


# <a id=section801></a>
# ### 8.1 Setup LIME Algorithm

# - We now **create** our **explainer**. 
# 
# 
# - The **categorical_features parameter** lets it know **which features** are **categorica**l (in **this case**, **all of them**). 
# 
# 
# - The **categorical names parameter** gives a **string representation** of **each categorical feature's numerical value**.

# In[128]:


from lime.lime_tabular import LimeTabularExplainer


# In[129]:


# Creating the LIME explainer object.

explainer = LimeTabularExplainer(X_train.values, mode='classification', class_names=['edible', 'poisonous'], 
                                 feature_names = feature_names, categorical_features=categorical_features, 
                                 categorical_names=categorical_names, kernel_width=3, verbose=True, random_state=0)


# ---

# <a id=section802></a>
# ### 8.2  Explore Key Features in Instance-by-Instance Predictions

# - **Start by choosing an instance** from the **test dataset**.
# 
# 
# - Use **LIME** to **estimate a local model** to use for **explaining our model's predictions**. The **outputs** will be:
# 
#   1. The **intercept** estimated for the local model.
#   2. The **local model's estimate** for the **Regression Forest's prediction**.
#   3. The **Regression Forest's actual prediction**.
# 
# 
# - Note, that the **actual value from the data does not enter into this** - the **idea of LIME** is to **gain insight** into **why the chosen model** - in our case the Random Forest regressor - **is predicting whatever it has been asked to predict**. Whether or not this prediction is actually any good, is a separate issue.

# In[253]:


# Selecting a random instance from the test dataset.

i = np.random.randint(0, X_test.shape[0])
print('i =', i)


# In[254]:


# Using LIME to estimate a local model. Using only 6 features to explain our model's predictions.

exp = explainer.explain_instance(X_test.values[i], predict_fn, num_features=6)


# - **Printing** the **DataFrame row** for the **chosen test instance**.

# In[255]:


# Here the index column is the original index as per the df dataframe and the number at the beginning the index after reset.

X_test.reset_index().loc[[i]]


# -  **LIME's interpretation** of our **Random Forest's prediction**.

# In[256]:


exp.show_in_notebook(show_table=True, show_all=False)


# - First, note that the **row** we **explained** is **displayed** on the **right side**, in **table** format. Since we had the **show_all parameter** set to **false**, only the **features used in the explanation are displayed**.
# 
# 
# - The **value column** displays the **original value for each feature**.

# - To get the **output generated above** in the **form of a list**.

# In[257]:


exp.as_list()


# **Obesrvations obtained from LIME's interpretation of our Random Forest's prediction**:
# 
# - The **values** shown after the condition is the **amount** by which the value is **shifted** from the **intercept** estimated for the local model. 
#   
#   
# - When all these values are **added** to the **intercept**, it gives us the **Prediction_local** (local model's estimate for the Regression Forest's prediction) calculated by **LIME**.

# In[258]:


print('Intercept =', exp.intercept[1])
print('Prediction_local =', exp.local_pred[0])


# In[259]:


# Calculating the Prediction_local by adding all the values obtained above for each condition into the intercept.
# The intercept can be obtained from the exp.intercept using the index 0.

intercept = exp.intercept[1]
prediction_local = intercept

for j in range(len(exp.as_list())):
    prediction_local += exp.as_list()[j][1]

print('Prediction_local =', prediction_local)


#     
#     
#     

# - Choosing **another instance** from the **test dataset**.

# In[264]:


# This time specifying a particular value of i in order to explain the working of LIME.

i = 515
print('i =', i, '\n')

exp = explainer.explain_instance(X_test.values[i], predict_fn, num_features=6)


# - **Printing** the **DataFrame row** for the **chosen test instance**.

# In[277]:


X_test.reset_index().loc[[i]]


# -  **LIME's interpretation** of our **Random Forest's prediction**.

# In[266]:


exp.show_in_notebook(show_table=True, show_all=False)


# In[267]:


exp.as_list()


# In[268]:


print('Intercept =', exp.intercept[1])
print('Prediction_local =', exp.local_pred[0])


# In[269]:


intercept = exp.intercept[1]
prediction_local = intercept

for j in range(len(exp.as_list())):
    prediction_local += exp.as_list()[j][1]

print('Prediction_local =', prediction_local)


# - By **changing** the chosen **i**, we observe that the **narrative provided by LIME** also **changes, in response to changes in the model** in the **local region** of the **feature space** in which it is working to **generate a given prediction**. 
# 
# 
# - This is clearly an **improvement on relying purely** on the **Regression Forest's (static) expected relative feature importance** and of **great benefit to models that provice no insight whatsoever**.

# - Now note that the **explanations** are **based not only on features**, **but** on **feature-value pairs**. 
# 
# 
# - **For example**, we are saying that **odor = foul** is **indicative of** a **poisonous mushroom**. 
# <br><br> 
#   - **In** the **context** of a **categorical feature**, **odor** could **take** many **other values**. 
# <br><br>
#   - Since we **perturb** each **categorical feature drawing samples** according to the **original training distribution**, the way to interpret this is: **if odor was not foul**, on **average**, this **prediction** would be **0.27 less 'poisonous'**. 
#   
# <br>   
# - Let's **check** if **this** is the **case**:

# In[270]:


# Checking the different categories in the odor feature.

odor_idx = list(feature_names).index('odor')
explainer.categorical_names[odor_idx]


# In[271]:


# Checking the feature frequencies of different categories in the odor feature.

explainer.feature_frequencies[odor_idx]


# In[272]:


# Setting foul_idx equal to the index of 'foul' category in the odor feature.
# Then creating non_foul array with different categories in the odor feature except foul category.

foul_idx = 4
non_foul = np.delete(explainer.categorical_names[odor_idx], foul_idx)
non_foul


# In[273]:


# Creating non_foul_normalized_frequencies array with feature frequencies of different categories in the odor feature.
# Setting feature frequency of foul category to 0. Then normalizing the feature frequencies to have a total sum of 1.

non_foul_normalized_frequencies = explainer.feature_frequencies[odor_idx].copy()
non_foul_normalized_frequencies[foul_idx] = 0
non_foul_normalized_frequencies /= non_foul_normalized_frequencies.sum()
non_foul_normalized_frequencies


# In[286]:


# Calculating the probabilies of mushroom being poisonous for different values of odor except foul.
# Finally calculating the probability of mushroom being poisonous if odor not equal to foul.

print('Making odor not equal foul')

temp = X_test.values[i].copy()
print('P(poisonous) before:', predict_fn(temp.reshape(1,-1))[0,1], '\n')

average_poisonous = 0

for idx, (name, frequency) in enumerate(zip(explainer.categorical_names[odor_idx], non_foul_normalized_frequencies)):
    if name == 'foul':
        continue
    temp[odor_idx] = idx
    p_poisonous = predict_fn(temp.reshape(1,-1))[0,1]
    average_poisonous += p_poisonous * frequency
    print('P(poisonous | odor=%s): %.2f' % (name, p_poisonous))

print ('\nP(poisonous | odor != foul) = %.2f' % average_poisonous)


# - **Probability of poisonous when odor equals foul** = **1 - P(poisonous | odor != foul)** = **1 - 0.58** = **0.42**
# 
# 
# - We see that **in this** particular **case**, the **linear model** is **pretty close**: it **predicted** that **on average odor = foul increases** the **probability of poisonous by 0.27**, when **in fact it is by 0.42**. 
# 
# 
# - Notice though that **we only changed one feature (odor)**, when the **linear model takes into account perturbations of all** the **features at once**.