H2Oによる機械学習：アンサンブル決定木

　新しい h2oml コマンド群により、従来の統計モデルだけでは捉えきれないデータの知見を、 H2Oを介した機械学習で明らかにできます。勾配ブースティングマシン（GBM）とランダムフォレストというアンサンブル決定木を使用して、 分類や回帰を実行できます。
　ハイパーパラメータの調整、検証データまたは交差検証（CV）の利用、モデル性能の評価、予測の説明などを行えます。

　Stataユーザーは長年、線形回帰、ロジスティック回帰、従来型の統計モデルを利用してデータから知見を得てきました。 しかし、応答変数と予測変数の関係が複雑で非線形になる状況では、こうした古典的な手法では複雑なデータパターンを十分に捉えられないことがあります。

• 線形回帰（またはロジスティック回帰）モデルでは表現できない形で予測変数が相互作用している場合
• 慎重に変数選択をしてもモデルの精度が頭打ちになる場合
• 欠損データや多重共線性に対して頑健で、観測済みデータの範囲を超えても汎化しやすいモデルが必要な場合
• 予測力のために予測の説明可能性を犠牲にする必要がない場合

　GBMとランダムフォレストは、Stataでのデータ分析方法を大きく変えます。
　Stata内から H2O の機械学習アルゴリズムにシームレスにアクセスできるため、使い慣れたStata環境を離れることなく、高性能な予測モデルを利用できます。 直感的なStata構文のコマンドを使用するだけで、従来手法を上回る高度なアンサンブル学習モデルを学習できます。
　機械学習は「ブラックボックス」でしょうか。もうそうではありません。 SHAP（Shapley additive explanations）値、部分依存プロット（PDP）、変数重要度ランキングなどのツールにより、 GBMとランダムフォレストは強力な予測を提供しながら説明可能性も維持します。

StataでH2Oを使用するための準備

例題データセット：Telco顧客解約データセット

　h2oml コマンド群のデモを行うために、バイナリ分類モデルの GBM に焦点を当てます。 ランダムフォレストや回帰用の GBM などの他のモデルは、構文とワークフローが似ているため、同様の手順で実行できます。
　カリフォルニアで自宅電話とインターネットサービスを提供する Telco という架空の通信会社のデータを分析します。 使用するデータセットには、26 の変数を持つ 7,043 人の顧客に関する情報が含まれています。
　分析の目的は、どの顧客に解約のリスクがあるか (つまり、Telco とのサービスを終了する可能性があるか) を識別できる予測モデルを開発することです。 バイナリ応答である解約は、顧客が過去 1 か月間にサービスを中止したか否かを示します。
　予測子 (特徴) には、顧客の人口統計 (性別や年齢など)、アカウントの詳細 (契約タイプや契約期間など)、 サービス サブスクリプション (インターネット、電話、セキュリティ サービスなど) が含まれます。
　潜在的な解約の兆候となる顧客行動パターンを予測し、通信事業者が貴重な顧客を維持するための積極的な措置を講じるのに役立ちます。

　Stata で H2O 機械学習用にデータを準備します。まず、データセットを読み込み、H2O クラスターを開始します。

　次に、現在のデータセットを H2O フレーム churn に入れて、現在の H2O フレームに設定します。

　フレームを 80% と 20% の比率で分割し、トレーニング フレームとテスト フレームとして使用します。 後ほど、トレーニングフレームに CV （交差検証法）を使用して、オーバーフィッティング（過学習）を制御します。

　便宜上、予測子の名前を格納するグローバルマクロを作成します。

ベースとなる GBM モデル

　

　再現性を担保するため h2orseed(19) を使用し、3 層別 CV (SCV) で GBM モデルを適合させます。その他の設定はデフォルト値を使用します。 SCV は、不均衡なクラス (解約レベル) の分類タスクに特に役立ちます。

　出力結果の上部にはモデルの詳細が示されており、h2oml gbbinclass がベルヌーイ損失を使用していることがわかります ( h2oml gbregress による回帰には他の損失関数も使用できます)。
　トレーニングフレームには 5,643 の観測値が含まれており、トレーニングデータに対して 3 層 SCV が使用されています。
ハイパーパラメータセクション (モデルパラメータ) には、ハイパーパラメータのユーザー指定値と実際の値の両方が報告されます。 早期停止が指定されている場合は、これらの値が異なる場合があります。

