JMDC データサイエンティストの齋藤です。
データ分析の第一歩、EDA（探索的データ分析）にどう取り組んでいますか？ 予測のための機械学習の話はよく聞きますが、EDAのための機械学習はあまり目にしないと感じるので、 今回は実務における「XGBoost+SHAPによるEDA」の実践例を取り上げてみたいと思います。 題材は２０２１年７月にリリースした「新型コロナウイルス感染時の重症化リスクファクターに関する分析結果」です。 https://www.jmdc.co.jp/wp-content/uploads/2021/07/news20210709_2.pdf このブログの内容はテクニカル中心ですが、分析結果自体も面白いのでレポートもご覧いただけると嬉しいです。
XGBoost+SHAPでEDAする理由
EDAはデータ分析の第一歩ですが意外と面倒な事が多いです。
XGBoost+SHAPの枠組みを活用することで非線形性・交互作用・データ欠損に対応したうえで、
自動的かつ網羅的に特徴量とターゲットの関係性を把握することができます。
機械学習で重要な特徴量を絞り込んだうえで、その特徴量について掘り下げていくことで、効率的かつ的確なEDAが可能になります。
なお今回は医学系論文での使用例も多いXGBoostを採用していますが、もちろんLightGBMでも同様の分析が可能です。
分析デザインの概要
新型コロナ入院患者を母集団、重症化したか否かをターゲットに、リスクファクターと考えられる特徴量を調査しました。
データベース
弊社DPC調査データに基づく約1,400万人分の医療機関データ
分析対象者
2021年2月末までに新型コロナを理由とした入院患者 7,373 名
ターゲット
入院患者のうちICU（集中治療室）入室者185名（正例2.5%の２値分類）
特徴量
下記の19特徴量
Feature Importance
特徴量とターゲットの関係性の把握といえば、まずはFeature Importanceです。
算出方法は大きく３種類ありますが（gain, permutation, SHAP）、それぞれ適した場面が違うと考えています。
gain
- モデルの学習時その特徴量により損失関数をどれだけ減らしたか。
- XGBoostアルゴリズムに沿った算出方法。ロジスティック回帰などのGLM系では同じ枠組は使えない。
- 予測精度改善にはこの重要度が高い特徴量を中心にFeature Engineeringすると良いか。
permutation
- その特徴量がランダムシャッフルされたらどれだけ精度が落ちるか。
- モデル種類によらず計算可能（tree系でもGLMでも）であり比較可能性に優れている。
- 特徴量どうしの組み合わせを考慮しないので計算量も少ない。
- 算出結果の定性的な解釈が難しいのが課題か。
SHAP
- 各レコードで予測値と平均値の差を各特徴量に分解し、それを全レコードで積みあげたもの。
- 構築されたモデルの振る舞い（特徴量と予測値の関係性）の把握という点では最もストレート。
- 計算量が多いという課題があるが、GBDT系モデルにおいてはtree SHAPと呼ばれる高速アルゴリズムが考案されている。
今回はXGBoostで特徴量と予測値の関係性を把握したいのでSHAPによるFeature Importanceを採用しています。
実際にFeature Importanceを確認してみるとBMIが最も高く、次点がDRUG_DM_FLG（糖尿病治療薬の服薬有無）となっています。
しかしFeature Importanceだけでは特徴量がどのような値をとったときに重症化リスクが上昇するのか、またその程度がどれくらいかは分かりません。
そこで特徴量と予測値の関係をより具体的に把握できるSHAP plotを確認してみます。
SHAP
SHAPは上述のとおり予測値と平均値の差を各特徴量に分解したものです。
各レコードについて各特徴量のSHAP値の和(+bias)をとれば予測確率のlogit値（log[p/(1-p)]）になります。
SHAP dependence plot
各特徴量について自身の値を横軸に、SHAP値を縦軸にplotしたものがSHAP dependence plotで、構築されたモデルにおける特徴量と予測値の関係性を一目で把握できます。
Feature Importance上位のBMIとDRUG_DM_FLGについてdependence plotを確認してみます。
- BMIは25前後でSHAP値が大きく上昇しています。BMI=25は人間ドッグ学会の定める肥満の判定ラインです。
- BMI=30あたりでSHAP値が少し低下していますが、それでも標準体(BMI:18.5～24.9)に比べればSHAP値は大きく、肥満はコロナ重症化のリスク因子になっていそうです。
- 一方で低体重(BMI<18.5)でもSHAP値の上昇が見られます。低体重による免疫力低下の影響等が仮説として考えられそうです。
- BMIのdependence plotでは右端の欠損レコードでSHAP値が非常に高くなっていることにも注目です。入院時に呼吸困難等の重篤な状態に陥っている場合に測定が省略されている等の仮説が考えられます。
- DRUG_DM_FLGは「0：服薬なし、１：服薬あり」の２値カテゴリ変数なので、糖尿病治療薬の服薬がコロナ重症化のリスク因子になっていそうです。
SHAP summary plot
SHAP summary plotは複数の特徴量のSHAP分布どうしを比較しやすくしたplotです。
縦軸が特徴量、横軸がSHAP値で、各レコードは特徴量自体の値の大小で色づけしてplotされています。
なお、このグラフで各特徴量についてSHAP値の絶対値の平均をとったものがSHAP Feature Importanceになります。
DRUG_DM_FLGは２値カテゴリ変数ということもありplotが綺麗に２分されています。 一方でSHAP値のレンジはそこまで広くなく、たとえば特徴量重要度４位のADM_DATE（入院開始日）の方がレンジは広いです。 summary plotを確認することでFeature Importanceでは見えない「特定の場面で予測値に強く影響する特徴量」も把握することができます。
XGBoost+SHAPによるEDAの注意点
XGBoost+SHAPによるEDAは便利な手法ですが注意点もあります。
- あくまで構築済モデルの特徴量と予測値の関係を表現する手法であり、学習データにおける特徴量とターゲットの関係性を直接表現したものではないこと。
- モデルが過学習していたらFeature ImpotanceやSHAPもおかしくなる可能性があるため、ハイパーパラメーターのチューニングやCV精度検証は必要であること。
- モデルの精度が低い場合は、Feature Impotance上位の特徴量でもターゲットとの関係性がそこまで強くないこと。ただし、予測精度が低くても上位特徴量を２カテゴリ区分して検定すれば有意差が出る場合もある。
- データ欠損有無とターゲットに強い相関がある場合、「特徴量の大小がターゲットに与える影響」ではない要因でFeature Impotanceが大きくなる可能性があること（今回のBMI）。
- 特徴量どうしに強い相関がある場合はSHAP値を分け合ってしまうこと。予測が目的でないなら注目したい特徴量と相関が強い特徴量は除いた方が解釈しやすい。
統計的検証への橋渡し
SHAP dependence plotだけでもEDAのゴールになり得ますが、グラフではなく定量的な表現や、学習データにおける特徴量とターゲットの関係性の直接的な検証が必要な場面もあるかと思います。
そのための統計的検証の代表は群間比較ですが、今回の分析では群を構築する部分でtreeモデルの情報を活かすアプローチをとりました。
XGBoostモデルの情報を活用した群の閾値の設定
群間比較では群を区分するための閾値をどう設定するかが重要です。
そのために今回は構築されたtreeモデルの閾値を直接調べてみました。
分岐に使用された閾値の出現回数を参考にすることでXGBoostモデルの木構造に沿った群を構築することができます。
BMIの場合には、閾値の出現回数の分布に加えて、SHAP dependence plotで「やせ」と「肥満」の双方にリスクがあると示唆されていることを踏まえ、
閾値分布の結果と整合する人間ドック学会の「やせ」と「肥満」の基準で群の閾値を設定しました。
さらに年齢調整も加味して多変量解析を行った結果、「肥満」は重症化リスク（オッズ比）1.8倍の因子となりました（p-value=0.013）。
まとめ
XGBoost+SHAPの枠組みをEDAに活用することで、手軽かつ網羅的に特徴量とターゲットの関係性を調べることができました。 GBDTモデルの振る舞いは複雑でSHAPも部分的な表現をしたものに過ぎませんが、多次元データからインサイトを得るための大きな武器になるかと思います。
