Using more sophisticated modeling packages (like XGBoost) to understand more complicated dynamics in the data
How to approach data exploration to understand more complicated relationships between the variables in your data
Why the "coherence" of a model is important - arguably, on the same level as its predictive performance
How to assess the "coherence" of a model using ICE plots
How to ascribe "reasons" to individual predictions
How to "consolidate" features to make the reasons more coherent and understandable
Using visualizations independently and from the SHAP package
In [1]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
plt.style.use('classic') # I prefer the appearance of the "old-school" color scheme
%matplotlib inline
from sklearn.metrics import roc_auc_score, confusion_matrix, log_loss, accuracy_score, r2_score, roc_curve, precision_recall_curve
import xgboost as xgb
pd.set_option("display.max_rows",999)
pd.set_option("display.max_columns",999)
np.set_printoptions(edgeitems=30, linewidth=100000)
In [2]:
# must have ml_insights installed ('pip install ml_insights')
import ml_insights as mli
In [3]:
mli.__version__
Out[3]:
You can get the data at:
https://drive.google.com/open?id=1wS-nOWQtZ3DA057or43E-zpm7duoqEMM
or
Download the file and place in the same directory as the notebook.
In [4]:
df_ins = pd.read_csv('ins_data.csv')
df_ins.info()
In [ ]:
This data set contains information about different insurance customers who have both home insurance and auto insurance with an agency. The information includes their premiums, coverages, vehicles, drivers, state of residence, and zip-code demographics.
In this exercise we will be building a model to predict who is likely to have an umbrella policy. Imagine, that our goal is to be able to recommend to new agencies (which currently don't sell umbrella policies) which customers are the best candidates. Our binary outcome variable is has_umbrella
We will focus on techniques for exploring, understanding, and interpreting data and models in a binary classification framework.
In [ ]:
In [5]:
df_ins.has_umbrella.value_counts(), np.mean(df_ins.has_umbrella)
Out[5]:
In [6]:
plt.hist(df_ins.homeown_premium)
Out[6]:
In [7]:
# Generally a good idea to put in custom bin widths
plt.hist(df_ins.homeown_premium, bins=np.linspace(0,10000,101));
In [8]:
bins_1 = np.linspace(0,2000,9)
bins_2 = np.array([2000,3000, 6000, 10000])
bins_final = np.unique(np.concatenate((bins_1, bins_2)))
mli.histogram_pair(df_ins.homeown_premium, df_ins.has_umbrella, bins=bins_final);
In [9]:
plt.hist(df_ins.auto_premium);
In [10]:
plt.hist(df_ins.auto_premium,np.linspace(0,15000,30+1));
In [ ]:
In [11]:
bins_1 = np.linspace(0,4000,8+1)
bins_2 = np.array([4000,6000, 8000, 15000])
bins_final = np.unique(np.concatenate((bins_1, bins_2)))
mli.histogram_pair(df_ins.auto_premium, df_ins.has_umbrella, bins=bins_final);
In [12]:
## Could do this on all the variables to understand...
In [ ]:
In [13]:
chosen_fold_variant = 'fold_num_1'
test_fold_num = 0
In [14]:
features_1 = ['auto_premium','homeown_premium']
In [ ]:
In [ ]:
In [15]:
# Define train and test sets
X_train_1 = df_ins.loc[df_ins[chosen_fold_variant]!=test_fold_num,features_1]
X_test_1 = df_ins.loc[df_ins[chosen_fold_variant]==test_fold_num,features_1]
y_train = df_ins.has_umbrella[df_ins[chosen_fold_variant]!=test_fold_num]
y_test = df_ins.has_umbrella[df_ins[chosen_fold_variant]==test_fold_num]
In [ ]:
In [16]:
xgb_def1 = xgb.XGBClassifier()
In [17]:
xgb_def1.fit(X_train_1,y_train)
pred_probs_def1 = xgb_def1.predict_proba(X_test_1)[:,1]
In [18]:
roc_auc_score(y_test, pred_probs_def1), log_loss(y_test, pred_probs_def1)
Out[18]:
In [19]:
mli.plot_pr_curves([y_test],[pred_probs_def1])
In [ ]:
In [ ]:
In [ ]:
In [20]:
plt.figure(figsize=(12,4))
bins1 = np.linspace(0,.2,8+1)
bins2 = np.linspace(.2,1,8+1)
bins_final = np.unique(np.concatenate((bins1, bins2)))
mli.plot_reliability_diagram(y_test,pred_probs_def1,size_points=False, bins=bins_final);
In [ ]:
In [21]:
### Define an example set
X_explore_1 = X_test_1.sample(100)
### "Mutate" the data points and make predictions on them
mxr_def_1 = mli.ModelXRay(xgb_def1,X_explore_1)
In [22]:
mxr_def_1.feature_dependence_plots(num_pts=1, y_scaling='logit')
Out[22]:
In [23]:
### Visualize the responses
indices = mxr_def_1.feature_dependence_plots(y_scaling='logit')
In [ ]:
In [ ]:
In [24]:
xgb_complex1 = xgb.XGBClassifier(n_estimators=1000, learning_rate=.01, gamma=5,
max_depth=2, subsample=.9,
reg_lambda=3)
In [25]:
xgb_complex1.fit(X_train_1,y_train, eval_set=[(X_test_1, y_test)],
eval_metric='logloss', early_stopping_rounds = 50)
pred_probs_complex_1 = xgb_complex1.predict_proba(X_test_1)[:,1]
In [26]:
## Performance of this model
roc_auc_score(y_test, pred_probs_complex_1), log_loss(y_test, pred_probs_complex_1)
Out[26]:
In [27]:
## Performance of previous model
roc_auc_score(y_test, pred_probs_def1), log_loss(y_test, pred_probs_def1)
Out[27]:
In [28]:
mli.plot_pr_curves([y_test,y_test],[pred_probs_def1,pred_probs_complex_1])
In [ ]:
In [29]:
bins1 = np.linspace(0,.2,8+1)
bins2 = np.linspace(.2,1,8+1)
bins_final = np.unique(np.concatenate((bins1, bins2)))
plt.figure(figsize=(10,4))
mli.plot_reliability_diagram(y_test,pred_probs_complex_1,size_points=False,
bins=bins_final);
In [30]:
### Examine Model More Deeply
mxr_complex_1 = mli.ModelXRay(xgb_complex1,X_explore_1)
In [31]:
## Result for more regularized model
mxr_complex_1.feature_dependence_plots(pts_selected=indices, y_scaling='logit')
Out[31]:
In [32]:
## Results for less regularized model
mxr_def_1.feature_dependence_plots(pts_selected=indices, y_scaling='logit')
Out[32]:
In [ ]:
In [ ]:
In [33]:
features_1a = ['auto_premium','homeown_premium'
,'avg_homeown_rate_in_state']
In [34]:
X_train_1a = df_ins.loc[df_ins[chosen_fold_variant]!=test_fold_num,features_1a]
X_test_1a = df_ins.loc[df_ins[chosen_fold_variant]==test_fold_num,features_1a]
In [35]:
xgb_complex1a= xgb.XGBClassifier(n_estimators=3000, learning_rate=.01, gamma=5, max_depth=3,
subsample=.8, reg_lambda=3)
In [36]:
xgb_complex1a.fit(X_train_1a,y_train, eval_metric='logloss',
eval_set=[(X_test_1a, y_test)], early_stopping_rounds=50)
Out[36]:
In [37]:
pred_probs_complex_1a = xgb_complex1a.predict_proba(X_test_1a)[:,1]
In [38]:
roc_auc_score(y_test, pred_probs_complex_1a), log_loss(y_test, pred_probs_complex_1a)
Out[38]:
In [39]:
mli.plot_pr_curves([y_test,y_test, y_test],[pred_probs_def1,pred_probs_complex_1, pred_probs_complex_1a])
In [ ]:
In [40]:
bins1 = np.linspace(0,.1,5)
bins2 = np.linspace(.2,1,9)
bins_final = np.unique(np.concatenate((bins1, bins2)))
plt.figure(figsize = (10,4))
mli.plot_reliability_diagram(y_test,pred_probs_complex_1a,size_points=False,
bins=bins_final);
In [41]:
X_explore_1a = X_test_1a.sample(100)
mxr1a = mli.ModelXRay(xgb_complex1a,X_explore_1a)
In [42]:
mxr1a.feature_dependence_plots(y_scaling='logit')
Out[42]:
In [ ]:
In [ ]:
In [ ]:
In [43]:
features_2 = ['auto_premium','homeown_premium', 'home_dwell_cov', 'home_pers_prop_cov',
'num_home_pol',
'yob_policyholder','min_vehicle_year', 'max_vehicle_year', 'num_vehicles',
'max_driver_yob', 'min_driver_yob',
'median_household_income','median_house_value'
,'avg_homeown_rate_in_state'
]
In [44]:
X_train_2 = df_ins.loc[df_ins[chosen_fold_variant]!=test_fold_num,features_2]
X_test_2 = df_ins.loc[df_ins[chosen_fold_variant]==test_fold_num,features_2]
In [45]:
xgb_complex2 = xgb.XGBClassifier(n_estimators = 3000, learning_rate=.01,
max_depth=3, subsample=.8, gamma=5,
reg_alpha=2, reg_lambda=5)
xgb_complex2.fit(X_train_2,y_train, eval_metric='logloss',
eval_set=[(X_test_2, y_test)], early_stopping_rounds=50)
pred_probs_complex2 = xgb_complex2.predict_proba(X_test_2)[:,1]
In [46]:
roc_auc_score(y_test, pred_probs_complex2), log_loss(y_test, pred_probs_complex2)
Out[46]:
In [47]:
mli.plot_pr_curves([y_test, y_test],[pred_probs_complex_1a,pred_probs_complex2])
In [48]:
bins1 = np.linspace(0,.1,5)
bins2 = np.linspace(.2,1,9)
bins_final = np.unique(np.concatenate((bins1, bins2)))
plt.figure(figsize = (10,4))
mli.plot_reliability_diagram(y_test,pred_probs_complex2,size_points=False,
bins=bins_final);
In [49]:
X_explore_2 = X_test_2[X_test_2.median_house_value<1000000].sample(100)
mxr_complex2 = mli.ModelXRay(xgb_complex2,X_explore_2)
In [50]:
mxr_complex2.feature_dependence_plots(pts_selected=indices, y_scaling='logit')
Out[50]:
In [ ]:
In [51]:
xgb_def2 = xgb.XGBClassifier()
xgb_def2.fit(X_train_2,y_train)
pred_probs_def2 = xgb_def2.predict_proba(X_test_2)[:,1]
In [52]:
roc_auc_score(y_test, pred_probs_def2), log_loss(y_test, pred_probs_def2)
Out[52]:
In [53]:
mli.plot_pr_curves([y_test, y_test, y_test],[pred_probs_complex_1a,pred_probs_complex2, pred_probs_def2])
In [54]:
bins1 = np.linspace(0,.1,5)
bins2 = np.linspace(.2,1,9)
bins_final = np.unique(np.concatenate((bins1, bins2)))
plt.figure(figsize = (10,4))
mli.plot_reliability_diagram(y_test,pred_probs_def2,size_points=False,
bins=bins_final);
In [55]:
mxr_def2 = mli.ModelXRay(xgb_def2,X_explore_2)
In [56]:
indices = mxr_def2.feature_dependence_plots(y_scaling='logit')
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
### In the previous, we focused our attention to understanding the overall dynamics of the model.
### However, in many cases, we want to understand the reasons behind a specific prediction.
#### arbitrarily order the features $X_1, X_2, \ldots, X_k$.
#### Compare $P(Y=1|X_1 = x_1, X_2 = x_2, \ldots, X_{j-1}=x_{j-1})$ to $P(Y=1|X_1 = x_1, X_2 = x_2, \ldots, X_j=x_j)$ (actually, compare the log odds) and attribute the difference to feature $X_j$ being equal to $x_j$.
#### Average the "attributions" to a feature over all possible orderings of the features
The "magic" of SHAP is that step (3) (averaging over all possible orderings) is typically computationally prohibitive, but for tree-based models they found a clever way to compute it by exploiting the structure.
Given a particular data point and its prediction, we can assign numerical values to each feature such that the sum of those values represents the distance between the "baseline" average and the predicted value.
A Unified Approach to Interpreting Model Predictions Scott Lundberg, Su-In Lee https://arxiv.org/abs/1705.07874
Consistent Individualized Feature Attribution for Tree Ensembles Scott M. Lundberg, Gabriel G. Erion, Su-In Lee https://arxiv.org/abs/1802.03888
In [ ]:
In [57]:
reas_df = mli.get_reason_score_matrix(xgb_complex2, X_test_2)
reas_df.round(decimals=2).head()
Out[57]:
In [58]:
## get_reason_matrix is essentially the following 3 lines
# X_test_dmat = xgb.DMatrix(X_test_2)
# reas_mat = xgb_complex2.get_booster().predict(X_test_dmat, pred_contribs=True)
# reas_df_2 = pd.DataFrame(reas_mat, columns = list(X_test_2.columns)+['Intercept'])
# reas_df_2.head()
In [59]:
## Demonstrate that these numbers sum up to the overall log_odds
log_odds_vec = np.sum(reas_df, axis=1)
pv = 1/(1+np.exp(-log_odds_vec))
np.sum(np.abs(pv- pred_probs_complex2)>.0001)
Out[59]:
In [ ]:
In [60]:
## Let's examine a specific case
row_num = 72
case_df = pd.DataFrame(X_test_2.iloc[row_num])
case_df['feat_importance'] = pd.DataFrame(reas_df.iloc[row_num])
case_df.columns = ['feat_value', 'feat_importance']
In [61]:
case_df.sort_values('feat_importance', ascending=False)
Out[61]:
In [62]:
# Get the average "impact" of each variable
reas_df.abs().mean().sort_values(ascending=False)
Out[62]:
In [63]:
## Look at correlations between the reasons
pd.DataFrame(np.round(np.corrcoef(reas_df.values[:,:-1].T), decimals=2), columns = reas_df.columns[:-1], index = reas_df.columns[:-1])
Out[63]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [64]:
def analyze_effect(feature_name, reas_df, test_df, fillna_val=0):
plt.subplot(2,2,1)
plt.scatter(test_df[feature_name].fillna(fillna_val), reas_df[feature_name], alpha=.1)
plt.subplot(2,2,2)
plt.scatter(reas_df[feature_name], test_df[feature_name].fillna(fillna_val), alpha=.1)
plt.subplot(2,2,3)
plt.hist(test_df[feature_name].fillna(fillna_val))
plt.subplot(2,2,4)
plt.hist(reas_df[feature_name])
In [65]:
plt.figure(figsize=(10,8))
analyze_effect('yob_policyholder',reas_df, X_test_2, fillna_val=1900)
In [66]:
plt.figure(figsize=(10,8))
analyze_effect('num_home_pol',reas_df, X_test_2, fillna_val=0)
In [67]:
# plt.figure(figsize=(10,8))
# analyze_effect('auto_premium',reas_df, X_test_2, fillna_val=-1000)
In [68]:
# plt.figure(figsize=(10,8))
# analyze_effect('median_house_value',reas_df, X_test_2, fillna_val=-1000)
In [69]:
# plt.figure(figsize=(10,8))
# analyze_effect('median_household_income',reas_df, X_test_2, fillna_val=-10000)
In [ ]:
Often, ascribing value to individual features may yield confusing results. For example, many individual features may represent the same concept or otherwise be highly correlated. In these cases, it may be somewhat arbitary which of the correlated features gets "credit" for the impact. In other cases, there may be many similar variables, each of which has a tiny impact, but collectively have a much greater impact. For this reason, we may want to "consolidate" the impact of features.
In [70]:
# Create a dictionary mapping the "group name" to the list of features included in that group
reason_mapping_umb = {
'Value_of_Real_Estate':['homeown_premium', 'home_dwell_cov','home_pers_prop_cov','num_home_pol'],
'State_Specific_Factors':['avg_homeown_rate_in_state'],
'Value_of_Automobiles':['min_vehicle_year', 'max_vehicle_year', 'num_vehicles','auto_premium'],
'Age_of_Policyholder_and_Family':['yob_policyholder', 'min_driver_yob', 'max_driver_yob'],
'Zipcode_Wealth':['median_household_income','median_house_value']}
In [71]:
cons_df = mli.consolidate_reason_scores(reas_df, reason_mapping_umb)
In [72]:
cons_df.head()
Out[72]:
In [73]:
cons_df.abs().mean().sort_values(ascending=False)
Out[73]:
In [ ]:
In [74]:
## Look at correlations between the reasons
pd.DataFrame(np.round(np.corrcoef(cons_df.values.T), decimals=2), columns = cons_df.columns, index = cons_df.columns)
Out[74]:
In [75]:
reason_string_vector = mli.predict_reason_strings(xgb_complex2, X_test_2, reason_mapping_umb, .3)
In [ ]:
In [76]:
pd.Series(reason_string_vector[pred_probs_complex2>.5]).value_counts()
Out[76]:
In [ ]:
In [ ]:
In [77]:
# if needed
# !pip install shap
In [78]:
import shap
shap.initjs()
#X_test_2_sm = X_test_2
In [79]:
explainer = shap.TreeExplainer(xgb_complex2)
In [80]:
shap_values = explainer.shap_values(X_test_2)
In [81]:
row_num=72
shap.force_plot(explainer.expected_value, shap_values[row_num,:], X_test_2.iloc[row_num,:], matplotlib=True)
In [82]:
cons_shap_df = mli.consolidate_reason_scores(pd.DataFrame(shap_values, columns=X_test_2.columns), reason_mapping_umb)
In [83]:
shap.force_plot(explainer.expected_value, cons_shap_df.iloc[row_num,:].values, feature_names=cons_shap_df.columns,matplotlib=True)
In [ ]:
In [ ]:
In [84]:
shap_interaction_values = shap.TreeExplainer(xgb_complex2).shap_interaction_values(X_test_2)
In [85]:
shap.summary_plot(shap_interaction_values, X_test_2)
In [86]:
shap.dependence_plot(
("median_household_income", "median_household_income"),
shap_interaction_values, X_test_2,
display_features=X_test_2
)
In [87]:
shap.dependence_plot(
("median_household_income", "yob_policyholder"),
shap_interaction_values, X_test_2,
display_features=X_test_2
)
In [ ]:
In [ ]:
In [ ]:
pip install structureboost
In [88]:
#!pip install structureboost
In [89]:
import structureboost as stb
In [90]:
features_3 = ['auto_premium','homeown_premium', 'home_dwell_cov', 'home_pers_prop_cov',
'num_home_pol',
'yob_policyholder','min_vehicle_year', 'max_vehicle_year', 'num_vehicles',
'max_driver_yob', 'min_driver_yob',
'median_household_income','median_house_value'
,'state_abbrev'
]
In [91]:
X_train_3 = df_ins.loc[df_ins[chosen_fold_variant]!=test_fold_num,features_3]
X_test_3 = df_ins.loc[df_ins[chosen_fold_variant]==test_fold_num,features_3]
In [92]:
US_50_graph = stb.graphs.US_50_and_DC_graph()
In [93]:
US_50_graph.vertices
Out[93]:
In [94]:
US_50_graph.edges
Out[94]:
In [95]:
def_config = stb.default_config_dict()
my_feat_config = stb.get_basic_config(X_train_3, def_config)
my_feat_config
Out[95]:
In [96]:
my_feat_config['state_abbrev']['graph'] = US_50_graph
In [97]:
my_stboost = stb.StructureBoost(num_trees=2000, learning_rate=.03,
feature_configs=my_feat_config,
max_depth=2, mode='classification')
In [98]:
my_stboost.fit(X_train_3, y_train, eval_set=(X_test_3, y_test), early_stop_past_steps = 5, eval_freq=20)
In [99]:
pred_probs_stb = my_stboost.predict(X_test_3)
In [100]:
# StructureBoost performance
roc_auc_score(y_test, pred_probs_stb), log_loss(y_test, pred_probs_stb)
Out[100]:
In [101]:
# Previous XGBoost performance
roc_auc_score(y_test, pred_probs_complex2), log_loss(y_test, pred_probs_complex2)
Out[101]:
In [ ]:
In [ ]:
In [ ]: