.ipynb .pdf

repository open issue

Colab

Contents

12. 予測を説明する

ニューラルネットワークの予測は一般に解釈可能ではありません。この章ではその予測を説明する方法を探ります。これは説明可能なAI（XAI）という広範なトピックの一部であり、なぜ特定の予測がなされるのかを理解するにの役立つはずです。モデルの予測を理解できるようになれば、実用的、理論的、かつ規制の観点で正当化することにつながります。そのため重要なトピックとなっています。またその根拠を理解できれば、モデルの予測を利用する可能性が高くなることが示されており[LS04]、確かに実用的です。もう一つの実用的な関心事は、モデルがどのようにその予測に至ったのかを理解できれば、方法を正確に実装することがはるかに容易になるということです。透明性の理論的な正当性は、モデルドメインの不完全性（すなわち、共変量シフト）[DVK17]を識別するのに役立ちます。最近、欧州連合[GF17]とG20[Dev19]は、機械予測に関する説明を推奨または要求するガイドラインを採択しました。そのため、これは現在、コンプライアンスの問題になってきています。欧州連合はさらに踏み込んで、より厳しい法律案を検討しています。

読者層と目的

この章はStandard LayersとDeep Learning on Sequencesに基づいています。また、条件付き確率を含む確率論に関する十分な知識があることを前提としています。そうでない場合は、私のノートや入門的な確率のテキストを読んで概要を理解することをお勧めします。この章を読むと、以下のことができるようになると想定されます。

• なぜ説明が重要かを正当化できる

• 正当化、解釈、説明の区別ができる

• 特徴量重要度とシャープレイ値を計算できる

• 反実仮想を定義し、それを計算できる

• どのモデルが解釈可能で、どのように解釈可能なサロゲートモデルを適合させるかを理解する

説明可能なAIの必要性に関する有名な例として、肺炎でERに運ばれた患者の死亡リスクを評価するML予測器を構築したCaruanaらの[CLG+15]が挙げられます。この論文のアイデアは、肺炎の患者をこのツールでスクリーニングすれば、医師がどの患者がより死亡リスクが高いかを知るのに役立つというものでした。そのモデルはかなり正確でしたが、その予測の解釈を調べたところ、医学的におかしな推論がなされていました。驚くべきことに、そのモデルは喘息持ちの患者（喘息患者という）が肺炎でERに来院した場合、死亡リスクが低下することを示唆していました。喘息というのは呼吸が困難になる病気であるにも関わらず、肺炎患者の死亡率が低くなるということがわかったのです。

この結果は偶発的なものでした。喘息は実際には肺炎での死亡リスクが高いのですが、医師はそのことを実感しているため、彼らに対してより積極的で丁寧な対応をしていたのです。すなわち、喘息患者に対して医師がより丁寧なケアと配慮をしていたため、死亡者数が少なくなっていたのです。経験則から言えば、モデルの予測は正しいです。しかし、もしこのモデルが実用化されれば、喘息患者を死亡リスクが低いと誤って判断し、喘息患者は本来受けられたはずのケアを受けられずに命を落とす可能性がありました。幸いにも研究者はモデルの解釈可能性によってこの問題を特定し、喘息患者を命の危険にさらすリスクを回避することができました。このように、解釈は常に予測モデルを構築する際に考慮すべきステップであることがわかります。

12.1. 説明とは何か

ここでは、Miller[Mil19]の説明の定義を使います。Millerは、解釈可能性、正当化、説明を以下の定義で区別しています。

• 解釈可能性：観察者が判断の原因を理解できる度合いを指します。Millerはこれを説明可能性と同義とみなしました。これは一般的にモデルの特性です。

• 正当化：モデルのテストエラーや正確度のように、なぜその決定が良いのかについての証拠や説明です。これはモデルの特性です。

• 説明：結果の文脈と原因を与える、人間を対象とした情報の提示です。これらが本章の主要な論点です。これは一般的にモデルの特性ではなく、新たに生成する追加的な情報です。

説明の構成要素について詳しく説明しますが、説明は予測を正当化することとは異なることに注意してください。正当化とは既に見たように、モデルの予測が正確であると信じるべき理由についての経験的な証拠です。一方、説明とは予測の原因を明らかにすることであり、最終的に人に理解されることを目的としています。

ディープラーニングはそれだけではブラックボックス的なモデリング手法です。解釈可能性も説明可能性もありません。重みやモデル式を調べても、なぜその予測がなされるのかについての洞察はほとんど得られません。解釈可能性はディープラーニングに対する追加的なタスクであり、モデルの予測に説明を加えることを意味します。しかしこれは難しい問題です。ディープラーニングのブラックボックス的な性質に加えて、モデルの予測に対する説明とは何かについてのコンセンサスが取れていないためです[DVK17]。ある人は解釈可能性に各予測を正当化する自然言語の説明を期待しますが、ある人はどの特徴が予測に最も貢献したかを示すだけで十分と考えます。

MLモデルの解釈には2つのアプローチがあります。説明による事後解釈と自己説明型モデル[MSK+19]です。自己説明型モデルは，専門家がモデルの出力を見て，論理的に特徴量と結びつけることができるように構築されており、本質的に解釈可能です。ただしタスクモデルに強く依存します[MSMuller18]。身近な例では、分子動力学や一点量子エネルギー計算のような物理学に基づくシミュレーションがあります。分子動力学の軌跡を調べ、出力された数値を見て、例えば薬物分子がタンパク質に結合すると予測する理由を説明することができます。