　Metric summary の表には、トレーニングデータと CV データのバイナリ分類パフォーマンス指標が表示されます。 クラスが不均衡であるため、指標として、適合率 - 再現率曲線の下の領域 (AUCPR) を使用します。 AUCPR の範囲は 0 ～ 1 で、1 は完璧なパフォーマンスを意味します。 CV の指標が主な焦点ですが、トレーニング指標もチェックして、多少の過学習は許容しつつ、学習不足（アンダーフィッティング）になっていないことを確認します。 トレーニングと CV における AUCPR の差は正 (0.8024 - 0.6585 = 0.1439) であることが予想されますが、 大きなギャップは過学習を示している可能性があり、モデルの汎化性能が失われている可能性があります。

　後で比較できるように、推定結果を保存します。

　3 分割交差検証（3-fold CV）において、各foldの評価指標のばらつきを確認することで、モデルの性能が特定のデータ分割に依存していないかを検証できます。 CV 指標の変動が大きい場合、汎化性能が不十分である可能性があります。 これは、h2omlestat cvsummary コマンドを使用して確認できます。

モデルの選択とハイパーパラメータの調整

　

　ベースラインモデルの CV AUCPR は 0.6585 でした。モデルの性能を向上させるために、ハイパーパラメータを調整します。 GBM には 11 個の調整可能なパラメータがありますが、簡潔にするために、 ntrees() (ツリーの数) と predsamprate() (予測子のサンプリングレート) のみを調整します。 ハイパーパラメータの調整は反復的な手順であり、この例は説明のみを目的としています。

　出力結果の上部には、チューニング方法 (Cartesian)、チューニング指標 (AUCPR)、ハイパーパラメータのグリッド検索範囲などのチューニングの詳細が表示されます。 選択された値 (ntrees() は 100、predsamprate() は 0.15) は、最も性能が良いモデルに対応します。 出力の残りの部分には、最適なモデルのハイパーパラメータ値と指標の概要が表示されます。
　 　チューニングにより、CV AUCPR が 0.6585 から 0.6739 に増加しました。この例ではハイパーパラメータ空間のごく一部しか探索していないため、改善はわずかです。

　この調整されたモデルを後で使用するために保存しましょう。

　h2omlestat gridsummary コマンドを使用して、グリッド検索の概要を取得できます。 このコマンドは、調整中のハイパーパラメータの構成を AUCPR でランク付けして一覧表示します。

　h2omlexplore コマンドを使用して、他の指標に基づいて最初の 2 つのモデルを比較することもできます。

　最適に調整されたモデル以外のモデル (性能がわずかに劣るが、使用するツリーが少ないモデルなど) を選択したい場合は、 h2omlselect コマンドを使用してそのモデルを選択できます。

　h2omlgof コマンドを使用して、他の指標に基づいて、最適なモデル gbm_tuned と 前のセクションのベースモデル gbm_default を比較してみましょう。

　出力結果には、トレーニングの結果と、その下に CV の結果が表示されます。 CV の結果を見ると、チューニングによって全ての指標において性能が向上したことがわかります。 チューニングしたモデルのログ損失、クラスごとのエラー、MSE、RMSE が低く、 AUC、AUCPR、ジニ係数が高いことは、 モデルの性能が向上したことを示しています。

　h2omlgraph varimp コマンドで変数の重要度順に予測子の一覧を表示できます。 初期設定で上位10変数を表示しますが、表示する変数の数を変更したい場合は top() オプションで指定します。

　変数重要度プロットによれば、tenuremonths と contract が予測に最も大きく寄与しています。 表示された予測子（上位10変数）の中では onlinebackup の重要度が最も低いことがわかります。

手法の選択 : GBM とランダムフォレスト

　

　前と同じ手順に従って、h2oml rfbinclass を使用してバイナリ分類用のランダムフォレストを実行します。 簡単にするために、ハイパーパラメータのチューニング後の最適なモデルを保存します。

　次に、テストフレーム (データ分割時に作成したテストデータセット) を使用して、調整されたランダムフォレスト (rf_tuned) と 調整された GBM メソッド (gbm_tuned) を比較します。 h2omlpostestframe コマンドを使用して、フレームの名前 (この場合は test) を指定します。 このフレームは、rf_tuned と gbm_tuned の両方の結果を計算する後続のすべての事後推定コマンドで使用されます。

　h2omlestat metrics を使って各手法のテストデータ（test frame）における性能指標を個別に表示する代わりに、 ここでは h2omlgof を使用して、複数のモデルの結果を並べて比較します。

　

　GBMはランダムフォレストよりもAUCPRが高いため、より優れた手法と判断されます。 さらに、ROC曲線を用いて両手法を比較することもできます。 ROC曲線では、曲線が左上隅に近いほど性能が高いことを示します。

