How to calculate joint feature contribution for XGBoost Classifier in python? I referred to this beautiful document to research about joint feature contibutions. But this works only for RandomForest algorithms because of treeinterpreter (does not work with xgboost). Is there a similar way out for XGBoost as well?
Basically what I want to achieve is to find out the joint contributions of all the combination of features towards the prediction. For example if I have a, b and c as my features, I want to know what is the effect of ab, bc and ca towards the prediction result. It is very similar to shap and lime, but for combinations of features.
I did some research and learnt about xgbfir package. It gives the joint contributions into an excel file. You can set the level of interaction with this. I wrote some code around it to generate a plot that solves the purpose. 

If the package is not installed

pip install xgbfir


After the installation:

import xgbfir
from matplotlib import pyplot as plt

xgbfir.saveXgbFI(model, feature_names=X.columns, OutputXlsxFile='FI.xlsx')

joint_contrib = pd.read_excel('FI.xlsx')

xls = pd.ExcelFile('FI.xlsx')
df1 = pd.read_excel(xls, 'Interaction Depth 0')
df2 = pd.read_excel(xls, 'Interaction Depth 1')
df3 = pd.read_excel(xls, 'Interaction Depth 2')

frames = [df1, df2, df3]
joint_contrib = pd.concat(frames)

joint_contrib=joint_contrib.sort_values(by='Gain', ascending=True)
joint_contrib=joint_contrib.head(20)

height = joint_contrib['Gain']
bars = joint_contrib['Interaction']
y_pos = np.arange(len(bars))

plt.barh(y_pos, height)
plt.yticks(y_pos, bars)
plt.show()


This will give the top 20 feature interactions in terms of gain. 

Thanks to Philip Cho who introduced me to xgbfir. 

Follow the link for more information regarding xgbfir