<jupyter_start><jupyter_text>Lambda School Data Science*Unit 2, Sprint 3, Module 4*---  Model Interpretation- Visualize and interpret **partial dependence plots**- Explain individual predictions with **shapley value plots**  SetupLibraries:- category_encoders- matplotlib- numpy- pandas- xgboost- scikit-learn- [**shap**](https://github.com/slundberg/shap)- [**pdpbox**](https://github.com/SauceCat/PDPbox)<jupyter_code># ! pipenv install shap pdpbox
import IPython
DATA_PATH = '../data/'

# Ignore this warning: https://github.com/dmlc/xgboost/issues/4300
# xgboost/core.py:587: FutureWarning: Series.base is deprecated and will be removed in a future version
import warnings
warnings.filterwarnings(action='ignore', category=FutureWarning, module='xgboost')<jupyter_output><empty_output><jupyter_text>I. Wrangle the data<jupyter_code># data analysis and wrangling
import pandas as pd
import numpy as np

# visualization
import matplotlib.pyplot as plt
%matplotlib inline

# encoders
from category_encoders import OneHotEncoder, OrdinalEncoder
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler

#metrics
from sklearn.metrics import r2_score

#pipeline
from sklearn.pipeline import make_pipeline

# machine learning
from sklearn.linear_model import LinearRegression
def wrangle(DATA_PATH):
  df = pd.read_csv(DATA_PATH+'lending-club/lending-club-subset.csv', parse_dates=['issue_d']).set_index('issue_d')

  # Just use 36 month loans
  df = df[df.term==' 36 months']

  # Index & sort by issue date
  df = df.sort_index()

  # Clean data, engineer feature, & select subset of features
  df = df.rename(columns=                     
      {'annual_inc': 'annual_income', 
      'fico_range_high': 'credit_score', 
      'funded_amnt': 'loan_amount', 
      'title': 'loan_purpose'})

  df['interest_rate'] = df['int_rate'].str.strip('%').astype(float)
  df['monthly_debts'] = df['annual_income'] / 12 * df['dti'] / 100

  columns = ['annual_income', 
            'credit_score', 
            'loan_amount', 
            'loan_purpose', 
            'monthly_debts', 
            'interest_rate']

  df = df[columns]
  df = df.dropna()


  return df
df = wrangle(DATA_PATH)
df<jupyter_output><empty_output><jupyter_text>- Turn the high cardinality column to lower cardinality column<jupyter_code>df['loan_purpose'].nunique()
# for loop
# apply a function to this column
top_5_purposes = df['loan_purpose'].value_counts()[:5].index
def top_5(purpose):
    if purpose in top_5_purposes:
        return purpose
    else:
        return 'Other_purpose'
df['loan_purpose'] = df['loan_purpose'].apply(top_5)
df['loan_purpose'].value_counts()<jupyter_output><empty_output><jupyter_text>II. Split the data<jupyter_code>X = df.drop(columns='interest_rate')
y = df['interest_rate']
cutoff1 = int(len(df)* 0.6)
cutoff2 = int(len(df)* 0.8)

X_train, y_train = X.iloc[:cutoff1], y.iloc[:cutoff1]
X_val, y_val = X.iloc[cutoff1:cutoff2], y.iloc[cutoff1:cutoff2]
X_test, y_test = X.iloc[cutoff2:], y.iloc[cutoff2:]<jupyter_output><empty_output><jupyter_text>III. Set the baseline<jupyter_code>from sklearn.metrics import r2_score
r2_score(y_train, [y_train.mean()]*len(y_train))<jupyter_output><empty_output><jupyter_text>IV.  Build the models  Fit Linear Regression model<jupyter_code>model_lr = make_pipeline(OneHotEncoder(use_cat_names=True),
                         SimpleImputer(),
                         LinearRegression())
model_lr.fit(X_train,y_train)
r2_score(y_val, model_lr.predict(X_val))
column = model_lr.named_steps['onehotencoder'].feature_names
coef = model_lr.named_steps['linearregression'].coef_
pd.Series(coef, index=column).sort_values()
# if the income increase 10000, the interest rate will decrease 0.02%
delta_income = 10000
delta_interest_rate = delta_income*(-0.000002)
delta_interest_rate<jupyter_output><empty_output><jupyter_text>Fit Gradient Boosting model<jupyter_code>from xgboost import XGBRegressor
model_gb = make_pipeline(OrdinalEncoder(),
                         XGBRegressor(n_estimators=200, 
                                      objective='reg:squarederror', 
                                      n_jobs=-1))
model_gb.fit(X_train, y_train)
y_pred = model_gb.predict(X_val)
print('Gradient Boosting R^2', r2_score(y_val, y_pred))<jupyter_output>Gradient Boosting R^2 0.2900706080657681<jupyter_text>Explaining Gradient Boosting???Linear models have coefficients, but trees do not.Instead, to see the relationship between individual feature(s) and the target, we can use partial dependence plots.  Visualize and interpret partial dependence plots OverviewPartial dependence plots show the relationship between 1-2 individual features and the target — how predictions partially depend on the isolated features. It's explained well by [PDPbox library documentation](https://pdpbox.readthedocs.io/en/latest/):>**The common headache**: When using black box machine learning algorithms like random forest and boosting, it is hard to understand the relations between predictors and model outcome. For example, in terms of random forest, all we get is the feature importance. Although we can know which feature is significantly influencing the outcome based on the importance calculation, it really sucks that we don’t know in which direction it is influencing. And in most of the real cases, the effect is non-monotonic. We need some powerful tools to help understanding the complex relations between predictors and model prediction.Let's also look at an [animation by <NAME>](https://twitter.com/ChristophMolnar/status/1066398522608635904), author of [_Interpretable Machine Learning_](https://christophm.github.io/interpretable-ml-book/pdp.htmlexamples):> Partial dependence plots show how a feature affects predictions of a Machine Learning model on average.> 1. Define grid along feature> 2. Model predictions at grid points> 3. Line per data instance -> ICE (Individual Conditional Expectation) curve> 4. Average curves to get a PDP (Partial Dependence Plot)  Partial Dependence Plots with 1 featurePDPbox- [Gallery](https://github.com/SauceCat/PDPboxgallery)- [API Reference: pdp_isolate](https://pdpbox.readthedocs.io/en/latest/pdp_isolate.html)- [API Reference: pdp_plot](https://pdpbox.readthedocs.io/en/latest/pdp_plot.html)<jupyter_code># Later, when you save matplotlib images to include in blog posts or web apps,
# increase the dots per inch (double it), so the text isn't so fuzzy
import matplotlib.pyplot as plt
plt.rcParams['figure.dpi'] = 72
from pdpbox.pdp import pdp_isolate, pdp_plot
feature = 'annual_income'

isolated = pdp_isolate(model=model_gb,
           dataset=X_val, 
           model_features=X_val.columns,
           feature= feature)

pdp_plot(isolated,feature_name=feature,plot_lines=True, frac_to_plot=0.01)
pdp_plot(isolated,feature_name=feature)
# just check the range that we are intested in 
pdp_plot(isolated,feature_name=feature)
plt.xlim((0, 50000));<jupyter_output><empty_output><jupyter_text>Partial Dependence Plots with 2 featuresSee interactions!PDPbox- [Gallery](https://github.com/SauceCat/PDPboxgallery)- [API Reference: pdp_interact](https://pdpbox.readthedocs.io/en/latest/pdp_interact.html)- [API Reference: pdp_interact_plot](https://pdpbox.readthedocs.io/en/latest/pdp_interact_plot.html)Be aware of a bug in PDPBox version <= 0.20 with some versions of matplotlib:- With the `pdp_interact_plot` function, `plot_type='contour'` gets an error, but `plot_type='grid'` works- This will be fixed in the next release of PDPbox: https://github.com/SauceCat/PDPbox/issues/40<jupyter_code>from pdpbox.pdp import pdp_interact, pdp_interact_plot
features = ['annual_income', 'credit_score']

interact = pdp_interact(model=model_gb,
           dataset=X_val, 
           model_features=X_val.columns,
           features= features)

pdp_interact_plot(interact, feature_names=features, plot_type='grid')<jupyter_output><empty_output><jupyter_text>BONUS: PDPs with categorical featuresJust for your future reference, here's a bonus example to demonstrate partial dependence plots with categorical features.1. I recommend you use Ordinal Encoder or Target Encoder, outside of a pipeline, to encode your data first. (If there is a natural ordering, then take the time to encode it that way, instead of random integers.) Then use the encoded data with pdpbox.2. There's some extra work to get readable category names on your plot, instead of integer category codes.<jupyter_code># retrain the model step by step, do not use pipeline
oe = OrdinalEncoder()
X_train_T = oe.fit_transform(X_train)
X_val_T = oe.transform(X_val)
model_gb = XGBRegressor(n_estimators=200, objective='reg:squarederror', n_jobs=-1)
model_gb.fit(X_train_T, y_train)
# Use Pdpbox
%matplotlib inline
import matplotlib.pyplot as plt
feature = 'loan_purpose'
isolated = pdp_isolate(model=model_gb, dataset=X_val_T, model_features=X_val_T.columns, feature=feature)
pdp_plot(isolated, feature);
# use the for loop below to map the numbers to class names
feature = 'loan_purpose'
for item in oe.mapping:
    if item['col'] == feature:
        feature_mapping = item['mapping']
        
feature_mapping = feature_mapping[feature_mapping.index.dropna()]
category_names = feature_mapping.index.tolist()
category_codes = feature_mapping.values.tolist()
feature_mapping = feature_mapping[feature_mapping.index.dropna()]
category_names = feature_mapping.index.tolist()
category_codes = feature_mapping.values.tolist()
pdp_plot(isolated, feature)

# Automatically change the xticks labels
plt.xticks(category_codes, category_names);<jupyter_output><empty_output><jupyter_text>Explain individual predictions with shapley value plots  Overview We’ll use TreeExplainer from an awesome library called [SHAP](https://github.com/slundberg/shap), for “additive explanations” — we can explain individual predictions by seeing how the features add up!<jupyter_code>oe = OrdinalEncoder()
X_train_T = oe.fit_transform(X_train)
model_gb = XGBRegressor(n_estimators=200, objective='reg:squarederror', n_jobs=-1)
model_gb.fit(X_train_T, y_train)
# get one observation from X_test
row = X_test.iloc[[1]]
row = oe.transform(row)
row
model_gb.predict(row)
import shap

explainer = shap.TreeExplainer(model_gb) # do not like pipeline
shap_values = explainer.shap_values(row)
shap_values 
# The expected value for target variable
y_train.mean(), model_gb.predict(X_train_T).mean()
shap.initjs()
shap.force_plot(
                base_value=explainer.expected_value,
                shap_values=shap_values,
                features=row)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09917978693315735, 'fraction_numerical': 0.01319883096068634, 'mean_word_length': 5.221700879765396, 'pattern_counts': {'":': 0, '<': 46, '<?xml version=': 0, '>': 52, 'https://': 14, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '147768', 'license': 'MIT', 'n_tokens_mistral': 3560, 'n_tokens_neox': 3327, 'n_words': 893, 'path': 'module4-model-interpretation/Module234_03_25.ipynb', 'repo_name': 'csslambda/Unit2-sprint3-ds25', 'size': '10617'}