causallibでIPWをやってみる
因果推論の勉強をはじめました。
IPW法を学んだのでPythonで実際にやってみます。
causallibというライブラリを使ってIPW法をやってみました。
こちらのサンプルに沿って進めていきます。
実行環境
- Googl Colaboratory
- causallib (0.6.0)
- sklearn (0.22.2)
使用するデータの確認
喫煙をやめたことによる体重の変化を検証します。
import causallib import sklearn data = load_nhefs() data.X.join(data.a).join(data.y).head()
データの内容はざっくりと以下のようになっています。
- data.X : 共変量 (年齢、性別、喫煙年数、運動習慣など)
- data.a : 割当（喫煙をやめていれば1、やめていなければ0）
- data.y : 目的変数 (体重の増減(ポンド))
割当を確認します。
#割当
data.a.value_counts()
0 1163 1 403 Name: qsmk, dtype: int64
1163人は喫煙のまま、403人が禁煙したようです。
目的変数を確認します。
#目的変数
data.y.describe()
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
体重は平均で2.6ポンド増加したようです。
禁煙した人としなかった人での体重の変化を比較してみます。
df = data.X.join(data.a).join(data.y) print("割当1 ： ", df[df.qsmk == 1]["wt82_71"].mean()) print("割当０ ： ", df[df.qsmk == 0]["wt82_71"].mean())
割当1 ： 4.525078989702233 割当０ ： 1.9844975347463472
禁煙した人の方が体重が増加していることがわかります。
IPW法
傾向スコアの予測
サンプルの通りに実装してみます。
#train learner = LogisticRegression(solver="liblinear") #傾向スコア予測に使うモデル ipw = IPW(learner) ipw.fit(data.X, data.a) #共変量から割当を予測 #傾向スコア probs = ipw.compute_propensity(data.X, data.a, treatment_values=1)#割当１になる確率 probs.head()
0 0.131687 1 0.168435 2 0.169179 3 0.481203 4 0.251937 Name: 1, dtype: float64
ここの値は割当１になる確率です。
重み付けの確認
傾向スコアの逆数で重み付けされていることを確認します。
ipw.compute_weights(data.X, data.a).head()
0 1.151658 1 1.202552 2 1.203629 3 1.927535 4 1.336785 dtype: float64
- 割当１のデータは 1/傾向スコア
- 割当０のデータは 1/(1 - 傾向スコア)
の値になっています。
結果を確認
#割当ごとに重みづけして結果を計算 outcomes = ipw.estimate_population_outcome(data.X, data.a, data.y) #ATE effect = ipw.estimate_effect(outcomes[1], outcomes[0]) effect
diff 3.490373 dtype: float64
禁煙によって3.4ポンド体重が増加したようです。
パラメータを変化させる
次にパラメータを変化させて実行してみます。
傾向スコアが0に近すぎると、逆数をとった時に極端に大きい値になります。
また、1に近すぎると、影響が小さくなってしまいます。
今回は傾向スコアを0.2から0.8の範囲に限定します。
learner = LogisticRegression(penalty='l1', C = 0.01, max_iter=500, solver = 'liblinear') #傾向スコアを0.2から0.8に限定 truncate_eps = 0.2 ipw = IPW(learner, truncate_eps = truncate_eps, use_stabilized=False) ipw.fit(data.X, data.a) #共変量から割当を予測
傾向スコアの確認
#傾向スコア probs = ipw.compute_propensity(data.X, data.a, treatment_values=1)#割当1になる確率 print(probs.min()) print(probs.max())
Fraction of values being truncated: 0.23499. 0.7047037999493072 0.2
ちゃんと0.2から0.8の間に収まっています。
結果の確認
#割当ごとに結果を重み付け outcomes = ipw.estimate_population_outcome(data.X, data.a, data.y) #ATE effects = ipw.estimate_effect(outcomes[1], outcomes[0], effect_types=["diff", "ratio"]) effects
Fraction of values being truncated: 0.23499. diff 3.376068 ratio 2.917687 dtype: float64
禁煙による体重の変化は3.37ポンドとなりました。
モデルの評価を行う
最後にモデルの評価をします。
evaluator = PropensityEvaluator(ipw)
results = evaluator.evaluate_simple(data.X, data.a, data.y, plots=None)
results.scores.prediction_scores
見慣れない指標について
- hinge : ヒンジ損失。(参考)
- matthews : マシューズの相関係数。-1から1。1に近いほど良い分類ができている。（参考）
- 0/1 : ???知ってる人教えてください
- brier : ブライアスコア。0から1で、0に近いほど正確であることを表す。(参考)
trainとvalidで表示されます。
from sklearn import metrics plots = ["roc_curve", "pr_curve", "weight_distribution", "calibration", "covariate_balance_love", "covariate_blance_slope"] metrics = {"roc_auc": metrics.roc_auc_score, "avg_precision": metrics.average_precision_score,} ipw = IPW(LogisticRegression(solver="liblinear")) evaluator = PropensityEvaluator(ipw) results = evaluator.evaluate_cv(data.X, data.a, data.y, plots=plots, metrics_to_evaluate=metrics)
- 左上：ROC曲線。AUCは0.7以上だと良いとされています。
- 右上：PrecisionとRecallのグラフ
- 左中：傾向スコアの分布。割当０は０に近いほど多く、割当１は１に近いほど多いと良いようです。
- 右中：キャリブレーションカーブ。横軸は傾向スコア。縦軸は割当１のデータの割合。x=yの直線に近いほど良い。
- 左下：共変量の平均差を標準誤差で割ったものの絶対値？？（ここがよくわかりませんでした）
- 右下：左下と同じ。
独自データに適用してみる
予め因果効果のわかっているデータを作成し、精度を検証したいと思います。
データの作成
年齢と性別によって目的変数yが変化するデータを作成します。
割当zは 年齢と性別からシグモイド関数によって決定しました。
from numpy.random import * import matplotlib.pyplot as plt import pandas as pd import numpy as np #データ数 num_data = 10000 #性別 gender = randint(0, 2, num_data) #年齢 age = randint(20, 81, num_data) #ノイズの生成,randnは標準正規分布 e_z = randn(num_data) #シグモイド関数に入れる部分 #これで割当を決める #ここを調節することで傾向スコアが満遍なく散らばるようにする z_base = 5*gender + age - 50 + 5*e_z #シグモイド関数 def sigmoid(x): return 1/(1+np.exp(-x)) #シグモイド関数を計算 z_prob = sigmoid(0.1*z_base) #割当の配列を作成 z = np.array([]) for i in range(num_data): tmp = np.random.choice([0,1], size=1, p=[1-z_prob[i], z_prob[i]])[0]#pは選ぶ確率 z = np.append(z, tmp)#配列の更新 #目的変数のノイズを生成 e_y = randn(num_data) #目的変数を決める #ここのzの係数が因果効果 y = 5000 + 10*gender + 30*age + 60*z + 10*e_y
割当zによる目的変数の変化は60にしました。
データフレームの作成
causallibで使えるようにデータを整形します。
#データフレームの作成 df = pd.DataFrame({'gender':gender, 'age': age, 'z':z, 'y':y}) df_x = df[['gender', 'age']] df_y = df['y'] df_z = df['z'] df.head()
| gender | age | z | y | |
|---|---|---|---|---|
| 0 | 0 | 76 | 1 | 7350.62 |
| 1 | 1 | 60 | 0 | 6784.8 |
| 2 | 1 | 34 | 0 | 6036.7 |
| 3 | 1 | 67 | 1 | 7089.36 |
| 4 | 0 | 50 | 0 | 6514.17 |
このようなデータを作成しました。
IPW法を適用する
これまで通りIPW法を適用します。
truncate_epsは0.01にしました。
#train learner = LogisticRegression(solver="liblinear") ipw = IPW(learner, truncate_eps=0.01) ipw.fit(df_x, df_z) #効果の確認 outcomes = ipw.estimate_population_outcome(df_x, df_z, df_y) #(割当１の目的変数 * 1/傾向スコア) - (割当0の目的変数 * 1/(1-傾向スコア))の合計 effect = ipw.estimate_effect(outcomes[1], outcomes[0], effect_types=["diff", "ratio"]) effect
出力
Fraction of values being truncated: 0.00000. diff 65.454931 ratio 1.010065 dtype: float64
割当の有無によるyの差は65.5と予測されました。
zの係数は60にしていたので、近い値ではありますね。
なお、同じ方法でデータを作成しても、予測値はばらつきました。
(55~80くらい)
何回か実行して平均をとった方が良いのかもしれないです。
評価
最後に評価を行います。
#評価 from sklearn import metrics plots=["roc_curve", "pr_curve", "weight_distribution", "calibration", "covariate_balance_love", "covariate_balance_slope"] metrics = {"roc_auc": metrics.roc_auc_score, "avg_precision": metrics.average_precision_score,} ipw = IPW(LogisticRegression(solver="liblinear"), truncate_eps=0.01) evaluator = PropensityEvaluator(ipw) results = evaluator.evaluate_cv(df_x, df_z, df_y, plots=plots, metrics_to_evaluate=metrics)
AUCが0.85と高く、 傾向スコアは理想的な形に分布しています。
データとしてはかなり予測しやすいものだったのではないでしょうか。
最後に
以上causallibによるIPW法でした！
因果推論初心者ですが、公式のサンプルがわかりやすかったのでなんとか理解できたかと。。。
評価方法やその基準についてわかっていないことが多いので、今後勉強していこうと思います。