　ROC の結果に基づくと、GBM はランダムフォレストよりわずかにパフォーマンスが優れています。
　 異なる方法やモデル間で分類予測を比較するもう 1 つの一般的なアプローチは、正しく予測された結果と誤って予測された結果の数を報告する混同行列を使用することです。 h2omlestat confmatrix コマンドにより、混同行列を表示することもできます。

新しいデータでの予測

　

　新しいデータ newchurn.dta が得られたとします。 GBM モデル gbm_tuned に基づいて、これらの新しい顧客の解約確率を予測したいと考えています。 新しいデータセットを H2O フレーム newchurn として読み込みます。

　h2omlpredict を使用して、解約確率とクラスを予測します。 デフォルトでは、クラス (「はい」または「いいえ」) を予測します。 確率を取得するには、pr オプションを指定します。
　前のセクションでは、h2omlposttestframe を使用して、test を事後推定フレームに設定しました。 新しいデータセットの予測を取得するには、frame(newchurn) を指定します。 以下では、GBM モデル gbm_tuned を使用して、クラスと確率の両方を予測します。

　frame(newchurn) を指定したため、クラス (churnhat) とクラス確率 (churnprob1 と churnprob2) 用に生成された変数は、 newchurn フレームで使用できます。予測クラスと確率の値をいくつか挙げてみましょう。

　

　たとえば、churnprob2 は最初の顧客の解約確率が 22%、2 番目の顧客の解約確率が 78% であることを示しています。 churnhat の予測クラス (「はい」または「いいえ」) は、churnprob2 が デフォルトの F1 スコアを最大化する閾値（ 0.2378 ）を超えているかどうかに基づいて割り当てられます。 この閾値は、 h2omlestat threshmetric コマンドを使用して確認できます。 このコマンドは、各評価指標が最も高くなる判定閾値を一覧表示します。 閾値を指定するには、h2omlpredict の threshold() オプションを使用します。

予測結果の説明

　

　機械学習では、モデルがなぜその予測結果を出したのかを理解することが重要です。 説明可能性（Explainability）を確保することで、予測結果の精度だけでなく、その判断根拠を人が理解し説明できるようになります。 グローバルな説明手法は、個々の予測ではなく、モデル全体がどのような特徴を重視し、どのような傾向で予測を行っているかを示します。 代表的な例として、以下のような手法があります。

• 変数の重要度
• グローバルサロゲートモデル (機械学習の予測を近似する、シンプルで解釈可能なモデル)
• PDP

　ローカルモデルは、単一の観測に対するモデルの動作を近似することで、個々の予測を説明します。 例としては、以下のような手法があります。

• ICE 曲線
• SHAP 値

グローバルな説明可能性の手法

　変数重要度グラフの例はすでに確認したので、ここではモデル全体の挙動を理解するためのサロゲートモデルの作成方法を説明します。 まず、前の手順で学習した最良のGBMモデル（gbm_tuned）を読み込みます。
　サロゲートモデルとは、複雑で解釈しにくいモデル（GBM、ランダムフォレスト、ニューラルネットワークなど）の振る舞いを、 より単純で解釈しやすいモデル（決定木など）で近似して説明するためのモデルです。

　次に、h2omlpredict コマンドの frame(churn) オプションを指定して、 GBM モデルを使用して、テストフレームだけでなくデータセット全体の解約を予測します。

　説明可能性を向上させるために、決定木を使用して、モデル gbm_tuned から churnhat の予測を近似するグローバルサロゲートモデルを構築します。

　まず、作業フレームを churn フレームに切り替えます。

　次に、説明のために、最大深度 3 (maxdepth(3)) の単一の決定木 (ntrees(1)) をグローバルサロゲートモデルとして学習します。

　この決定木は、複雑な GBM モデルの簡略化された近似として機能し、予測をより解釈しやすくしながら、有用な洞察が得られます。

　決定木は視覚的に解釈しやすい手法です。 h2omltree コマンドの dotsaving() オプションを使用すると、 Graphviz を使用してプロットできる DOT ファイルを生成して、より優れた視覚化を実現できます。

　終端ノードの値は、顧客が解約しない確率 (No) を表します。 たとえば、インターネットサービスがないか、オンラインバックアップとセキュリティを使用している 1 年または 2 年の契約を持つ顧客は、 留まる確率 (0.997) が最も高くなります。つまり、解約確率は 0.003 (1 - 0.997) に過ぎず、離脱する可能性が最も低くなります。