自己説明型モデルはディープラーニングの解釈には役に立たない/関連付かないように思えるかもしれません。しかし、後の節で、自己説明型のサロゲートモデル (またはプロキシモデル）を作り、ディープラーニングモデルと一致するように訓練すればよいことがわかるでしょう。最初からサロゲートモデルを使わずにディープラーニングモデルを介するのは訓練コストを減らせるからですが、それはなぜでしょうか。学習済みニューラルネットワークは任意の点をラベル付けできる、すなわち学習データを無限に生成することができるからです。サロゲートモデルの他に、Attention機構[BCB14]のように自己説明的な特徴を内包したディープラーニングモデルを構築することもできます。Attention機構に基づいて入力特徴量と予測値を結びつけることができます。また機械学習にはシンボリック回帰というものがあり、直接解釈できる数式を扱うことで自己説明的なモデルを構築しようとします[AGFW21, BD00, UT20]。その特性から、シンボリック回帰はサロゲートモデルを生成するために用いられます[CSGB+20]。

説明による事後解釈には様々なアプローチがあります。代表的なものは学習データの重要度、特徴量重要度、反実仮想的な説明[WSW22, RSG16a, RSG16b, WMR17]です。データの重要度に基づく事後解釈の例は、予測を説明する最も影響力のある学習データを特定することです[KL17]。それによって説明がつくかどうかは議論の余地がありますが、どのデータが予測に関連しているかを理解するのに役立つことは確かです。特徴量重要度はおそらく最も一般的なXAIアプローチで、コンピュータビジョンの研究に頻繁に登場し、例えば画像の分類にとって最も重要なピクセルをハイライトします。

反実仮想的な説明は事後解釈の新しい方法です。反実仮想は説明として機能する新しいデータ点です。反実仮想は、その特徴量がどれほど重要で敏感であるかについての洞察を与えます。例として融資を勧めるモデルがあるとします。そのモデルは以下の反実仮想的な説明を生成することができます（[WMR17]より）。

あなたは、年収、郵便番号、資産に基づいてローンを拒否されました。もし、あなたの年収が45,000ドルであれば、あなたはローンを提供されたでしょう。

2番目の文が反実仮想であり、特徴量をどのように変えればモデルの結果に影響を与えるかを示しています。反実仮想は複雑さと説明力を良いバランスで提供します。

以上が広範なXAI分野に関する概観でした。解釈可能なディープラーニングについての最近のレビューはSamekらの[SML+21]を見てください。またディープラーニングを含む解釈可能な機械学習に関する網羅的な情報はChristopher Molnarがオンラインブックで公開しています。予測誤差や予測の信頼性は正当化の意味合いが強いのでここでは扱いませんが、Regression & Model Assessmentの手法が適用できるので参照してください。

12.2. 特徴量重要度

特徴量重要度は、機械学習モデルを解釈する上で最もわかりやすく、最も一般的な方法です。特徴量重要度の出力は各特徴量に対するランキングまたは数値であり、通常は単一の予測に対するものです。モデル全体の特徴量重要度は大域的特徴量重要度と呼ばれ、単一の予測に対しては局所的特徴量重要度と呼ばれます。大域的な特徴量重要度と解釈可能性を持つことは比較的まれです。正確なディープラーニングモデルは特徴空間の位置によって重要な特徴量が変わるためです。

まずは線形モデルで特徴量重要度を見てみましょう。

(12.1)\[\begin{equation} \hat{y} = \vec{w}\vec{x} + b \end{equation}\]

ここで \(\vec{x}\)は特徴量ベクトルです。特徴量重要度を評価する簡単な方法は、特定の特徴量\(x_i\)に関する重み\(w_i\)を単に見ることです。この重み\(w_i\)は、他のすべての特徴が一定で、\(x_i\)が1増加した場合にどの程度変化するかを示しています。もし、特徴量の大きさが同程度であれば、この方法は特徴量の順位付けとして機能するでしょう。しかし、特徴量が単位を持つ場合、単位の選択と特徴量の相対的な大きさに影響されます。例えば、気温が摂氏から華氏に変更された場合、1度上昇した時の影響は小さくなります。

特徴量の大きさや単位の影響を排除して特徴量重要度を評価する少し良い方法は、\(w_i\)を特徴量の標準偏差で割ることです。標準偏差とは、予測値の二乗誤差の総和を偏差平方和で割った値です。すなわち標準偏差は予測の正確度と特徴量の分散の比です。標準偏差で割った\(w_i\)はt-分布と比較できるため、\(t\)-統計量と呼ばれます。

(12.2)\[\begin{equation} t_i = \frac{w_i}{S_{w_i}},\; S^2_{w_i} = \frac{1}{N - D}\sum_j \frac{\left(\hat{y}_j - y_j\right)^2}{\left(x_{ij} - \bar{x}_i\right)^2} \end{equation}\]

ここで、Nは例数、 Dは特徴量数、 \(\bar{x}_i\)はi番目の特徴量の平均値です。\(t_i\)値は、特徴量の順位付けと仮説検定に利用できます。もし\(P(t > t_i) < 0.05\)であれば、その特徴量は統計的に有意で、\(P(t)\)はStudent’s \(t\)-分布です。特徴量の有意性はモデルに存在する他の特徴量に依存することに注意してください。つまり新しい特徴量を追加すると、一部が冗長になる可能性があります。

次に非線形の場合を見てみましょう。非線形学習関数\(\hat{f}(\vec{x})\)では、特徴量が1増加した場合に予測がどのように変化するかを微分近似で計算します。

\[ \frac{\Delta \hat{f}(\vec{x})}{\Delta x_i} \approx \frac{\partial \hat{f}(\vec{x})}{\partial x_i} \]

1だけ変化させると

(12.3)\[\begin{equation} \Delta \hat{f}(\vec{x}) \approx \frac{\partial \hat{f}(\vec{x})}{\partial x_i}. \end{equation}\]

実際にはこの式を少し変形します。0を中心としたテイラー級数を使う代わりに、他のルート（関数が0となる点）を中心にします。そうすることで、決定境界（ルート）で系列を「接地」し、予測されるクラスを決定境界から「遠ざけ」たり「近づけ」たりすることができます。もう1つの方法は、テイラー級数の1次の項を使用して線形モデルを構築するというものです。そして、その線形モデルに上記と同様のことを行い、その係数を特徴量の「重要度」として使用します。具体的には、\(\hat{f}(\vec{x})\)に対して以下のようなサロゲート関数を使用します。

(12.4)\[\begin{equation} \require{cancel} \hat{f}(\vec{x}) \approx \cancelto{0}{f(\vec{x}')} + \nabla\hat{f}(\vec{x}')\cdot\left(\vec{x} - \vec{x}'\right) \end{equation}\]

ここで、\(\vec{x}'\)は\(\hat{f}(\vec{x})\)のルートです。たいてい自明なルート\(\vec{x}' = \vec{0}\)を選択するかもしれませんが、近傍のルートが理想的です。このルートはしばしばベースライン入力と呼ばれます。上記の線形の例とは対照的に、部分的な\(\frac{\partial \hat{f}(\vec{x})}{\partial x_i}\)の積とベースラインからの増分\((x_i - x_i')\)を考えます。

12.2.1. ニューラルネットワークの特徴量重要度

ニューラルネットワークでは、偏導関数は出力に対する実際の変化を近似するには不十分です。入力に対する小さな変化が不連続な場合（ReLUのような非線形性のため）、ほとんど説明力を持たなくなることがあります。これはshattered gradients問題[BFL+17]と呼ばれています。また個々の特徴量に分けると、特徴量間の相関も欠落してしまいます。これは線形モデルにはない問題です。したがって、微分近似は局所的な線形モデルでは十分に機能しますが、ディープニューラルネットワークでは機能しません。

ニューラルネットワークにおけるshattered gradients問題を回避する方法はいろいろあります。よく使われるのはintegrated gradients [STY17] とSmoothGrad[STK+17]の2つの方法です。integrated gradientsは\(\vec{x}'\)から\(\vec{x}\)まで直線で結ぶ経路を考え、この経路上で対象となる変数の微分値を積分で統合します。

(12.5)\[\begin{equation} \textrm{IG}_i = \left(\vec{x} - \vec{x}'\right) \int_0^1\left[\nabla\hat{f}\left(\vec{x}' + t\left(\vec{x} - \vec{x}'\right)\right)\right]_i\,dt \end{equation}\]

\(t\)は経路に沿ったある増分で、\(t = 0\) のとき \(\vec{x}' + t\left(\vec{x} - \vec{x}'\right) = \vec{x}'\)、\(t = 1\) のとき \(\vec{x}' + t\left(\vec{x} - \vec{x}'\right) = \vec{x}\) です。この式により各特徴量\(i\)のintegrated gradientを得ます。integrated gradientは各特徴量の重要度ですが、shattered gradientsの持つ複雑さはありません。またモデル \(f(\vec{x})\) がほとんど至るところで微分可能であれば、 \(\sum_i \textrm{IG}_i = f(\vec{x}) - f(\vec{x}')\) という式が成立します。これはintegrated gradientsで計算された各特徴量重要度の合計値が、ベースラインと予測値の差に等しくなることを意味しています。すなわちベースラインから予測値の変化量を完全に分離してくれます[STY17]

integrated gradientsの実装は比較的簡単です。経路を入力特徴量 \(\vec{x}\)とベースライン \(\vec{x}'\)の間にある離散入力の集合に分割することにより、リーマン和で経路の積分を近似します。これらの入力の勾配をニューラルネットワークで計算します。そして、ベースラインからの特徴量の変化量\(\left(\vec{x} - \vec{x}'\right)\)を乗じます。

SmmothGradはintegrated gradientsと同様の考え方です。しかし経路にそった勾配を合計するのではなく、予測の近くにあるランダムな点から勾配を計算します。式は以下の通りです。

(12.6)\[\begin{equation} \textrm{SG}_i = \sum_j^M\left[\nabla\hat{f}\left(\vec{x}' + \vec{\epsilon}\right)\right]_i \end{equation}\]

\(M\)はサンプル数の選択であり、\(\vec{\epsilon}\)は\(D\)ゼロ平均ガウシアンからサンプリングされます[STK+17]。 ここでの実装上の唯一の変更点は、経路を一連のランダムな摂動に置き換えることです。

これらの勾配ベースの方法以外にも、Layer-wise Relevance Propagation (LRP)はニューラルネットワークにおける特徴量重要度の解析の一般的な方法です。LPRは、1つの層の出力値を入力特徴量に分割するニューラルネットワークを介した逆伝播を行うことで機能します。これは「関連性を分散させる」ということです。LPRの変わったところは、各層の種類毎に独自の実装が必要なことです。解析的な導関数に頼らず、層の方程式のテイラー級数展開で対応します。GNNやシーケンスモデル用のLRPもあり、LRPは材料や化学のほとんどの場面で使うことができます[MBL+19]。

12.2.2. シャープレイ値

モデル非依存的に特徴量重要度を扱う方法として、シャープレイ値があります。シャープレイ値はゲーム理論に由来するもので、協力的なプレーヤーに、その貢献度に応じて報酬を支払う方法についての解決策です。各特徴量がプレーヤーであり、予測値への貢献度に応じて「支払う」ことを想定しています。シャープレイ値 \(\phi_i(x)\)は、インスタンス\(x\)の特徴量\(i\)に対する支払いです。予測関数値 \(\hat{f}(x)\)をシャープレイ値に分割して、その和が関数値\(\sum_i \phi_i(x) = \hat{f}(x)\)となるようにします。つまり、ある特徴量のシャープレイ値は予測に対する数値的な貢献度と解釈できます。シャープレイ値の強力な利点は、モデルに依存せず、予測値を各特徴量に分割でき、予測の説明に必要な属性（対称性、線形性、順序不変性など）を持つことです。欠点は、厳密な計算には特徴量の組み合わせの数だけコストがかかること、スパース性を持たないことであり、結果的に特徴量数の増加に伴って有用性が低くなります。ここで紹介する手法もスパース性を持たないものがほとんどです。L1正則化(Standard Layers参照)のように、常にモデルをスパースにすることでスパースな説明を実現することができます。

シャープレイ値は次のように計算されます。

(12.7)\[\begin{equation} \phi_i(x) = \frac{1}{Z}\sum_{S \in N \backslash x_i}v(S\cup x_i) - v(S) \end{equation}\]

\[ Z = \frac{|S|!\left(N - |S| - 1\right)!}{N!} \]

\(S \in N \backslash x_i\)は特徴量\(x_i\)を除いた全ての特徴量の集合を意味し、\(S\cup x_i\)は特徴量 \(x_i\)を集合に戻すことを意味します。また\(v(S)\)は\(S\)に含まれる特徴量のみを使用した場合の\(\hat{f}(x)\)の値であり、\(Z\)は正規化用の値です。この式は、特徴量\(i\)を追加/削除することによって形成される\(\hat{f}\)の取りうるすべての差の平均と解釈することができます。

しかし、特徴量\(i\)をモデル式からどのように「取り除く」ことができるでしょうか。特徴量\(i\)を無用のものとして扱う（周辺化する）ことでできます。周辺化とは確率変数\(P(x) = \int\, P(x,y)\,dy\)を積分する方法であることを思い出してください。これは取りうるすべての値\(x\)を積分します。周辺化は確率変数の関数にも使うことができます。それは明らかに確率変数でもあるのですが、期待値\(E_y[f | X = x] = \int\,f(X=x,y)P(X=x,y)\, dy\)を取ることによって使うことができます。積分では確率変数\(X\)が固定されているため、\(E_y[f]\)は\(x\)の関数であることを強調しましたが、\(x\)が固定されている場所（関数の引数）\(f(x,y)\)の期待値を計算することによって除かれます。本質的には、すべての取り得る値\(y\)の平均である\(f(x,y)\)を新しい関数\(E_y[f]\)に置き換えています。ここまでかなり詳細に説明してきましたが、下のコードを見れば直感的に理解できます。もう一つ付け加えるとすれば、値が\(\hat{f}\)の平均値に対する相対的な変化であるということです。余分な項は無視してもかまいませんが、念のため入れておきます。したがって、値の方程式は [vStrumbeljK14]となります。

(12.8)\[\begin{equation} v(x_i) = \int\,f(x_0, x_1, \ldots, x_i,\ldots, x_N)P(x_0, x_1, \ldots, x_i,\ldots, x_N)\, dx_i - E\left[\hat{f}(\vec{x})\right] \end{equation}\]

周辺化\(\int\,f(x_0, x_1, \ldots, x_i,\ldots, x_N)P(x_0, x_1, \ldots, x_i,\ldots, x_N)\, dx_i\)はどのように計算するのでしょうか。既知の確率分布はありません。その場合、データを経験的な分布として考えることで\(P(\vec{x})\)からサンプリングできます。すなわち、データ点をサンプリングすることで\(P(\vec{x})\)からサンプリングできます。ただし、\(\vec{x}\)を共にサンプリングする必要があるため少し複雑です。除かれる特徴量との間に相関がある場合、個々の特徴をランダムに混ぜることはできません。

Strumbeljら[vStrumbeljK14]は\(i\)番目のシャープレイ値を直接推定できることを示しました。

(12.9)\[\begin{equation} \phi_i(\vec{x}) = \frac{1}{M}\sum^M \hat{f}\left(\vec{z}_{+i}\right) - \hat{f}\left(\vec{z}_{-i}\right) \end{equation}\]

\(\vec{z}\)は、実際の例\(\vec{x}\)とランダムな漫然とした例\(\vec{x}'\)から成る「キメラ」の例です。\(\vec{x}\)と\(\vec{x}'\)からランダムに選択し、\(\vec{z}\)を構成します。\(\vec{z}_{+i}\)は例\(\vec{x}\)の\(i\)番目の特徴量を持ち、\(\vec{z}_{-i}\)はランダムな例\(\vec{x}'\)の\(i\)番目の特徴量を持ちます。\(M\)はこの値に対して良いサンプルを得るために十分に大きく選ばれます。[vStrumbeljK14]は\(M\)の選択方法に関する指針を示していますが、基本的には計算可能で妥当な範囲で大きな\(M\)を選択します。この近似の一つの変更点は、期待値（ときには\(\phi_0\)と表記される）を表す明示的な項を使っていることで、「完全性」を有する方程式は次のようになります。

(12.10)\[\begin{equation} \sum_i \phi_i(\vec{x}) = \hat{f}(\vec{x}) - E[\hat{f}(\vec{x})] \end{equation}\]

期待値を\phi_0\(として明示的に含める場合、それは\)\vec{x}$に依存しません。

(12.11)\[\begin{equation} \phi_0 + \sum_{i=1} \phi_i(\vec{x}) = \hat{f}(\vec{x}) \end{equation}\]

特徴量を周辺化することは、特徴量をその平均に置き換えることとは異なります。

この効率的な近似方法、強力な理論、モデル非依存生により、シャープレイ値は予測値に対する特徴量重要度を記述するのに優れた選択肢となります。

12.3. Notebookの実行

このページ上部の    を押すと、このノートブックがGoogle Colab.で開かれます。必要なパッケージのインストール方法については以下を参照してください。

Tip

必要なパッケージをインストールするには、新規セルを作成して次のコードを実行してください。

!pip install dmol-book

もしインストールがうまくいかない場合、パッケージのバージョン不一致が原因である可能性があります。動作確認がとれた最新バージョンの一覧はここから参照できます

import haiku as hk
import jax
import tensorflow as tf
import jax.numpy as jnp
import numpy as np
import matplotlib.pyplot as plt
import urllib
from functools import partial
from jax.example_libraries import optimizers as opt
import dmol

np.random.seed(0)
tf.random.set_seed(0)

ALPHABET = [
    "-",
    "A",
    "R",
    "N",
    "D",
    "C",
    "Q",
    "E",
    "G",
    "H",
    "I",
    "L",
    "K",
    "M",
    "F",
    "P",
    "S",
    "T",
    "W",
    "Y",
    "V",
]

まず、アミノ酸配列とone-hotベクトルを相互に変換する関数を定義しましょう。

def seq2array(seq, L=200):
    return np.pad(list(map(ALPHABET.index, seq)), (0, L - len(seq))).reshape(1, -1)


def array2oh(a):
    a = np.squeeze(a)
    o = np.zeros((len(a), 21))
    o[np.arange(len(a)), a] = 1
    return o.astype(np.float32).reshape(1, -1, 21)


urllib.request.urlretrieve(
    "https://github.com/whitead/dmol-book/raw/master/data/hemolytic.npz",
    "hemolytic.npz",
)
with np.load("hemolytic.npz", "rb") as r:
    pos_data, neg_data = r["positives"], r["negatives"]

12.4. 特徴量重要度の例

ペプチドが赤血球を破壊するかどうか（溶血性）を予測するペプチド予測タスクで、特徴量重要度法の例を見てみましょう。 これはStandard Layersの溶解度予測の例に似ています。データは[BW21]を利用します。 モデルはペプチド配列（例：DDFRD）を取り込み、そのペプチドが溶血性である確率を予測します。 ここでの特徴量重要度法の目標は、どのアミノ酸が溶血活性に最も重要であるかを特定することです。 下の閉じたセルはデータをロードし処理してデータセットにします。

# create labels and stich it all into one
# tensor
labels = np.concatenate(
    (
        np.ones((pos_data.shape[0], 1), dtype=pos_data.dtype),
        np.zeros((neg_data.shape[0], 1), dtype=pos_data.dtype),
    ),
    axis=0,
)
features = np.concatenate((pos_data, neg_data), axis=0)
# we now need to shuffle before creating TF dataset
# so that our train/test/val splits are random
i = np.arange(len(labels))
np.random.shuffle(i)
labels = labels[i]
features = features[i]
L = pos_data.shape[-2]

# need to add token for empty amino acid
# dataset just has all zeros currently
features = np.concatenate((np.zeros((features.shape[0], L, 1)), features), axis=-1)
features[np.sum(features, -1) == 0, 0] = 1.0

batch_size = 16
full_data = tf.data.Dataset.from_tensor_slices((features.astype(np.float32), labels))

# now split into val, test, train
N = pos_data.shape[0] + neg_data.shape[0]
split = int(0.1 * N)
test_data = full_data.take(split).batch(batch_size)
nontest = full_data.skip(split)
val_data, train_data = nontest.take(split).batch(batch_size), nontest.skip(
    split
).shuffle(1000).batch(batch_size)

Jax（Haikuを使用）で畳み込みモデルを再構築し、勾配をもう少し簡単に扱えるようにします。またその他にもいくつかモデルに変更を加えています。畳み込みに加えて、配列の長さとアミノ酸の割合も追加情報として渡しています。

def binary_cross_entropy(logits, y):
    """Binary cross entropy without sigmoid. Works with logits directly"""
    return (
        jnp.clip(logits, 0, None) - logits * y + jnp.log(1 + jnp.exp(-jnp.abs(logits)))
    )


def model_fn(x):
    # get fractions, excluding skip character
    aa_fracs = jnp.mean(x, axis=1)[:, 1:]
    # compute convolutions/poolings
    mask = jnp.sum(x[..., 1:], axis=-1, keepdims=True)
    for kernel, pool in zip([5, 3, 3], [4, 2, 2]):
        x = hk.Conv1D(16, kernel)(x) * mask
        x = jax.nn.tanh(x)
        x = hk.MaxPool(pool, pool, "VALID")(x)
        mask = hk.MaxPool(pool, pool, "VALID")(mask)
    # combine fractions, length, and convolution ouputs
    x = jnp.concatenate((hk.Flatten()(x), aa_fracs, jnp.sum(mask, axis=1)), axis=1)
    # dense layers. no bias, so zeros give P=0.5
    logits = hk.Sequential(
        [
            hk.Linear(256, with_bias=False),
            jax.nn.tanh,
            hk.Linear(64, with_bias=False),
            jax.nn.tanh,
            hk.Linear(1, with_bias=False),
        ]
    )(x)
    return logits


model = hk.without_apply_rng(hk.transform(model_fn))


def loss_fn(params, x, y):
    logits = model.apply(params, x)
    return jnp.mean(binary_cross_entropy(logits, y))


@jax.jit
def hemolytic_prob(params, x):
    logits = model.apply(params, x)
    return jax.nn.sigmoid(jnp.squeeze(logits))


@jax.jit
def accuracy_fn(params, x, y):
    logits = model.apply(params, x)
    return jnp.mean((logits >= 0) * y + (logits < 0) * (1 - y))

rng = jax.random.PRNGKey(0)
xi, yi = features[:batch_size], labels[:batch_size]
params = model.init(rng, xi)

opt_init, opt_update, get_params = opt.adam(1e-2)
opt_state = opt_init(params)


@jax.jit
def update(step, opt_state, x, y):
    value, grads = jax.value_and_grad(loss_fn)(get_params(opt_state), x, y)
    opt_state = opt_update(step, grads, opt_state)
    return value, opt_state

WARNING:absl:No GPU/TPU found, falling back to CPU. (Set TF_CPP_MIN_LOG_LEVEL=0 and rerun for more info.)

epochs = 32
for e in range(epochs):
    avg_v = 0
    for i, (xi, yi) in enumerate(train_data):
        v, opt_state = update(i, opt_state, xi.numpy(), yi.numpy())
        avg_v += v
opt_params = get_params(opt_state)


def predict(x):
    return jnp.squeeze(model.apply(opt_params, x))


def predict_prob(x):
    return hemolytic_prob(opt_params, x)

もし、コードを追うのが大変でも大丈夫です！この章のゴールはモデルの説明を得る方法を示すことであり、必ずしもモデルを構築する方法ではありません。ですから、次の数行に注目してください。ここでは、予測値を得てそれを説明するために、どのようにモデルを使うかを説明します。モデルは、logitsの場合は predict(x) 、確率の場合は predict_prob を介して呼び出されます。

Sequence Models

one-hotエンコーディングと配列モデルについてはDeep Learning on Sequencesを参照してください。

アミノ酸配列、ペプチドを試して、モデルのイメージをつかんでみましょう。モデルはlogits（オッズの対数）を出力し、これをシグモイド関数にかけると確率が得られます。ペプチドは配列からone-hotベクトルの行列に変換する必要があります。ここでは、2つの既知の配列を試してみましょう。Qは溶血生残基としてよく知られており、2番目の配列はポリGで、これは最も単純なアミノ酸です。

s = "QQQQQ"
sm = array2oh(seq2array(s))
p = predict_prob(sm)
print(f"Probability {s} of being hemolytic {p:.2f}")

s = "GGGGG"
sm = array2oh(seq2array(s))
p = predict_prob(sm)
print(f"Probability {s} of being hemolytic {p:.2f}")

Probability QQQQQ of being hemolytic 1.00
Probability GGGGG of being hemolytic 0.00

これは妥当な結果に見えます。モデルの出力は、これらの2つの配列に関する私たちの直感に一致しています。

さて、このモデルの正確度を計算しましょう。非常に良い結果が得られます。

acc = []
for xi, yi in test_data:
    acc.append(accuracy_fn(opt_params, xi.numpy(), yi.numpy()))
print(jnp.mean(np.array(acc)))

0.95208335

12.4.1. 勾配

では、なぜある配列が溶血性となるのか調べてみましょう。まず、入力に対する勾配を計算することから始めます。これは素朴なアプローチでshattered gradientsになりやすいですが、この後計算するintegrated gradientsとsmooth gradientsのプロセスの一部であり、無駄ではありません。それでは溶血性であることが知られている、より複雑なペプチド配列を使って、より興味深い解析をしてみましょう。

def plot_grad(g, s, ax=None):
    # g = np.array(g)
    if ax is None:
        plt.figure()
        ax = plt.gca()
    if len(g.shape) == 3:
        h = g[0, np.arange(len(s)), list(map(ALPHABET.index, s))]
    else:
        h = g
    ax.bar(np.arange(len(s)), height=h)
    ax.set_xticks(range(len(s)))
    ax.set_xticklabels(s)
    ax.set_xlabel("Amino Acid $x_i$")
    ax.set_ylabel(r"Gradient $\frac{\partial \hat{f}(\vec{x})}{\partial x_i}$")

s = "RAGLQFPVGRLLRRLLRRLLR"
sm = array2oh(seq2array(s))
p = predict_prob(sm)
print(f"Probability {s} of being hemolytic {p:.2f}")

Probability RAGLQFPVGRLLRRLLRRLLR of being hemolytic 1.00

コードは至ってシンプルです。それでは勾配を計算してみましょう。

gradient = jax.grad(predict, 0)
g = gradient(sm)
plot_grad(g, s)

モデルはlogitsを出力することを忘れないでください。勾配が正の値であれば、そのアミノ酸が溶血性である確率を高め、負の値であればそのアミノ酸配列が非溶血性である確率を高めることを意味します。興味深いことに、ロイシン(L)とアルギニン(R)は強い位置依存性があることがわかります。

12.4.2. Integrated Gradients

次にintegrated gradients法を実装します。3つの基本的なステップを踏みます。

1. ベースラインから入力ペプチドに向かう入力配列（経路）を作ります。

2. 各入力に対して勾配を計算します。

3. 勾配の合計を計算し、それにベースラインとペプチドの差を乗じます。

ベースラインはすべて0であり、0.5の確率を与えます (logits = 0, a model root)。このベースラインはまさに決定境界上にあります。すべてのグリシンやアラニンのように他のベースラインを使うこともできますが、それらは0.5の確率かそれに近いものであるべきです。ベースライン選択の詳細とインタラクティブな探索は[SLL20]を見てください。

def integrated_gradients(sm, N):
    baseline = jnp.zeros((1, L, 21))
    t = jnp.linspace(0, 1, N).reshape(-1, 1, 1)
    path = baseline * (1 - t) + sm * t

    def get_grad(pi):
        # compute gradient
        # add/remove batch axes
        return gradient(pi[jnp.newaxis, ...])[0]

    gs = jax.vmap(get_grad)(path)
    # sum pieces (Riemann sum), multiply by (x - x')
    ig = jnp.mean(gs, axis=0, keepdims=True) * (sm - baseline)
    return ig


ig = integrated_gradients(sm, 1024)

plot_grad(ig, s)

位置依存性がより顕著になり、アルギニンは位置に対して非常に敏感であることがわかります。先程の通常の勾配と比較すると、質的な変化はあまりありません。

12.4.3. SmoothGrad

SmmotGradを行うステップはintegrated gradientsとほとんど同じです。

1. 入力ペプチドにランダムな摂動を加えた入力配列（経路）を作成します。

2. 各入力の勾配を計算します。

3. 勾配の平均を計算します。

ハイパーパラメータとして、\(\sigma\)があり、原理的にはモデルの出力を変化させながら、できるだけ小さくする必要があります。

def smooth_gradients(sm, N, rng, sigma=1e-3):
    baseline = jnp.zeros((1, L, 21))
    t = jax.random.normal(rng, shape=(N, sm.shape[1], sm.shape[2])) * sigma
    path = sm + t
    # remove examples that are negative and force summing to 1
    path = jnp.clip(path, 0, 1)
    path /= jnp.sum(path, axis=2, keepdims=True)

    def get_grad(pi):
        # compute gradient
        # add/remove batch axes
        return gradient(pi[jnp.newaxis, ...])[0]

    gs = jax.vmap(get_grad)(path)
    # mean
    ig = jnp.mean(gs, axis=0, keepdims=True)
    return ig


sg = smooth_gradients(sm, 1024, jax.random.PRNGKey(0))
plot_grad(sg, s)

通常の勾配の結果に酷似しているように見えます。これは、1次元の入力と浅いネットワークが、shattered gradientsに対してそれほど敏感ではないためと思われます。

12.4.4. シャープレイ値

次に、式10.9を使って各特徴量に対するシャープレイ値を近似してみましょう。シャープレイ値の計算は勾配を必要としないため、これまでのアプローチとは異なります。基本的なアルゴリズムは次のようになります。

1. ランダムな点x’を選択します。

2. xとx’を組み合わせて点\(z\)を作ります。

3. 予測関数の変化を計算します。

効率化のために行った1つの工夫は、パディングで配列を変更しないようにすることです。基本的に配列を長くするようなことはしないようにしています。

def shapley(i, sm, sampled_x, rng, model):
    M, F, *_ = sampled_x.shape
    z_choice = jax.random.bernoulli(rng, shape=(M, F))
    # only swap out features within length of sm
    mask = jnp.sum(sm[..., 1:], -1)
    z_choice *= mask
    z_choice = 1 - z_choice
    # construct with and w/o ith feature
    z_choice = z_choice.at[:, i].set(0.0)
    z_choice_i = z_choice.at[:, i].set(1.0)
    # select them via multiplication
    z = sm * z_choice[..., jnp.newaxis] + sampled_x * (1 - z_choice[..., jnp.newaxis])
    z_i = sm * z_choice_i[..., jnp.newaxis] + sampled_x * (
        1 - z_choice_i[..., jnp.newaxis]
    )
    v = model(z_i) - model(z)
    return jnp.squeeze(jnp.mean(v, axis=0))


# assume data is alrady shuffled, so just take M
M = 4096
sl = len(s)
sampled_x = train_data.unbatch().batch(M).as_numpy_iterator().next()[0]
# make batched shapley so we can compute for all features
bshapley = jax.vmap(shapley, in_axes=(0, None, None, 0, None))
sv = bshapley(
    jnp.arange(sl),
    sm,
    sampled_x,
    jax.random.split(jax.random.PRNGKey(0), sl),
    predict,
)

# compute global expectation
eyhat = 0
for xi, yi in full_data.batch(M).as_numpy_iterator():
    eyhat += jnp.mean(predict(xi))
eyhat /= len(full_data)

シャープレイ値に関する一つの良いチェックは、それらの合計がモデル関数の値からすべてのインスタンスにわたる期待値を引いたものに等しいことを確認することです。ただし [vStrumbeljK14]の式を近似して使用しているので、完全な一致は期待できません。この値は次のように計算されます。

print(np.sum(sv), predict(sm))

6.7373457 8.068422

予想通り、いくらか違います。これは今回使っている近似法の影響です。サンプル数がシャープレイ値の合計にどのように影響するかを調べることで、それを確認することができます。

Fig. 12.1 シャープレイ値近似における、シャープレイ値の総和とサンプル数の関数の間数値の比較

徐々に収束しています。最後に、個々のシャープレイ値を表示して見ましょう。それが予測に対する説明となります。

plot_grad(sv, s)

ここまでに見てきた4つの手法、勾配法、Integrated Gradient法、SmoothGrad法、シャープレイ値の結果を並べて示します。

heights = []
plt.figure(figsize=(12, 4))
x = np.arange(len(s))
for i, (gi, l) in enumerate(zip([g, ig, sg], ["Gradient", "Integrated", "Smooth"])):
    h = gi[0, np.arange(len(s)), list(map(ALPHABET.index, s))]
    plt.bar(x + i / 5 - 1 / 4, h, width=1 / 5, edgecolor="black", label=l)
plt.bar(x + 3 / 5 - 1 / 4, sv, width=1 / 5, edgecolor="black", label="Shapley")
ax = plt.gca()
ax.set_xticks(range(len(s)))
ax.set_xticklabels(s)
ax.set_xlabel("Amino Acid $x_i$")
ax.set_ylabel(r"Importance [logits]")
plt.legend()
plt.show()

普段からペプチドを扱っている者として、ここではシャープレイ値が最も正確だと思います。LとRのパターンが重要だとは考えていませんでしたが、シャープレイ値はそう示しています。また他の手法の結果と異なり、シャープレイ値はフェニルアラニン(F)が重要な効果を持つとは示していません。

この結果から何を結論づけることができるでしょうか。おそらく次のような説明を加えることができるでしょう。「この配列は、主にグルタミン、プロリン、そしてロイシンとアルギニンの配列によって溶血性であると予測されています」。

12.5. 特徴量重要度は何のためにあるのか？

特徴量重要度は、実用的な予測や洞察を与える明確な説明につながることはほとんどありません。因果関係がないため、存在しない特徴量の説明に意味を見出すことに繋がりかねません[CK18]。もう一つの注意点は、実際の化学物質の体系ではなく、モデルを説明しているということです。例えば、「溶血活性は5番目のグルタミンによるものです」と解釈するのは避けましょう。代わりに「モデルは5番目にグルタミンが位置するため溶血活性であると予測しました」としてください。

実用的な説明は、特徴量をどのように変えれば結果に影響するかを示すもので、結果の原因を知っていることに似ています。したがって、上述の理由から、特徴量重要度に説明性があるかどうかについては議論が続いています[Lip18]。参考までに、特徴量重要度を人の概念に結びつけようとする研究分野は、testing with concept activation vectors（TCAV）[KWG+18]と呼ばれています。ちなみに私自身はXAIのために特徴量重要度をあまり使っていません。それは、説明が実用的でも因果関係を示すものでもなく、しばしば他の混乱を招くからです。

12.6. 学習データの重要度

もう一つの私たちが期待する説明や解釈は、どの学習データ点が予測に最も貢献しているかということです。これは次の質問に対するより直接的な回答になります。「なぜ私のモデルはこれを予測したのでしょうか」。ニューラルネットワークは学習データの結果であり、なぜその予測がなされたのかに対する答えは学習データを辿ることで得られます。ある予測に対して学習データをランク付けすることで、どの学習データ点がニューラルネットワークの予測に影響を与えているのかに関する洞察を得ることができます。これは影響関数\(\mathcal{I}(x_i, x)\)のようであり、学習データ点\(i\)と入力\(x\)に対する影響度スコアを与えます。影響度を計算する最も簡単な方法は、ニューラルネットワークに\(x_i\)がある場合（つまり\(\hat{f}(x)\)）とない場合（つまり\(\hat{f}_{-x_i}(x)\)）を学習して、影響度を以下のように定義します。

(12.12)\[\begin{equation} \mathcal{I}(x_i, x) = \hat{f}_{-x_i}(x) - \hat{f}(x) \end{equation}\]

例えば、学習データから学習データ点\(x_i\)を除いた後、予測値が高くなれば、その点は正の影響力を持っているということになります。この影響関数の計算は通常データ点の数だけモデルを繰り返し学習する必要がありますが、通常は計算できません。 [KL17] show a way to approximate this by looking at infinitesimal changes to the weights of each training point. これらの影響関数を計算するには損失関数に関するHessianを計算する必要があるため、一般的には使用されません。しかし、JAXを使っている場合は、その計算を簡単に行うことができます。

カーネルモデルを利用する場合、学習データが特徴量となります。intergrated gradientsのような上記の方法は、学習データの重要度を与えます。

学習データの重要度はディープラーニングの専門家にとって有用な解釈を提供します。ある予測に対してどの学習データ点が最も影響力を持っているのかを教えてくれます。これはデータに関する問題に対処する場合や偽陽性に対する説明を辿るのに役立ちます。しかし、ディープラーニングモデルの予測結果を利用する一般の利用者は、おそらく学習データのランク付けだけでは満足しないでしょう。

12.7. サロゲートモデル

解釈可能性におけるより一般的な考え方の一つは、解釈可能なモデルをブラックボックスモデルに特定の例の近傍で適合させることでしょう。なぜなら、解釈可能なモデルはたいてい大域的にブラックボックスモデルに適合させることはできないからです。そうでなければ、最初から解釈可能なモデルを使い、ブラックボックスモデルは使わないでしょう。しかし、解釈可能なモデルは興味ある例の周辺の小さな領域にだけ当てはめれば、局所的に正しい解釈可能なモデルを使って説明を与えることができます。この解釈可能なモデルをローカルサロゲートモデルと呼びます。解釈可能なローカルサロゲートモデルには、決定木、線形モデル、（簡潔な説明のための）スパース線形モデル、ナイーブベイズ分類器などがあります。

ローカルサロゲートモデルとして一般的に知られているアルゴリズムはLocal Interpretable Model-Agnostic Explanations (LIME) [RSG16a]と呼ばれています。LIMEは、元のブラックボックスモデルを学習させた損失関数を利用して、ローカルサロゲートモデルを興味ある例の近傍にフィットさせます。ローカルサロゲートモデルの損失関数は、サロゲートモデルを回帰する際に、興味ある例に近い点を評価するよう重み付けされます。LIMEの論文では、サロゲートモデルのスパース化を表記に含めていますが、それはローカルサロゲートモデルの特性ではないため、ここでは一旦省きます。よって、サロゲートモデルの損失は次のように定義されます。

(12.13)\[\begin{equation} \mathcal{l^s}\left(x'\right) = w(x', x)\mathcal{l}\left(\hat{f}_s(x'), \hat{f}(x')\right) \end{equation}\]

\(w(x', x)\)は興味ある例\(x\)の近くにある点に重みを付ける重みカーネル関数、 \(\mathcal{l}(\cdot,\cdot)\)は元のブラックボックスモデルの損失、\(\hat{f}(\cdot)\) はブラックボックスモデル、\(\hat{f}_s(\cdot)\)はサロゲートモデルを表します。

[RSG16a]で定式化されているLIMEは特徴量の重要度に関する記述を与えますが、サロゲートモデルによっては解釈可能な場合もあります。例えば決定木のようなものです。

重み関数は少しアドホックです。つまりデータ型に依存します。スカラーラベルの回帰タスクでは、カーネル関数を使いますが、様々な選択肢があります。ガウシアン、コサイン、エパネチコフなどです。テキストデータでは、LIMEの実装では、ハミング距離を使っています。これは単に2つの文字列の間で一致しないテキストトークンの数をカウントするものです。画像もハミング距離を使いますが、スーパーピクセルは例と同じか空白とします。

点 \(x'\)はどのように生成されるのでしょうか。連続値の場合、\(x'\)は一様にサンプリングされますが、特徴空間はしばしば閉じていないため、これは非常に難しいことです。重み付き関数に従って\(x'\)をサンプリングし、重み付けを省略すれば、（それに従ってサンプリングされたので）閉じていない特徴空間のような問題を避けることができます。一般に、連続ベクトル特徴空間では、LIMEは少し主観的です。画像やテキストの場合、\(x'\)はトークン（単語）をマスキングする、スーパーピクセルをゼロ化（黒化）することにより形成されます。これは、シャープレイ値にかなり近い説明となるはずで、実際、LIMEがシャープレイ値と同等であることを、いくつかの小さな表記法の変更で示すことができます。

12.8. 反実仮想

反実仮想

最適化の定式化はXAIで使われているものですが、他の文脈の反実仮想は”近さ”の基準を持っていません。

反実仮想は最適化問題の解です。\(x\)と異なるラベルを持ち、\(x\)にできるだけ近い例\(x'\)を見つけます[WMR17]。これは次のように定式化できます。

(12.14)\[\begin{split}\textrm{minimize}\qquad d(x, x')\\ \textrm{such that}\qquad \hat{f}(x) \neq \hat{f}(x')\end{split}\]

\(\hat{f}(x)\)がスカラーを出力する回帰問題では、制約条件を\(\hat{f}(x)\)からある\(\Delta\)だけ離すように修正する必要があります。この最適化問題を満たす\(x'\) は反実仮想（発生しなかった条件、異なる結果を導いたであろう条件）と呼ばれます。通常、\(x'\)を求めることは、無微分最適化として扱われます。\(\frac{\partial \hat{f}}{\partial x'}\)を計算して制約付き最適化しますが、実際にはモンテカルロ最適化のように \(\hat{f}(x) \neq \hat{f}(x')\)までランダムに\(x\)を摂動させた方が速い場合があります。教師なし学習で新しい\(x'\)を提案できる生成モデルを使用することもできます。分子に関する普遍的な反実仮想生成器については[WSW22]を参照してください。分子のグラフニューラルネットワークに特化した手法については[NB20]を参照してください。

距離の定義は、LIMEの説明の中でも述べたように、重要な主観的関心ごとです。分子構造の文脈で使われる一般的な距離の例は、Moragnフィンガープリント[RH10]のような分子フィンガープリント/記述子のタニモト係数（またはJaccard係数）です。

反実仮想はシャープレイ値と比較して一つ欠点があります。それは完全な説明を与えてはくれないことです。シャープレイ値は予測値の合計であり、説明のどのような部分も見逃していないことを意味しています。一方、反実仮想はできるだけ少ない特徴量を変える（距離を最小化する）ため、予測に寄与している一部の特徴量についての情報を見逃してしまうことがあります。またシャープレイ値の利点は実用的であることですが、反実仮想は直接使用することができます。

12.8.1. 例

上記のペプチドの例でこのアイデアを素早く実装することができます。距離はハミング距離と定義します。そして\(x'\)は一つのアミノ酸置換です。これを列挙してラベルの置換ができるかどうか試してみましょう。まず1回の置換を行う関数を定義します。

def check_cf(x, i, j):
    # copy
    x = jnp.array(x)
    # substitute
    x = x.at[:, i].set(0)
    x = x.at[:, i, j].set(1)
    return predict(x)


check_cf(sm, 0, 0)

DeviceArray(8.552943, dtype=float32)

次に、jnp.meshgridで可能なすべての置換を作り、vmapで先ほど定義した関数を適応します。ravel()<jax.numpy.ravel>`はインデックスの配列を一次元にするため、複雑なvmapを行う必要はありません。

ii, jj = jnp.meshgrid(jnp.arange(sl), jnp.arange(21))
ii, jj = ii.ravel(), jj.ravel()
x = jax.vmap(check_cf, in_axes=(None, 0, 0))(sm, ii, jj)

次に、予測値が負になった（すなわちlogitsが0より小さい）アミノ酸置換をすべて表示します。

from IPython.core.display import display, HTML

out = ["<tt>"]
for i, j in zip(ii[jnp.squeeze(x) < 0], jj[jnp.squeeze(x) < 0]):
    out.append(f'{s[:i]}<span style="color:red;">{ALPHABET[j]}</span>{s[i+1:]}<br/>')
out.append("</tt>")
display(HTML("".join(out)))

RAGL-FPVGRLLRRLLRRLLR
RAGLQF-VGRLLRRLLRRLLR
RAGLAFPVGRLLRRLLRRLLR
RAGLQFAVGRLLRRLLRRLLR
RAGLCFPVGRLLRRLLRRLLR
RAGLQFCVGRLLRRLLRRLLR
RAGLQFPCGRLLRRLLRRLLR
RAGLIFPVGRLLRRLLRRLLR
RAGLQFIVGRLLRRLLRRLLR
RAGLLFPVGRLLRRLLRRLLR
RAGLQFLVGRLLRRLLRRLLR
RAGLFFPVGRLLRRLLRRLLR
RAGLQFFVGRLLRRLLRRLLR
RAGLQFPFGRLLRRLLRRLLR
RAGLPFPVGRLLRRLLRRLLR
RAGLTFPVGRLLRRLLRRLLR
RAGLWFPVGRLLRRLLRRLLR
RAGLVFPVGRLLRRLLRRLLR
RAGLQFVVGRLLRRLLRRLLR

解釈はいくつかありますが、基本的にはグルタミンを疏水基と交換するか、プロリンをV、F、A、Cに置き換えることでペプチドを非溶血性にするという解釈です。反実仮想として述べると、「もしグルタミンを疎水性アミノ酸に交換すれば、そのペプチドは非溶血性になるでしょう」ということになります。

12.9. 特定のアーキテクチャの説明

上記と同じ原則がGNNにも適用されますが、これらのアイデアをグラフ上で動作するように変換する最適な方法については様々なアイデアがあります。GNNに特化した解釈可能性の理論については[AZL21]を、GNNで説明を構築するために利用できる手法については[YYGJ20]を参照してください。

NLPは説明と解釈を構築するための特別なアプローチが存在するもう一つの分野です。最近の調査として[MRC21]を参照してください。

12.10. モデル非依存的な分子の反実仮想の説明

化学における反実仮想に関連する主な課題は(12.14)の微分を計算することの難しさです。したがって、このタスクに焦点を当てたほとんどの手法は、これまで見てきたようにモデルのアーキテクチャに特化しています。Wellawatteら[WSW22]はモデルのアーキテクチャに関係なく分子に対してこれを行うMolecular Model Agnostic Counterfactual Explanations（MMACE）という方法を導入しています。

MMACE法はexmolパッケージで実装されています。分子とモデルを与えると、exmolは局所的な反実仮想的説明を生成することができます。MMACE法には2つの主要なステップがあります。まず、与えられた基本分子を中心に局所的な化学空間を展開します。次に、各サンプル点に、ユーザが指定したモデルアーキテクチャのラベルを付けます。これらのラベルは、局所的な化学空間における反実仮想を特定するために使用されます。MMACE法はモデル非依存的で、exmolパッケージは分類と回帰の両方のタスクに対して反実仮想を生成することができます。

それでは、exmolを使ってどのように反実仮想を生成するのか見てみましょう。この例では、分子の臨床毒性を予測するランダムフォレストモデルを学習します。この二値分類タスクでは、MoleculeNetグループ[WRF+18]が発表したClassification章で使用したのと同じデータセットを使います。

12.11. Notebookの実行

上の　   をクリックしてインタラクティブなGoogle Colabでこのページを開始しましょう。ご自身の環境でもGoogle Colabでも、パッケージのインストール関する詳細は以下を参照してください。

Tip

パッケージをインストールするには、新規セルを作成して次のコードを実行してください。

!pip install exmol jupyter-book matplotlib numpy pandas seaborn sklearn mordred[full] rdkit-pypi

import pandas as pd
import matplotlib as mpl
import matplotlib.pyplot as plt
import rdkit, rdkit.Chem, rdkit.Chem.Draw
from rdkit.Chem.Draw import IPythonConsole
import numpy as np
import mordred, mordred.descriptors
import warnings
from sklearn.model_selection import train_test_split

import tensorflow as tf
from tensorflow import keras
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import Dense
from tensorflow.keras.wrappers.scikit_learn import KerasClassifier
from sklearn.model_selection import StratifiedKFold
from sklearn.model_selection import cross_val_score
import exmol

IPythonConsole.ipython_useSVG = True


toxdata = pd.read_csv(
    "https://github.com/whitead/dmol-book/raw/master/data/clintox.csv.gz"
)

# make object that can compute descriptors
calc = mordred.Calculator(mordred.descriptors, ignore_3D=True)
# make subsample from pandas df
molecules = [rdkit.Chem.MolFromSmiles(smi) for smi in toxdata.smiles]

# view one molecule to make sure things look good.
molecules[0]

データをインポートしたら、Mordredパッケージで入力記述子を生成します。

# Get valid molecules from the sample
valid_mol_idx = [bool(m) for m in molecules]
valid_mols = [m for m in molecules if m]
# Compute molecular descriptors using Mordred
features = calc.pandas(valid_mols, quiet=True)
labels = toxdata[valid_mol_idx].FDA_APPROVED
# Standardize the features
features -= features.mean()
features /= features.std()

# we have some nans in features, likely because std was 0
features = features.values.astype(float)
features_select = np.all(np.isfinite(features), axis=0)
features = features[:, features_select]
print(f"We have {len(features)} features per molecule")

We have 1478 features per molecule

この例では、Kerasで実装されたシンプルで密なニューラルネットワーク分類器を使用します。まず、このシンプルな分類器を学習し、それを使ってexmolの反実仮想のラベルを生成します。学習済みモデルの性能を改善することで、より正確な結果を期待することができますが、exmolの仕組みを理解するには、今のところ以下の例で十分です。

# Train and test spit
X_train, X_test, y_train, y_test = train_test_split(
    features, labels, test_size=0.2, shuffle=True
)
ft_shape = X_train.shape[-1]

# reshape data
X_train = X_train.reshape(-1, ft_shape)
X_test = X_test.reshape(-1, ft_shape)

それではモデルを構築して実行してみましょう。ディープラーニングの概要 章に密なモデルに関する詳しいイントロダクションがあります。

model = tf.keras.models.Sequential()
model.add(tf.keras.Input(shape=(ft_shape,)))
model.add(tf.keras.layers.Dense(32, activation="relu"))
model.add(tf.keras.layers.Dense(32))
model.add(Dense(1, activation="sigmoid"))
model.compile(loss="binary_crossentropy", optimizer="adam", metrics=["accuracy"])

# Model training
model.fit(X_train, y_train, epochs=50, batch_size=32, verbose=0)
_, accuracy = model.evaluate(X_test, y_test)
print(f"Model accuracy: {accuracy*100:.2f}%")

 1/10 [==>...........................] - ETA: 1s - loss: 0.1194 - accuracy: 0.9688


10/10 [==============================] - 0s 1ms/step - loss: 0.2520 - accuracy: 0.9392

Model accuracy: 93.92%

私たちが作ったモデルの正確度は良さそうですね！

次に、SMILES及び/またはSEFLIESの分子表現を取り込み、学習済み分類器からラベルの予測を出力するラッパー関数を書きます。SELFIESの詳しい説明はDeep Learning on Sequencesの章にあります。このラッパー関数は exmolのexmol.sample_spaceに入力として与えられ、与えられたベースとなる分子の周りに局所的な化学空間を作ります。exmolは、Superfast Traversal, Optimization, Novelty, Exploration and Discovery (STONED)アルゴリズム[NPK+21]を生成アルゴリズムとして使用して、局所空間を拡張していきます。ベースとなる分子が与えられると、STONEDアルゴリズムは分子のSELFIES表現をランダムに変異させます。これらの変異は文字列置換、挿入、欠損です。

def model_eval(smiles, selfies):
    molecules = [rdkit.Chem.MolFromSmiles(smi) for smi in smiles]
    features = calc.pandas(molecules)
    features = features.values.astype(float)
    features = features[:, features_select]
    labels = np.round(model.predict(np.nan_to_num(features).reshape(-1, ft_shape)))
    return labels

次に、STONEDを使って、exmol.sample_spaceで局所的な化学空間をサンプリングしてみます。この例では、引数num_samplesでサンプル空間の大きさを変更します。ここで選択したベースとなる分子はFDA非承認分子です。

space = exmol.sample_space("C1CC(=O)NC(=O)C1N2CC3=C(C2=O)C=CC=C3N", model_eval);

いったんサンプル空間を作成したら、exmol.sample_space関数を使って局所的な化学空間の反実仮想を特定できます。各反実仮想は、付加情報を含むpythonのdataclassです。

exps = exmol.cf_explain(space, 2)
exps[1]

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Cell In [32], line 2
      1 exps = exmol.cf_explain(space, 2)
----> 2 exps[1]

IndexError: list index out of range

生成された反実仮想はexmolのプロットコードexmol.sample_spaceとexmol.sample_spaceを使って簡単に可視化できます。ベースと反実仮想の間の類似度はECFP4フィンガープリントのタニモト係数です。上位3つの反実仮想をここに示します。

exmol.plot_cf(exps, nrows=1)

ここで選択したベースとなる分子はFDA非承認です。生成された反実仮想を見ると、複素環式基は毒性に影響を与えると結論づけることができます。したがって、我々のモデルによると、複素環式基を変更することでベースとなる分子は非毒性化されるかもしれません。このことは、反実仮想の説明がどのように修正を加えることができるかについての実用的な洞察を与える理由も示しています。

最後に、生成した化学空間も可視化してみましょう！

exmol.plot_space(space, exps)

12.12. まとめ

• ディープラーニングモデルの解釈は、モデルの正確性を保証し、予測を人にとって有用なものにするために必要不可欠です。法令順守のために要求されることもあります。

• ニューラルネットワークの解釈可能性は、より広範なトピックであるAIにおける説明可能性（XAI）の一部であり、このトピックはまだ初期段階です。

• 説明はまだ定義が曖昧ですが、多くの場合、モデルの特徴量で表現されます。

• 説明の戦略としては、特徴量重要度、学習データの重要度、反実仮想、局所的に正確なサロゲートモデルなどがあります。

• ほとんどの説明は例ごとに（推論時に）生成されます。

• 最も体系的ですが、計算コストのかかる説明はシャープレイ値です。

• 反実仮想は最も直感的で満足のいく説明を提供しますが、完全な説明にはならないかもしれないという意見があります。

• exmolはモデル非依存的な分子の反実仮想の説明を生成するソフトウェアです。

12.13. Cited References

WRF+18

Zhenqin Wu, Bharath Ramsundar, Evan N Feinberg, Joseph Gomes, Caleb Geniesse, Aneesh S Pappu, Karl Leswing, and Vijay Pande. Moleculenet: a benchmark for molecular machine learning. Chemical science, 9(2):513–530, 2018.

LS04

John D Lee and Katrina A See. Trust in automation: designing for appropriate reliance. Human factors, 46(1):50–80, 2004.

DVK17(1,2)

Finale Doshi-Velez and Been Kim. Towards a rigorous science of interpretable machine learning. arXiv preprint arXiv:1702.08608, 2017.

GF17

Bryce Goodman and Seth Flaxman. European Union regulations on algorithmic decision-making and a “right to explanation”. AI Magazine, 38(3):50–57, 2017.

Dev19

Organisation for Economic Co-operation and Development. Recommendation of the Council on Artificial Intelligence. 2019. URL: https://legalinstruments.oecd.org/en/instruments/OECD-LEGAL-0449.

CLG+15

Rich Caruana, Yin Lou, Johannes Gehrke, Paul Koch, Marc Sturm, and Noemie Elhadad. Intelligible models for healthcare: Predicting pneumonia risk and hospital 30-day readmission. In Proceedings of the 21th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, 1721–1730. ACM, 2015.

Mil19

Tim Miller. Explanation in artificial intelligence: insights from the social sciences. Artificial intelligence, 267:1–38, 2019.

MSK+19

James W Murdoch, Chandan Singh, Karl Kumbier, Reza Abbasi-Asl, and Bin Yu. Interpretable machine learning: definitions, methods, and applications. eprint arXiv, pages 1–11, 2019. URL: http://arxiv.org/abs/1901.04592.

MSMuller18

Grégoire Montavon, Wojciech Samek, and Klaus-Robert Müller. Methods for interpreting and understanding deep neural networks. Digital Signal Processing, 73:1–15, 2018.

BCB14

Dzmitry Bahdanau, Kyunghyun Cho, and Yoshua Bengio. Neural machine translation by jointly learning to align and translate. arXiv preprint arXiv:1409.0473, 2014.

AGFW21

Mehrad Ansari, Heta A Gandhi, David G Foster, and Andrew D White. Iterative symbolic regression for learning transport equations. arXiv preprint arXiv:2108.03293, 2021.

BD00

Lynne Billard and Edwin Diday. Regression analysis for interval-valued data. In Data analysis, classification, and related methods, pages 369–374. Springer, 2000.

UT20

Silviu-Marian Udrescu and Max Tegmark. Ai feynman: a physics-inspired method for symbolic regression. Science Advances, 6(16):eaay2631, 2020.

CSGB+20

Miles Cranmer, Alvaro Sanchez Gonzalez, Peter Battaglia, Rui Xu, Kyle Cranmer, David Spergel, and Shirley Ho. Discovering symbolic models from deep learning with inductive biases. Advances in Neural Information Processing Systems, 33:17429–17442, 2020.

WSW22(1,2,3)

Geemi P Wellawatte, Aditi Seshadri, and Andrew D White. Model agnostic generation of counterfactual explanations for molecules. Chem. Sci., pages –, 2022. URL: http://dx.doi.org/10.1039/D1SC05259D, doi:10.1039/D1SC05259D.

RSG16a(1,2,3)

Marco Tulio Ribeiro, Sameer Singh, and Carlos Guestrin. " why should i trust you?" explaining the predictions of any classifier. In Proceedings of the 22nd ACM SIGKDD international conference on knowledge discovery and data mining, 1135–1144. 2016.

RSG16b

Marco Tulio Ribeiro, Sameer Singh, and Carlos Guestrin. Model-agnostic interpretability of machine learning. arXiv preprint arXiv:1606.05386, 2016.

WMR17(1,2,3)

Sandra Wachter, Brent Mittelstadt, and Chris Russell. Counterfactual explanations without opening the black box: automated decisions and the gdpr. Harv. JL & Tech., 31:841, 2017.

KL17(1,2)

Pang Wei Koh and Percy Liang. Understanding black-box predictions via influence functions. In International Conference on Machine Learning, 1885–1894. PMLR, 2017.

SML+21

Wojciech Samek, Grégoire Montavon, Sebastian Lapuschkin, Christopher J. Anders, and Klaus-Robert Müller. Explaining deep neural networks and beyond: a review of methods and applications. Proceedings of the IEEE, 109(3):247–278, 2021. doi:10.1109/JPROC.2021.3060483.

BFL+17

David Balduzzi, Marcus Frean, Lennox Leary, J. P. Lewis, Kurt Wan-Duo Ma, and Brian McWilliams. The shattered gradients problem: if resnets are the answer, then what is the question? In Doina Precup and Yee Whye Teh, editors, Proceedings of the 34th International Conference on Machine Learning, volume 70 of Proceedings of Machine Learning Research, 342–350. PMLR, 06–11 Aug 2017. URL: http://proceedings.mlr.press/v70/balduzzi17b.html.

STY17(1,2)

Mukund Sundararajan, Ankur Taly, and Qiqi Yan. Axiomatic attribution for deep networks. In International Conference on Machine Learning, 3319–3328. PMLR, 2017.

STK+17(1,2)

Daniel Smilkov, Nikhil Thorat, Been Kim, Fernanda Viégas, and Martin Wattenberg. Smoothgrad: removing noise by adding noise. arXiv preprint arXiv:1706.03825, 2017.

MBL+19

Grégoire Montavon, Alexander Binder, Sebastian Lapuschkin, Wojciech Samek, and Klaus-Robert Müller. Layer-Wise Relevance Propagation: An Overview, pages 193–209. Springer International Publishing, Cham, 2019. URL: https://link.springer.com/chapter/10.1007%2F978-3-030-28954-6_10.

vStrumbeljK14(1,2,3,4)

Erik Štrumbelj and Igor Kononenko. Explaining prediction models and individual predictions with feature contributions. Knowledge and information systems, 41(3):647–665, 2014.

BW21

Rainier Barrett and Andrew D. White. Investigating active learning and meta-learning for iterative peptide design. Journal of Chemical Information and Modeling, 61(1):95–105, 2021. URL: https://doi.org/10.1021/acs.jcim.0c00946, doi:10.1021/acs.jcim.0c00946.

SLL20

Pascal Sturmfels, Scott Lundberg, and Su-In Lee. Visualizing the impact of feature attribution baselines. Distill, 2020. https://distill.pub/2020/attribution-baselines. doi:10.23915/distill.00022.

CK18

Kangway V Chuang and Michael J Keiser. Comment on “predicting reaction performance in c–n cross-coupling using machine learning”. Science, 362(6416):eaat8603, 2018.

Lip18

Zachary C Lipton. The mythos of model interpretability: in machine learning, the concept of interpretability is both important and slippery. Queue, 16(3):31–57, 2018.

KWG+18

Been Kim, Martin Wattenberg, Justin Gilmer, Carrie Cai, James Wexler, Fernanda Viegas, and others. Interpretability beyond feature attribution: quantitative testing with concept activation vectors (tcav). In International conference on machine learning, 2668–2677. PMLR, 2018.

NB20

Danilo Numeroso and Davide Bacciu. Explaining deep graph networks with molecular counterfactuals. arXiv preprint arXiv:2011.05134, 2020.

RH10

David Rogers and Mathew Hahn. Extended-connectivity fingerprints. Journal of chemical information and modeling, 50(5):742–754, 2010.

AZL21

Chirag Agarwal, Marinka Zitnik, and Himabindu Lakkaraju. Towards a rigorous theoretical analysis and evaluation of gnn explanations. arXiv preprint arXiv:2106.09078, 2021.

YYGJ20

Hao Yuan, Haiyang Yu, Shurui Gui, and Shuiwang Ji. Explainability in graph neural networks: a taxonomic survey. arXiv preprint arXiv:2012.15445, 2020.

MRC21

Andreas Madsen, Siva Reddy, and Sarath Chandar. Post-hoc interpretability for neural nlp: a survey. arXiv preprint arXiv:2108.04840, 2021.

NPK+21

AkshatKumar Nigam, Robert Pollice, Mario Krenn, Gabriel dos Passos Gomes, and Alán Aspuru-Guzik. Beyond generative models: superfast traversal, optimization, novelty, exploration and discovery (stoned) algorithm for molecules using selfies. Chem. Sci., 12:7079–7090, 2021. URL: http://dx.doi.org/10.1039/D1SC00231G, doi:10.1039/D1SC00231G.

previous

11. Modern Molecular NNs

next

13. アテンションレイヤー