#はじめに
前回、『Uplift-Modelingで介入効果を最適化する』では、RCTによって得られたデータからある介入の効果を平均的に求めるのではなくて、ある特徴量を持った集団にのみ介入するという戦略を試しました。その結果、「男性向け」・「女性向け」プロモーションメールを平均的によかった方に全振りするのではなくて、約65%の人には「男性向け」メール、残りの人には「女性向け」メールを送る戦略がコンバージョンを最大化することを示しました。
その後、この問題をRCTとは言えない観察データへと拡張するための勉強をしていたのですが、その途中で傾向スコア分析という手法に行き着き、それもまとめておこうと思いました。
*本記事で使用したコードはこちら (propensity-score-basics.ipynb) からご覧いただけます。
#目次
#統計的因果推論の基礎
統計的因果推論は、ある介入の因果効果をできるだけ正しく推定したいというモチベーションで発展してきた分野です。例えば、あるサプリを飲んだ時に体重が減るのか否かやある広告を打った時にその商品の売上が上がるかどうかを厳密に推定したいということを議論します。このような問題においてよく用いられるのが、ルービンによる因果モデルです。
とりあえずある介入を2値変数 $Z$ ($Z$=1の時、介入有りとする。)、介入がなされた目的であるアウトカムを $Y$とします。ここでルービン因果モデルでは、$Y$ は単に1変数ではなくて$Z$=1の時のアウトカムである $Y_{1}$ と $Z$=0の時のアウトカムである$Y_{0}$の2変数が潜在すると考えます。
この時、個人レベルの介入効果は $Y_{1} - Y_{0}$ と表せますが、ある個人についてこのどちらかの変数は観測し得ないため、個人レベルの介入効果を通常知ることはできません。これが前回紹介した「因果推論における根本問題」でした。
ここで、個人レベルの介入効果の代わりに取り扱われて来たのが、平均因果効果(Average Treatment Effect)で、これは母集団における $Y_{1} - Y_{0}$ の期待値、すなわち、$ATE = E(Y_{1} - Y_{0}) = E(Y_{1}) - E(Y_{0})$ として表せます。
この $ATE$を求めるための有効な方法がRCTでした。RCTでは母集団を介入群 ($Z$= 1)と非介入群($Z= 0$)にランダムに割り当てることで、介入の有無以外の変数が介入群と非介入群で均質になることが期待されるため、RCTによって生じた群間差が $ATE$ の推定量として用いられます。
しかし、いつでもRCTを行うことができるとは限りません。例えば、タバコが健康に与える害の効果を推定する時は、ランダムにあなたはタバコを吸って、あなたはタバコを吸わないでというような割り当ては現実的ではないため、このよう場合はランダムに割り当てられたわけではないデータ（観察データ）を用いざるを得ません。
#RCTができない時の問題点
それでは、観察データを用いるときは何が問題になるのでしょうか。多くの場合、交絡因子の存在が因果効果の推定を難しくします。
例えばのちのCase Study②で取り上げるテレビCMがスマホアプリ利用時間等に与える因果効果を観察データから推定したいという時に、単に介入群と非介入群におけるスマホアプリ利用の程度の平均の差を求めてしまうと、CMは逆にアプリ利用を減少させてしまうという結果になります。これは、例えばCMを見ている人はテレビ視聴時間も長く、テレビ視聴時間が長い人は年齢層も高いためそもそものスマホアプリ利用頻度が低いと考えられ、これを考慮しないためにCMの効果が正しく推定されていないと考えることができます。
このように、介入の有無をランダムに割り当てることができない場合は、事前に存在する介入群と非介入群間に差が見られる特徴量（テレビCMの例では対象者の年齢）が割り当て確率 $P(Z=1)$ とアウトカム $Y$ の両方に影響を与える場合、この影響を除去してから介入の因果効果を推定する必要があるのです。
#傾向スコアとIPW推定量
このような交絡因子（介入の割り当て確率とアウトカムの両方に影響を与える特徴量）が存在する場合によく用いられる手法の一つに傾向スコア (Propensity Score)を用いた逆確率重み付け法 (IPW)があります。
傾向スコアは、説明変数ベクトル $x_i$ を持つある個人 $i$ が介入群に存在する確率 $e(x_i) = P(Z_i = 1 | x_i)$ として表せます。そして、この傾向スコア、すなわち、個人個人について本来ならこう割り当てられるべきであると考えられる確率を用いてランダムではない割り当てを補正するのです。
その補正の方法の一つがIPWという方法で、$Z=1$である時のアウトカムの期待値の推定量を、$\hat{E}(Y_1) := (\sum_{i=1}^{N}z_iy_i,/,e(x_i)),/,(\sum_{i=1}^{N}z_i,/,e(x_i))$ として、$Z=0$である時のアウトカムの期待値の推定量を、$\hat{E}(Y_0) := (\sum_{i=1}^{N}(1-z_i)y_i,/,(1-e(x_i))),/,(\sum_{i=1}^{N}(1-z_i),/,(1-e(x_i)))$ とすると、$\hat{E}(Y_1) - \hat{E}(Y_0)$が$ATE$ の一致推定量となることが知られています。(注)
(注)しかし、「強い意味での無視可能性」として、$(y_0, y_1) \perp z,|,\boldsymbol{x}$が成り立つという仮定が必要になります。つまり介入の割り当て変数$Z$とアウトカム$Y_0, Y_1$は、共変量$\boldsymbol{x}$ が与えられると独立という仮定で、これは「割り当て変数$Z$は共変量のみに依存し、結果変数には（直接的に）依存しない。」と言い換えることができます。例えば、英語教育プログラムの英語の能力に対する因果効果を推定したいというときに、「英語の能力が低いと英語教育プログラムに割り当てられていた」というような状況があってはならず、居住地や親の収入などの共変量$\boldsymbol{x}$によって、英語教育プログラムへの割り当てが説明できる必要があるということです。この条件が成り立っていることを確かめる方法の1つとして、以下2つのCase Studyでは、c統計量を用いた確認を行っています。
#Case Study① [2]
「Right Heart Catheterization Dataset」という心臓カテーテルに関するデータを用いて心臓カテーテルがその後の患者の180日以内死亡率に与える因果効果をIPWにより推定してみます。（以後、心臓カテーテルはRHCと表記します。）
# データの読み込み
data_df = pd.read_csv("http://biostat.mc.vanderbilt.edu/wiki/pub/Main/DataSets/rhc.csv")
data_df.head()
# 死亡率とRHCの有無のクロス集計
pd.crosstab(data_df.death, data_df.swang1)
# RHCの有無での死亡率の差を計算
(1486 / (698 + 1486)) - (2236 / (1315 + 2236))
------------------
0.05072 ...
------------------
これにより、RHCを受けていたほうが、約5%死亡率が高い結果となりました。しかし、このRHCの有無はランダムではないため、その他に観測されている変数を用いて患者ごとに傾向スコアを計算し、この結果を補正IPWにより補正してみます。
# 用いる説明変数群、詳細はデータセットの記述を参照
cols = ["cat1", "sex", "race", "edu", "income",
"resp", "card", "neuro", "gastr", "renal", "meta", "hema", "seps", "trauma", "ortho",
"das2d3pc", "dnr1", "ca", "surv2md1", "aps1", "scoma1", "wtkilo1", "temp1",
"resp1", "hrt1", "pafi1", "paco21", "ph1", "wblc1", "hema1", "sod1", "pot1", "crea1",
"bili1", "alb1", "cardiohx", "chfhx", "dementhx", "psychhx", "chrpulhx", "renalhx",
"liverhx", "gibledhx", "immunhx", "transhx", "amihx",
"age", "meanbp1"]
# 説明変数中のカテゴリカル変数
categorical_columns = ["cat1", "sex", "race", "edu", "income", "ca", "dnr1",
"resp", "card", "neuro", "gastr", "renal", "meta", "hema", "seps", "trauma", "ortho"]
# カテゴリカル変数のダミー化
X = data_df[cols]
dummy = pd.get_dummies(X[categorical_columns], drop_first=True)
X = pd.concat([X, dummy], axis=1)
X = X.drop(categorical_columns, axis=1)
# 切片の導入
X.loc[:, "Intercept"] = 1
# RHC有無のダミー変数
z1 = pd.get_dummies(data_df["swang1"])["RHC"]
# 目的変数
y = pd.get_dummies(data_df["death"])["Yes"]
# StatsModelsのLogitにより傾向スコアを推定
glm = sm.Logit(z1, X)
result = glm.fit()
ps = result.predict(X)
# c統計量としてAUCをを計算
fpr, tpr, thresholds = roc_curve(z1, ps)
auc(fpr, tpr)
# 傾向スコアのc統計量が0.80以上あることが望ましいとされているが、とりあえずは大丈夫そう
------------------
0.79631 ...
------------------
# IPWによりATEを推定
ipwe1 = sum((z1 * y) / ps) / sum(z1 / ps)
ipwe0 = sum(((1 - z1) * y) / (1 - ps)) / sum((1 - z1) / (1 - ps))
ATE = ipwe1 - ipwe0
ATE
------------------
0.05907 ...
------------------
このように、RHCに割り当てられるであろう確率である傾向スコアを用いて半年以内死亡率の推定量を補正してみたところ約1%差が広がる結果となりました。つまり、RHCは死亡率を上げているということになりますね...。参考にさせていただいた記事でもほぼ同じ結果となっていますので、大きな間違いはないと思いますが、この分野の知識を私は持っていないので、知り合いの医学部生にこの結果の妥当性を聞いてみようと思います。
#Case Study② [3]
次に、先ほども少し出てきたテレビCMがあるスマホアプリの利用率や利用回数、利用時間に与える因果効果を推定してみます。データセットはこちらを利用させていただきました。
Case Study②と同様の流れで進めてみます。
# データの読み込み
data_df = pd.read_csv('https://github.com/iwanami-datascience/vol3/raw/master/kato%26hoshino/q_data_x.csv')
data_df.head()
# アプリ利用ダミーとCM視聴有無のクロス集計
pd.crosstab(data_df.gamedummy, data_df.cm_dummy)
# CM視聴有無でのアプリ利用率の差を計算
(312 / (3832 + 312)) - (428 / (5428 + 428))
--------------
0.00220 ...
--------------
# CM視聴有無でのアプリ利用回数の差を計算
data_df[data_df.cm_dummy == 1].gamecount.mean() - data_df[data_df.cm_dummy == 0].gamecount.mean()
--------------
-1.48454 ...
--------------
# CM視聴有無でのアプリ利用時間の差を計算
data_df[data_df.cm_dummy == 1].gamesecond.mean() - data_df[data_df.cm_dummy == 0].gamesecond.mean()
--------------
-629.64...
--------------
これを見ると、CMを視聴してもアプリ利用率はほとんど変わらず、利用回数や利用時間については逆にCM視聴により減少しているという直感に反する結果が計算されています。しかしこれは、前述したように対象者の年齢に加え、居住地域などによる影響があるかもしれないため、これらを傾向スコアによるIPWで補正してATEを推定してみましょう。
# 説明変数
cols = ["age", "sex", "TVwatch_day", "marry_dummy", "child_dummy", "inc", "pmoney",
"area_kanto", "area_tokai", "area_keihanshin",
"job_dummy1", "job_dummy2", "job_dummy3", "job_dummy4", "job_dummy5", "job_dummy6",
"fam_str_dummy1", "fam_str_dummy2", "fam_str_dummy3", "fam_str_dummy4"]
X = data_df[cols].copy()
# 切片の導入
X.loc[:, "Intercept"] = 1
# CM視聴有無ダミー
z1 = data_df.cm_dummy
# 目的変数群（1:アプリ利用ダミー, 2:アプリ利用回数、3:アプリ利用時間）
y1 = data_df.gamedummy
y2 = data_df.gamecount
y3 = data_df.gamesecond
# StatsModelsのLogitにより傾向スコアを推定
glm = sm.Logit(z1, X)
result = glm.fit()
ps = result.predict(X)
ps
-----------
0 0.046217
1 0.255983
2 0.177427
3 0.227593
4 0.242373
-----------
# c統計量としてAUCをを計算
fpr, tpr, thresholds = roc_curve(z1, ps)
auc(fpr, tpr)
# 今回もとりあえずは大丈夫そう
-----------
0.79170 ...
-----------
# IPWによりアプリ利用ダミーへのATEを推定
ipwe11 = sum((z1 * y1) / ps) / sum(z1 / ps) # Treated
ipwe10 = sum(((1 - z1) * y1) / (1 - ps)) / sum((1 - z1) / (1 - ps)) # Control
ATE1 = ipwe11 - ipwe10
ATE1
-----------
0.03231 ...
-----------
# IPWによりアプリ利用回数へのATEを推定
ipwe21 = sum((z1 * y2) / ps) / sum(z1 / ps) # Treated
ipwe20 = sum(((1 - z1) * y2) / (1 - ps)) / sum((1 - z1) / (1 - ps)) # Control
ATE2 = ipwe21 - ipwe20
ATE2
-----------
5.3490 ...
-----------
# IPWによりアプリ利用時間へのATEを推定
ipwe31 = sum((z1 * y3) / ps) / sum(z1 / ps)
ipwe30 = sum(((1 - z1) * y3) / (1 - ps)) / sum((1 - z1) / (1 - ps))
ATE3 = ipwe31 - ipwe30
ATE3
-----------
1513.69 ...
-----------
傾向スコアを用いた補正により、今度はCM視聴がアプリ利用率を平均3.2%、アプリ利用時間を平均5.3回、アプリ利用時間を平均1513秒押し上げるという直感に整合的な結果が推定できました。また、これらの結果は今回参考にさせていただきた本の結果ともだいたい同じとなっています。（誤差はダミー変数のベースラインを何にする等の違いかと思われます。）
#さいごに
今回は、傾向スコアによりRCTではない観察データから因果効果を推定する方法であるIPWという手法を紹介し、2つのCase Studyで実験をしてみました。傾向スコアとしては今後IPWによる推定量が$ATE$の一致推定量となる証明を追ってみたり、推定値の信頼区間を計算する方法を勉強しようと思っています。しかし、直近では、Case Study②で使用したデータでUplift Modelingを試し、有効な介入戦略を立ててみようと思います。
#参考
[1] 統計的因果推論(2): 傾向スコア(Propensity Score)の初歩をRで実践してみる
[2] 傾向スコア（Propensity score）をRで実践 マッチングとIPWの結果を比較
[3] 加藤諒, 星野宗宏 因果効果推定の応用 CM接触の因果効果と調整効果 岩波データサイエンスVol.3 岩波書店