　次に、重要な予測子 (以前に h2omlgraph varimp で特定) が churn にどのような影響を与えるかを分析します。 選択した予測子の予測に対する限界効果を示すグローバルな説明可能性手法である PDP を使用します。 最初に GBM モデルの推定結果を復元して、後続の事後推定コマンドが最適な GBM モデル (gbm_tuned) に適用されるようにします。

　h2omlpostestframe を notest オプションとともに使用して、churn フレームを テストフレームとして扱わずに、事後推定分析のアクティブフレームとして設定します。

　次に、主要な予測子の PDP (Partial Dependence Plot; 部分依存プロット) を生成します。
　PDPは他の変数を平均的な値に固定したとき、その変数だけを変化させると予測値（ここでは解約確率）がどう変わるかを示します。

　

　縦軸の Partial dependence は、ここではモデルが予測した解約確率となり、赤い折れ線に対応します。
　PDP パターン (赤い線) は、サロゲートツリーの結果と一致します。 たとえば、解約の確率 (y 軸に表示) は、1 年または 2 年の契約 ( contract ) を持つ顧客と、 会社のサービスを長期間 ( tenuremonths ) 使用する顧客で低下します。
　青い棒（ヒストグラム）はデータの分布を示しています。 棒が高いほどその領域の観測データ数が多く、PDPの結果を信頼しやすいことを意味します。 逆に棒が低い領域はデータ数が少ないため、赤線の変動は慎重に解釈する必要があります。

ローカルな説明可能性の手法

　

　ローカルな説明可能性では、各予測因子の個々の予測への貢献を推定する SHAP 値を使用します。
　SHAP（SHapley Additive exPlanations）値 は、機械学習モデルの予測結果に対して、各変数がどれだけ予測値を押し上げたり押し下げたりしたかを定量的に示す指標です。 SHAP 値を確認することで、特定の顧客が解約すると予測される理由を説明するのに役立ち、機械学習の決定を透明化します。

　ここでは、h2omlgraph shapvalues を使用して、観測行 19 に注目し、 上位 10 の SHAP 重要予測因子の SHAP 値を生成します。 観測行 19 は、月単位の契約サービスを 9 か月間使用し、観測された解約と予測された解約の両方の値が Yes である女性顧客です。

　青いバーは解約確率を高める予測因子を示し、赤いバーは解約確率を下げる予測因子を示します。
　グラフに示された [E(f(x))] = -1.593 は、全顧客の平均的な予測値を示し、この値がベースラインとなります。
　SHAP 値は、前述の分析結果と一致しています。特に、月単位の契約（ contract ）、短い契約期間（ tenuremonths ）、 オンラインセキュリティサービス（ internetserv ）を使用しないことが、 この特定の顧客の解約にプラスの影響を与えています。
扶養家族（ dependents ）がいるため離反しにくい傾向はありますが、月契約や光回線利用などの要因がそれを上回り、 最終的には離反（Yes）と予測されたことが読み取れます。
　バイナリ分類の SHAP 値はロジット スケールで報告されることに注意してください。 たとえば、生の予測値 f(x) = 0.2063 は、「観測値 19 が解約する」という予測確率のロジットです。 これを確率として解釈するには、逆ロジット変換する必要があり、約55.1%の確率となります。

　SHAP サマリー (beeswarm) プロットは、すべての観測値にわたる SHAP 値を視覚化し、予測因子の重要性と応答への影響の両方を示します。 説明のために、最も重要な上位 4 つの予測因子をプロットします。

　y 軸に表示される予測子は、SHAP の重要度によってランク付けされ、top(4) オプションにより上位4つが表示されます。 契約期間（contract）が月単位の顧客（青色）は正のSHAP値を示しており、解約予測を押し上げる方向に働いています。 一方、1年契約や2年契約（紫～赤色）は負のSHAP値を示し、解約予測を下げる方向に働いています。 よって、契約期間が短いほど顧客が解約する可能性が高くなることが伺えます。

H2O クラスターのシャットダウン

　

　分析が完了したら、次のコマンドを使用して H2O クラスターから切断できます。

　これにより、Stata-H2O セッションは終了しますが、クラスターはバックグラウンドで実行されたままになります。後で次のコマンドを使用して再接続できます。

　

　クラスターを完全にシャットダウンしてすべてのリソースを削除するには、次のコマンドを使用します。

　

参考

　さらに詳しい内容につきましては、下記のマニュアルをご覧ください。

開発元マニュアル