はじめに

SHAPを実務に使うためにやってみたメモ

環境

本記事では以下を事前にインストールしておく必要がある。

Python
numpy
pandas
matplotlib
SHAP
scikit-learn
LightGBM
RDKit

モデル構築

まずはRDKitを使ってモデルを作成してみよう。

まずはモジュールのインポート

import numpy as np
import pandas as pd
from rdkit import rdBase, Chem
from rdkit.Chem import AllChem, Descriptors
from rdkit.ML.Descriptors import MoleculeDescriptors
import lightgbm as lgb
from lightgbm import LGBMRegressor
from sklearn.feature_selection import VarianceThreshold
from sklearn.model_selection import train_test_split
import shap

続いて、smiles列のあるcsvファイルからデータを読み込み記述子計算および、目的変数の取得を行う。
今回はMoleculeNetのデータセットに含まれている溶解性の回帰モデルを試すこととする。

# ファイルの読み込み
import codecs
with codecs.open('../datas/raw/delaney-processed.csv', 'r', 'utf-8', 'ignore') as f:
    df = pd.read_csv(f)
# SMILESのデータを取り出す
smiles_list = df["smiles"].values

# SMILESのリストをMOLに変換
mols = []
for i, smiles in enumerate(smiles_list):
    mol = Chem.MolFromSmiles(smiles)
    if mol:
        mols.append(mol)
    else:
        print("{0} th mol is not valid".format(i))

# RDKitの記述子計算
descriptor_names = [descriptor_name[0] for descriptor_name in Descriptors._descList]
descriptor_calculator = MoleculeDescriptors.MolecularDescriptorCalculator(descriptor_names)
rdkit_descriptors_results = [descriptor_calculator.CalcDescriptors(mol) for mol in mols]
df_rdkit = pd.DataFrame(rdkit_descriptors_results, columns = descriptor_names, index=names_list)

# 目的変数の取得
df_t = df["measured log solubility in mols per litre"]

学習データとテストデータに分割しよう。

# 学習データとテストデータに分割
df_train_X, df_test_X, df_train_t, df_test_t = train_test_split(df_rdkit, df_t, test_size=0.3, random_state=0)

ここでのポイントはpandasのデータフレームの状態で分割することである。
pandasデータフレームのままSHAPのメソッドを呼び出すと、可視化の際にpandasの列名が特徴名として表示されるので大変便利である。

最後に学習してモデルを構築しよう。

model = LGBMRegressor()
model.fit(df_train_X, df_train_t)

これで準備は整った。

SHAPによるモデルの解釈

ほぼ、SHAPを使って機械学習モデルと対話するを参考にした。

下準備

下準備はこれだけだ。

shap.initjs()
explainer = shap.TreeExplainer(model=model)
shap_values = explainer.shap_values(X=df_test_X)

force_plot

まずはforce_plotによりテストデータの最初の1件だけのSHAPの結果を見てみよう。

shap.force_plot(explainer.expected_value, shap_values[0], df_test_X.iloc[0, :])

予測値は-1.43で、MolLogPという説明変数の影響により、予測値が平均(-3.046)よりも高くなっているということが見てとれる。
MolLogPが大きいから予測値が高いといっているわけではないことに注意しよう。

続いてテストデータの最初の100件のSHAPの結果をまとめて見てみよう。

shap.force_plot(explainer.expected_value, shap_values[0:100], df_test_X.iloc[0:100,:])

予測値が高い順から低い順に100件のSHAPの結果が描画される。
一番太い帯で表示されているのがMolLogPで、予測値が-5よりも小さい場合、特にMolLogPが大きな影響を与えているということが見て取れる。

次にY軸を「MolLogP effects」としてみよう。

こうすると、MolLogPの効果だけを取り出してみることができる。特定の説明変数の影響をみるにはこちらが便利だ。

sammary_plot

続いて sammary_plotを使ってみよう。

shap.summary_plot(shap_values, df_test_X)

これの何が便利かというと、説明変数の大きさが予測値にどう影響を与えているのかが分かる点だ。
例えばMolLogPについては、大きい(=ピンクが濃い)程、予測値にマイナスの影響を与えていることが分かる。

dependence_plot

続いてdependancy_plotを見てみよう。

shap.dependence_plot("MolLogP", shap_values, df_test_X, interaction_index=10)

これはある説明変数の予測値への影響が、他の説明変数に依存しているかを確認することができる。
予測のメカニズムを解析する際に、非常に役立ちそうである。

decision_plot

続いてdecision_plotを見てみよう。これはデータ毎の説明変数の影響を重ねてみることができる。
テストデータの最初の100件を可視化してみよう。

shap.decision_plot(explainer.expected_value, shap_values[0:100], df_test_X.iloc[0:100,:])

feature_order='hclust'を指定すると、階層的クラスタリングに基づき説明変数が表示され、似たような傾向の予測パスをグループ化して表示することができる。

shap.decision_plot(explainer.expected_value, shap_values[0:100], df_test_X.iloc[0:100,:], feature_order='hclust')

おわりに

ここで紹介しきれなかったものもあるので、リファレンスマニュアルを見て色々試してみたい。
予測結果の妥当性の確認にも使えそうな気がする。

物性予測モデルをSHAPで解釈してみる