プラスマイナスの影響を無視するために絶対値をとっています。この変数重要度を簡単に可視化するための関数として、summry_plot()が用意されています。Valueの意味で、LSTATとRMが非常に重要な変数であることが見て取れます。Valueがそのまま打たれます。X)上に来るほど先程の棒グラフの意味で重要な変数になります。Valueになります。Valueの分布に加えて、LSTATは低いほうがプラスの要因に、RMは高いほうがプラスの要因になりそうなことが見て取れます。Plotのような可視化も可能となります。dependence_plot()を使います。今回は特に重要度の高かったLSTATを見てみましょう。Valueが小さくなることが見て取れます。色付けは交互作用が見れるように他の変数の値でされています。デフォルトでは交互作用が一番はっきり現れる変数が自動で選ばれます。今回はDISが選ばれていて、DISの値が大きいほど赤く、小さいほど青い点が打たれます。centreからの距離を表しています。Debug可視化アルゴリズムSHAPを理解するために、ゲーム理論を調べてまとめました。データ分析関連のまとめ*1:この数値例は岡田章『新板ゲーム理論』の例をそのままお借りしています。*2:たとえば、山分けで3等分も考えられますが、この設定ではそれはうまくいきません。A君B君C君の3人で働いて3等分すると報酬は8万円になります。この場合、A君B君だけで働いて報酬を2等分するとA君とB君は10万円を稼ぐことができ、C君を外すインセンティブが生まれてしまいます。Valueは数学的に証明された望ましい性質をいくつかもっていますが、ここでは具体例と直感的な性質のみ説明しています。詳細を知りたい場合は協力ゲーム理論の教科書か論文をご確認下さい。*4:僕はこれの正式な名称を知りません。Plotまとめ参考文献はじめに本記事では、tidymodelsを用いて機械学習モデルを作成し、それをDALEXを用いて解釈する方法をまとめています。Modelsを目的とした複数のパッケージから成り立っています。具体的にどんなパッケージが所属しているかはリンクを確認して下さい。本記事では、特にメインパッケージであるDALEXとingredientsを用います。Effect(ALE)など、特徴量とターゲット変数の関係を見るための関数が用意されています。本記事ではtidymodelsの使い方には触れないので、過去の記事を参考にして頂ければと思います。dropout009.hatenablog.comdropout009.hatenablog.comdropout009.hatenablog.comまた、本記事では解釈手法そのものには詳しく触れません。ICEの詳細な説明は以下の記事も参考にして頂ければと思います。dropout009.hatenablog.comパッケージ本記事で用いるパッケージは以下になります。0)))}シミュレーション1データPDPとICEの関係を浮き彫りにするため、今回はシミュレーションデータを用意します。関係してないという状況を考えます。できるだけ簡単なほうが結果がわかりやすいと思うので、以下の線形で加法な形に特定します。平均0で標準偏差0.1の正規分布に従うノイズが乗るとしています。シミュレーションのために乱数を発生させます。Rには確率分布を扱うための関数がデフォルトで用意されていますが、distributions3を用いると、より直感的に確率分布を扱うことができます。関係を確認しておきましょう。無相関であることが見て取れます。モデルデータの準備ができたので、モデルを作成します。parsnipを使って作成します。Forestを使うことにします。Forestを利用することと、そのハイパーパラメータを指定set_engin()で利用パッケージとパッケージ固有のパラメータを指定fit()でデータを指定して学習までを一気に行っています。DALEXによる解釈モデルの学習が終わったので、DALEXを用いてモデルの振る舞いを解釈していきましょう。DALEXによる解釈は、DALEX::explain()でexplainerオブジェクトを作るところから始まります。嬉しいことに、parsnipで学習したモデルはexplain()にそのまま与えることができます。created!あとはexplainerにingredientsで用意された様々な解釈手法を適応するだけです。変数重要度まずは変数重要度を見てみましょう。これはexplainerをingredients::feature_importance()に与えるだけで計算できます。ForestDALEXは学習アルゴリズムに依存しない手法がベースになっているので、変数重要度の計算もpermutationベースのものが使われます。その変数をシャッフルしてどのくらい予測精度が落ちるのかがその変数の重要度として定義とされています。ingredientsによって作成されたオブジェクトはplot()で簡単に可視化することができます。ggplot2がベースになっているので、ggplot2の設定やレイヤーを+で重ねていくこともできます。重要ではないとしていることがわかります。これはそもそものシミュレーションの設定に沿っています。PDP次はPDPを見てみましょう。PDPはモデルの平均的な予測結果を可視化したものです。これもexplainerをingredients::partial_dependency()に与えるだけで計算できます。関係をうまく捉えられていることが見て取れます。シミュレーション2データの作成PDPは交互作用がない場合はうまくモデルを解釈する事ができますが、交互作用がある場合は可視化がうまく機能しないことが知られています。以下のようなシミュレーションを考えてみましょう。ときは負の効果があるという特定化になっています。theme_scatter()DALEXによる解釈PDPシミュレーション1と同じくモデルを作成、学習し、explainerオブジェクトを作ります。Forest")PDPも同様に作成します。関係をうまく可視化できるかにフォーカスすることにします。"none")青い線がPDP、グレーの線が本来の関係です。動いても予測値にはまるで効果がないように見えます。念の為ですが、モデルの予測自体はうまくいっています。ここで起きている現象は、「モデルは交互作用を学習できているが、PDPはそれをうまく可視化できていない」というものです。PDPは交互作用を平均化してしまうため、プラスの交互作用とマイナスの交互作用が相殺して効果がないように見えてしまっています。Plot交互作用をうまく捉える手法の一つがICEになります。ICEは単純にPDPの平均化する前の出力を全てプロットするというものです。平均化をしていないので交互作用が相殺されることを防ぐことができます。例によって、ingredients::ceteris_paribus()にexplainerを与えるだけで計算できます。tibbleのまま渡すと警告が出るのでdata.frameにしている。"none")青い線の一本一本が各インスタンスの予測結果に対応しています。PDPとは違って、ICEではモデルが交互作用を学習していることをうまく可視化することができています。大きさによるものです。値が全てのインスタンスで共通なら一直線上に並んだきれいなX印になります。PDPPDPは交互作用のある変数についても平均化してしまうのが問題でした。ICEはひとつの解決策でしたが、これはこれでやたら線が多くなってしまいますし、値の安定性も低くなります。効果が逆になることが問題でした。値でグループに分けてPDPを計算することで、交互作用の問題を解決できそうです。これを行うための関数がingredients::aggregate_profiles()として用意されています。グループを指定するとPDPをグループごとに求めることができる。"none")単純なPDPでは捉えられなかった交互作用を可視化することに成功しました。現実のシチュエーションでは、たとえば男女でグループ分けすることで、モデルが男女の効果の違いを捉えているのかを確認することができます。グループ分けすればいいことがわかっていましたが、実際にはどの変数でグループ化すればいいかわからない場合も多いかと思います。ingredients::cluster_profiles()を用いることで、似たようなICEをクラスター化してPDPを計算してくれます。クラスター数を指定する。グループ化した場合と同じく、単純なPDPでは捉えられなかった交互作用を可視化することに成功しました。まとめ本記事では、tidymodelsを用いて機械学習モデルを作成し、それをDALEXとingredientsを用いて解釈する方法をまとめました。一つの重要なパッケージであるiBreakDownは別の記事でまとめたいと思っています。本記事で使用したコードは以下にまとめてあります。searchでやっっていました。tidymodels配下のdialsとtuneを用いることで、より簡単にハイパーパラメータのサーチを行えるので、本記事ではその使い方を紹介したいと思います。パラメータサーチ以外のtidymodelsの使い方には本記事では言及しないので、以下の記事を参考にして頂ければと思います。dropout009.hatenablog.comdropout009.hatenablog.com前処理まずは前回の記事と同様、rsampleで訓練/テストデータの分割を行います。例によってデータはdiamondsを用います。そのためにTrain/Validationに分割されたデータ、前処理レシピ、学習用モデル、ハイパーパラメータの候補の4つが必要になります。最初の3つは以前の記事で触れているので、4つ目のハイパーパラメータに関する部分に詳しく触れていきます。Train/Validationデータこれはrsamplesの仕事です。Validationで評価するので、それに合わせます。Fold5前処理レシピ次に、前処理のレシピを作成します。Forestを使うので、前回の分析同様、前処理は最低限にしておきます。all_nominal学習用モデルparsnipでモデルを設定します。Forestを使いましょう。サーチしたいパラメータはここでは値を決めずtune::tune()を与えます。今回はmin_nとmtryをサーチすることにします。ハイパーパラメータここが本記事での新しい内容になります。dialesにはparsnipで指定できるハイパーパラメータに関して、探索するレンジを指定するための関数が用意されています。Forestのハイパーパラメータmin_nに対応する関数で、デフォルトだと2-40のレンジでハイパーパラメータが探索されます。これは最終ノードに最低でも必要なインスタンスの数を表していて、これを大きくするとより強い正則化がかかります。40]自分でレンジを決めたい場合は引数で指定することができます。10]同様に、もうひとつのハイパーパラメータmtryに関しても見てみましょう。こちらは各ツリーでの分割の際に用いる特徴量の数で、これを小さくするとより強い正則化がかかります。?]min_n()とは違って、レンジの最大値が指定されていません。モデルに投入する特徴量の数よりも大きいmtryを探索しても意味がない*1ので、こちらで直接指定してあげる必要があります。5))のようにレンジを指定することもできますが、実際にモデルに投入するデータフレームを与えてあげることでレンジを指定することもできます。9]dials::finalize()の第1引数にレンジを指定したい関数を、第2引数に特徴量のデータフレームを指定することで、適切なレンジが指定されます。特徴量の数は分析の途中で変動するので、このやり方は柔軟性があって良いんじゃないかと思います。これで探索範囲の指定ができるようになりました。次に、探索したいハイパーパラメータのリストを作ってtune::parameters()に与えるとparametersオブジェクトを作ることができます。nparam[+]このparametersオブジェクトをdials::grid_*()に渡すことで、実際に探索するハイパーパラメータの値を作ることができます。grid_regular()ならグリッドサーチ、grid_random()ならランダムサーチになります。今回はランダムサーチにしましょう。サーチの数はsizeで指定できます。これまでに作ったオブジェクトをtune::tune_grid()に渡します。T))object:recipeで作成された前処理レシピを渡します。resamplesに与えたデータがこの前処理を済ませたあとでモデルに投入されます。model:parsnipで定義した学習用のモデルです。resamples:rsamplesで作った学習/評価用のデータを渡します。Varidationで評価すると思うのでrsamples::vfold_cv()で作ったデータフレームを渡すのがいいと思います。grid:dialsとtuneで作ったハイパーパラメータの候補が格納されたデータフレームを渡します。metrics:精度の評価指標です。yardsticに準備されている関数を指定することができます。control:指定しなくても構いませんが、今回はログが出力されるようにしています。学習/評価が終わると、以下のようなデータフレームが手に入ります。1]>.metricsに予測精度が格納されています。unnest()でもとってこれますが、手っ取り早い関数としてtune::collect_metrics()が準備されています。rows特に精度の高いハイパーパラメータの候補が知りたい場合は、tune::show_best()で確認することができます。0.000679一番精度の高かったハイパーパラメータを使って全訓練データでモデルを再学習する場合は、tune::select_best()でハイパーパラメータをとってきてupdate()でモデルをアップデートできます。あとは通常通り学習・評価を行うだけです。まとめ本記事ではtidymodels配下のdialsとtuneを用いたハイパーパラメータのサーチを紹介しました。Validationデータにどうアクセスするか結構ややこしかったと思うのですが、tuneを用いることでよりスッキリと探索と評価を行うことができます。dialsとtuneは公式ドキュメントに個別の使い方はまとまっているのですが、目的がわかりにくい部分や、どの情報が最新かわからないところもあり、典型的なタスクに対しての用途を自分でまとめ直したという経緯になります。optimizationを用いたパラメータ探索も可能となっています。確認してみて下さい*2本記事で使用したコードは以下にまとめてあります。Controlのいいとこ取りになっています。パネルデータから因果推論を行う際に非常に強力な武器になると思ったので、本記事でコンセプトをまとめ、論文内の比較実験を再現しました。時点でのデータがとれているパネルデータを考えます。時点のみ処置を受けるケースを考えます。処置後(post)となります。書くことにしましょう。処置を受けた場合と受けなかった場合の差をみることで、処置がアウトカムに与えるインパクトを測定することができます。観測することができません。推定する必要が出てきます。これを反事実と呼びます。Difference(SDID)を比較していきます。推定できるのではというのがDIDの発想です。それぞれのグループの平均値を表しています。上式は以下のように解釈できます。処置前コントロール群のみを用いて予測を行うと(第一項)、コントロール群と処置群の違いでバイアスがかかり、次に処置前と処置後の違いでもバイアスがかかります。そこで、コントロール群<->処置群の差を第二項で、処置前<->処置後の差を第三項で補正しています。この意味において、DID推定量はバイアスを二重に補正していると言えます。差を見ることで処置のインパクトが推定できます。(SC)DIDは単純平均でしたが、コントロール群の加重平均で処置群を表現しようというのがSCの発想になります。SCは以下の2ステップを踏みます。うまく処置群を近似できるように決めます。\omega_{i}=1\right\}\end{align}は加重平均の重みなので0以上、足して1の制約がかかっています。推定を安定させています。見比べると:DIDでは単純平均でしたが、SCは加重平均を用いることで処置群の近似を改善しています。その一方で、SCにはDIDにあった処置前後のバイアス補正(第三項)が存在しません。(SDID)DIDとSCを見比べることで、双方の利点と欠点が見えました。そこで、SDIDでは両方のいいとこ取りをします。よい推定を目指します。同様に以下の最小化問題を解くことで求めます。これは重みとしては用いません。SDIDを比較した実験を行っているので、それを再現した結果を紹介します。Controlを用いて推定しています。使用データはカリフォルニアを含むアメリカ39州の年度別喫煙量のデータで、1970年から2000年にかけて31年分あります。禁煙法は1989年から実施されています。うまく予測できるかを確かめたいので、処置前(1988年以前)のデータに関して、各州の1期先の喫煙量を他の38州の喫煙量からDID/SC/SDIDで予測し、精度を検証します。正確に推定できると言えます。対して1980年から1988年までの9年分の予測をDID/SC/SDIDで行い、RMSEを計算します。\sum_{t=1980}^{1988}\left(Y_{it}-\hat{Y}_{it}\right)^{2}}\end{align}その結果が以下になります。予測することに成功しています。中央値で見るとDIDと比較して50%、SCと比較しても15%の改善です。具体的にどの部分で精度に差が出ているのかプロットを見て確認してみます。上図は各州の喫煙量の実測値とDID/SC/SDIDの推定値がプロットされています。横軸が年度で縦軸が喫煙量です。図が細かいのですが、基本的にSDID(赤)はDID(オレンジ)/SC(緑)よりも実測値Y(青)をうまく予測できています。差が顕著な部分をいくつか抜き出したものが下の図になります。HampshireやUtahなど、他の州の加重平均で表現できない州は予測がうまくいっていません。これに対して、SDIDは非常にロバストな予測ができています。Differencesのコンセプトをまとめ、実際に精度が改善することを確認しました。SDIDはRCTが行えない状況でパネルデータの因果推論を行う際、非常に強力な武器になると思います。無視することにします。勾配ブースティングまとめ参考文献はじめに今更ですが、XGboostの論文を読んだので、2章GBDT部分のまとめ記事を書こうと思います。*1この記事を書くにあたって、できるだけ数式の解釈を書くように心がけました。数式の意味をひとつひとつ追っていくことは、実際にXGBoost(またはLightGBMやCatBoostなどのGBDT実装)を使う際にも役立つと考えています。ハイパーパラメータがどこに効いているかを理解することでチューニングを効率化したり、モデルを理解することでよりモデルに合った特徴量のエンジニアリングができるのではないかと思います。この記事に限りませんが、記述に間違いや不十分な点などあればご指摘頂ければ嬉しいです。対する予測モデルを構築することを想定しましょう。*2今回はツリーをアンサンブルした予測モデルを構築します。ツリーを加法的に組み合わせた予測モデルは以下のように定式化できます。ひとつひとつのツリーを表しています。