causallibを使ってIPWを計算してみる
概要
IBMが公開している「causallib」というパッケージを使ってIPW(Inverse Probability Weighting)を推定してみました。
こちらのサンプルを実際にやってみます。
目次
1.最初に
今回は、「IPW(Inverse Probability Weighting)」をPythonで求めたいと思い、IBMが公開している「causallib」というパッケージを使ってみました。
詳細の解説は省略しますが、超ざっくり言うと、IPWは傾向スコアの値を利用してATE（Average Treatment Effect）を推定する方法です。(介入グループに所属するか否かを推論して、その値を重みとして利用してATEを推論する)
この辺りを勉強したい方は効果検証入門や岩波データサイエンスがオススメです。
特に、「効果検証入門」はとてもわかりやすく解説をしてくれているので、とっかかりとしてはこちらから読み始めるといいのかと思います。
なお、今回利用しているライブラリのバージョンは、下記の通りです。
import causallib import sklearn print(causallib.__version__) print(sklearn.__version__)
0.6.0 0.22.1
2.やってみた
今回利用するデータセットは「喫煙と体重の増減」についてのデータセットでした。
This loads an object in which data.X, data.a, and data.y respectively hold the features for each individual, whether they stopped-smoking, and their observed difference in weight between 1971 and 1983.
参照:Datasets
個人ごとに「共変量」、「介入の有無(喫煙を辞めたか否か)」、「目的変数（体重の変化）」がそれぞれ用意されているのですが、このデータを使って「喫煙を辞めることによる体重への影響」の効果検証をする、というものです。
以下、実際にサンプルを参考に一部を抜粋して記述していきます。
（自分の理解を深めるために、一部サンプルとは異なる記述をしています）
2-1.データの確認
まずは、今回使うデータがどんなデータなのか、をざっくり確認してみます。
データはPandasのDF、Seriesで保持しています。
共変量
print('●サイズ', '\n',data.X.shape)
print('●型', '\n',type(data.X))
print('●データを一部確認')
data.X.head()
介入
print('●サイズ', '\n',data.a.shape)
print('●型', '\n',type(data.a))
print('●count', '\n',data.a.value_counts())
●サイズ (1566,) ●型 <class 'pandas.core.series.Series'> ●count 0 1163 1 403 Name: qsmk, dtype: int64
目的変数
print('●サイズ', '\n',data.y.shape)
print('●型', '\n',type(data.y))
print('●統計情報', '\n',data.y.describe())
print('●データを一部確認')
print(data.y.head())
●サイズ (1566,) ●型 <class 'pandas.core.series.Series'> ●統計情報 count 1566.000000 mean 2.638300 std 7.879913 min -41.280470 25% -1.478399 50% 2.603811 75% 6.689581 max 48.538386 Name: wt82_71, dtype: float64 ●データを一部確認 0 -10.093960 1 2.604970 2 9.414486 3 4.990117 4 4.989251 Name: wt82_71, dtype: float64
2-2.とりあえず一連の流れを実行する
まずは最小要件として、一連の流れを実行してみます。
from causallib.datasets import load_nhefs
from causallib.estimation import IPW
from causallib.evaluation import PropensityEvaluator
from sklearn.linear_model import LogisticRegression
# 1.データの取得
data = load_nhefs()
# 2.スコアを推論するモデルの指定と初期化
learner = LogisticRegression(solver="liblinear")
ipw = IPW(learner)
# 3.学習
ipw.fit(data.X, data.a)
# 4.介入有無ごとに重み付けをした上で結果を計算
outcomes = ipw.estimate_population_outcome(X = data.X,
a = data.a,
y = data.y
)
# 5.介入有無間で差分を取得
effect = ipw.estimate_effect(outcome_1=outcomes[1], # 喫煙をやめた
outcome_2=outcomes[0], # 喫煙を続けた
effect_types="diff" # 差分を取得する
)
print(effect)
diff 3.47451 dtype: float64
喫煙をやめたことで、11年で体重が約3.5Lbs（ポンド）増加した、という解釈ができます。
(今回のデータセット全体に対して、「施策を実施した場合-しなかった場合」を求めたイメージです。ATEなので。)
当然、本来ならもっと細かく検証等をしないといけないのですが、とりあえず最低限の一連の流れが掴めたと思います。
2-3.パラメータを少しいじる
「2-2.とりあえず一連の流れを実行する」では、諸々デフォルトのパラメータを使っていましたが、もう少しパラメータを弄ってみます。
ここで注意したい点としては、IPWはその計算式の関係上、推論したスコアが「0に近すぎる」等の場合に、そのデータの影響度が大きくなりすぎる問題があることです。
この問題に対応するため、「推論したスコアの値が小さすぎる、もしくは大きすぎる場合に指定した閾値に設定」しています。
今回の例だと、推論した値が「0.3未満、もしくは0.7よりも大きい」場合はそれぞれ「0.3,0.7」に調整しています。（実際に利用する際は、こんな値に設定することは殆ど無いかと思いますが...、結果の確認をしやすくするためにこのような値を設定しています）
# 1.スコアを求めるモデルのパラメータ
learner = LogisticRegression(penalty="l1",
C=0.01,
max_iter=500,
solver='liblinear'
)
# 2.上記のモデルで求めたスコアが小さすぎたり大きすぎたりした場合に、指定した範囲の値にする
# truncate_eps ~ (1-truncate_eps)の範囲におさめる
truncate_eps = 0.3
ipw = IPW(learner=learner,
truncate_eps=truncate_eps
)
# 3.モデルの学習
# Xを元にaを予測するモデル
ipw.fit(X=data.X,
a=data.a
)
# 4.「2」で指定した値に基づいて、スコアが調整されていることを確認する
probs = ipw.compute_propensity(X=data.X,
a=data.a,
# 喫煙を辞めた確率を取得する。
# デフォルト（None）だと、aの最大値（今回の場合1）が利用される。これは、一般的に介入は「0:非介入」、「1:介入」としてデータを入れることが多いため
# 今回の場合は指定しなくても結果は同じものになります
treatment_values=1
)
probs.describe()
ちゃんと全てのスコアが「0.3~0.7の範囲内」におさまっていることが確認できます。
Fraction of values being truncated: 0.74266. count 1566.000000 mean 0.317264 std 0.042961 min 0.300000 25% 0.300000 50% 0.300000 75% 0.301921 max 0.700000 Name: 1, dtype: float64
また、「2-2.とりあえず一連の流れを実行する」では介入があるグループ、無いグループ間での単純な差分を結果として取得していましたが、他にも「比」や「オッズ比」も取得できます。
目的変数やタスクの内容によって、「どの結果を取得するか」を選択することができますね。
outcomes = ipw.estimate_population_outcome(X=data.X,
a=data.a,
y=data.y
)
effects = ipw.estimate_effect(outcomes[1], # 喫煙をやめた方
outcomes[0], # 喫煙を続けた方
effect_types=["diff", # 差分
"ratio", # 比
"or", # オッズ比
]
)
print(effects)
Fraction of values being truncated: 0.74266. diff 2.943005 ratio 2.584348 or 0.583127 dtype: float64
2-4.評価をする
分析の品質を担保するために、もう少し評価をしっかりやる必要があります。
今回参照したサンプルでは、下記の６つ*(学習/評価用データ)=12個の描画結果を確認していました。
from sklearn import metrics
plots=["roc_curve", # ROC曲線。
"pr_curve", # PR曲線。
"weight_distribution", # スコアの分布を介入有無グループごとに確認する
"calibration", # 「calibration-curve」
"covariate_balance_love", # 重み調整前後での共変量のバランス（介入有無間）を確認する
"covariate_balance_slope" # 重み調整前後での共変量のバランス（介入有無間）を確認する
]
metrics = {"roc_auc": metrics.roc_auc_score, # ROC曲線のAUC
"avg_precision": metrics.average_precision_score, # 適合率
}
ipw = IPW(LogisticRegression(solver="liblinear"))
evaluator = PropensityEvaluator(ipw)
results = evaluator.evaluate_cv(X=data.X,
a=data.a,
y=data.y,
plots=plots,
metrics_to_evaluate=metrics
)
上記のコードで出力された描画結果について、ざっくりですが下記に記述します。
また、今回は明示的に指定していないですが、「causallib=0.6.0」バージョンの時点ではデフォルトで「stratified foldsで５分割して交差検証した結果」を元に描画してくれています。
（上記の画像では「学習用データ」の結果のみが出力されていますが、同様に「評価用データ」での描画もなされます。）
- 上段左
- ROC曲線
- 予測確率の品質評価に用いる。
- 上段右
- PR曲線
- 予測確率の品質評価に用いる。
- 中段左
- 介入グループごとのスコアの分布。
- 予測確率の品質評価に用いる。
- 同じような分布をしていることが好ましい。
- 中段右
- CVごとのcaclibration-curveを出力。
- 予測確率の品質評価に用いる。
- 「x=y」と書いてあるグレーのラインに近いことが好ましい
- 厳密にやるならもう少し精査が必要な結果かもしれないが、一旦サンプルということで先に進める。
- 下段左
- 「介入グループごとの共変量の平均値」の差分を標準誤差で割ったもの、の絶対値。
- 介入グループ間における共変量のバランスが取れていることが好ましいため、値は0に近い方がいい。
- スコアに基づく重みの調整前後を表示しており、調整した後は0に近づいているのがわかる。
- 下段右
- 「下段左」と同じものを違う可視化方法で示したもの。
また、上記以外にもfoldごとの評価指標の値を取得することができます。
results.scores.prediction_scores
今回はとりあえず一連の流れで使ってみることを目的としているのでここまでとしますが、実際に利用する際は上記の評価項目に沿って適宜修正が必要です。
3.まとめ
分野的にPythonよりはRを利用する人が多い分野かもしれませんが、自分はPythonが好きなのでこのようなパッケージは大歓迎です。
「効果検証」というものはどんな事業をやるにしても大事なことだと思うので、引き続き頑張って勉強を続けていきたいと思います。
