因果効果の推定！Rで実践 - 傾向スコア，マッチング，IPW推定量 -

はじめに

先日，岩波データサイエンスvol.3が発売されました．私も，「傾向スコアを用いたバント効果の推定−−ノーアウト1塁のバントは，得点確率を有意に高めるか!?」と題した記事を寄稿させていただきました．興味がある方は，是非読んでください＆感想・コメントなどもTwitterなどで教えてください．

岩波データサイエンス Vol.3

作者: 岩波データサイエンス刊行委員会
出版社/メーカー: 岩波書店
発売日: 2016/06/10
メディア: 単行本（ソフトカバー）
この商品を含むブログ (1件) を見る

今回のテーマ

今回は，上記の記事に関連して，上記の本で用いているIPW推定量の推定方法を実際に実践してみます．
説明していく過程で，読んでくださる方に，
・2つの群の平均の比較で因果効果をきちんと推定できない例
・傾向スコアの推定方法（ロジスティック回帰を今回は用います）とその周辺知識
・マッチングの計算方法
・IPW推定量の計算方法
をご紹介したいと思います．ただ，詳しく書きすぎると長くなってしまうので，詳細な部分については参考文献をあげるにとどめ，本の内容にかぶってしまうところは省きます（買って読んでください！お願いします）．では，始めましょう．

データセットをつくる

今回は，適当な処置d=0,1が結果変数Yに与える影響を推定するという問題を考えます（d=1のときは処置を行い，d=0のときは処置を行わなかったとします）．ただ，普通の例では面白くないので，
・実際には「処置dはYを減少させる効果がある」のに「単純な2群の比較をするとYを増加させる影響がある」という結論を導いてしまうという例を作成します．

こういう例は，実は比較的簡単に生成できます．まずは，共変量を作りましょう．ここでは，適当で構いません．とりあえず正規分布に従う共変量を4つほど作っておきます（先に必要な関数を読み込んでおきます．IPW計算用の関数と，ロジット変換などを手軽にするための関数です．）

> ######functions######
> expit = function(x){
+ 	rr = exp(x)/(1+exp(x))
+ 	return(rr)
+ }
> 
> logit = function(x){
+ 	rr = log(x/(1-x))
+ 	return(rr)
+ }
> 
> IPW = function(n,p,d,y){
+ 	numerator = sum(d*y/p)
+ 	denominator = sum(d/p)
+ 	r = numerator/denominator
+ 	return(r)
+ }
>#############
> set.seed(20160619)
> 
> #------------------------------
> # 標本数
> #------------------------------
> 
> n = 800
> 
> #------------------------------
> # 共変量
> #------------------------------
> 
> x1 = rnorm(n)
> x2 = rnorm(n)
> x3 = rnorm(n)
> x4 = rnorm(n)

ここからがポイントです．この共変量を使って，傾向スコアの真値に対するモデルを作成します．改めて述べておきますが，傾向スコアというのは，共変量が観測されたもとで，処置を受ける確率 $p(d=1|x)$ です．今回は，この確率に対して，以下のモデルを仮定します．

$e = p(d=1|x) = \cfrac{\exp(-0.8\times x_1 + 0.8 \times x_2 + 0.25 \times x_3 - 0.2 \times x_4)}{1+\exp(-0.8\times x_1 + 0.8 \times x_2 + 0.25 \times x_3 - 0.2 \times x_4)}$

ここで，この式を見てもらうと， $x_{2}, x_{3}$ が大きいほど処置を受ける確率が高く， $x_{1}, x_{4}$ が小さいほど処置を受ける確率が高くなるということがわかります（と，作っただけですが）．そして，この傾向スコアの真値に基づいて割り付け変数 $d$ を，ベルヌーイ分布から発生させます．

> #------------------------------
> # 傾向スコアの真値
> #------------------------------
> 
> p = expit(-0.8*x1 + 0.8*x2 + 0.25*x3 - 0.2*x4)
> 
> #------------------------------
> # 各標本の割り付け変数(d=1,0)
> #------------------------------
> 
> d = rbinom(n,1,p)

次に，共変量と割り付け変数を用いて結果変数yを作ります．ここで，傾向スコアが大きいときほど，Yが大きな値をとるようにします．具体的には， $x_{2}, x_{3}$ が大きいほど処置を受ける確率が高く， $x_{1}, x_{4}$ が小さいほど処置を受ける確率が高くなるわけですから，yをつくるときには， $x_{2}, x_{3}$ の係数を正に， $x_{1}, x_{4}$ の係数を負にするとよいです．こうすると，傾向スコアが高いものほど，yが大きくなりますから，処置の効果がマイナスでも，係数さえ調整すれば，処置群と対照群を平均で比較すると，処置群のほうがyの平均は大きくなります．

ここでは，yに対するモデルとして，
$y = 80 + (-7) \times d - 7 \times x_{1} + 6\times x_{2} + 6 \times x_3 - 4 \times x_{4} + \mbox{誤差}$
とします．

> #------------------------------
> # 結果変数 d=1なら結果は7下がる
> #------------------------------
> 
> y = 80 + (-7)*d - 7*x1 + 6*x2 + 6*x3 - 4*x4 + rnorm(n)
>

最後に，観測データを作りますが，実際に観測できる変数は $x_{1}, x_{2},x_{3}, x_{4}$ と， $d, y$ ですから（傾向スコアの真値は，観測されません），次のような観測データ行列を作っておきます．

> #------------------------------
> # 観測データ
> #------------------------------
> 
> data = cbind(x1,x2,x3,x4,d,y) 
> head(data)
             x1         x2         x3         x4 d        y
[1,] -1.2163990 -1.9305651 -0.6880511  0.1844445 0 72.59452
[2,]  0.2341120  0.1233155  0.9854147 -1.6495960 1 84.68886
[3,]  0.4495428  0.2482259  0.3206955 -0.2648300 1 75.93403
[4,]  0.6136386 -0.1667101 -0.9992488 -0.2649914 0 70.71651
[5,] -0.9403283  0.2755147  1.4383773  0.2877559 1 89.74142
[6,] -0.7335436  0.3948032 -1.3212992  0.6402797 1 70.53397

単純な2群の比較

では，実際に観測されたデータでd=1の群（処置群）とd=0の群（対照群）の結果変数を比較します．

> #------------------------------
> # 群ごとの平均の差
> #------------------------------
> 
> mean(y[d==1]) - mean(y[d==0])
[1] 3.256249

結果は正になりました．つまり，2群の平均を比較するということを行うと，処置の効果によってYの値は上昇するという結論を導くことになります．しかし，先ほどのyに対する式を確認しますと，
$y = 80 + (-7) \times d - 7 \times x_{1} + 6\times x_{2} + 6 \times x_3 - 4 \times x_{4} + \mbox{誤差}$
dの係数はマイナスなので，d=1である，つまり処置を受けた場合には，Yの値は減少するはずなのです．つまり，2群の比較によって得られる結果は，処置の因果効果を正しく推定できていません．

そこで，傾向スコアの登場です．ここからは，傾向スコアの推定と，それを用いた因果効果の推定法について述べます．

傾向スコアの推定と，傾向スコアを利用した推定量

まず，傾向スコア推定に焦点を当てます．傾向スコアとは，共変量が与えられた元で，処置群に割り付けられる確率でした．少し，わかりにくいかもしれませんが，言い換えると「共変量xの標本は，確率 $e(x)$ で，表が出るコインを振って，表なら処置群に，裏なら対照群に入る」ということです．このような問題は分布の言葉で書くと，ベルヌーイ分布を用いて，次のように書けます．
$(d|x) \sim \mbox{bernoulli}(e(x_{i}))$

気づいた方もいるかもしれませんが，これは一般化線形モデルの形に似ています．しかし， $e(x)$ が共変量 $x$ のどのような関数型で表現できるかは実際にはわからないので，このままでは，以下のようになります
$e(x) = f(x)$

ただ，関数型 $f$ がわからないと，傾向スコアを推定できません，そこで，ここでは，ロジスティック回帰モデルを選択するということにします．つまり，以下のように関数の形に仮定を置くことにします．
$\log\cfrac{e(x)}{1-e(x)} = x'\beta$

荒い言い方になりますが，一般的には傾向スコアが，どんな形で $x$ に依存するかはわからないので，試す必要がある場合もあります．また，どの変数 $x$ に依存しているのかもわからないという問題もあります．その場合も，とりあえず関係のありそうな変数を入れて，また対数変換，交互作用，指数変換などを用いて，傾向スコアのもっともらしいモデルを探します．

その際には，モデルのただしさを何で判断するの？という問題があります．これについては各分野の人が，「目安」を用意していますが，この基準さえ満たしておけば良いという判断基準はないので，あくまで目安であることに注意した方がよいと思います（実際，うまく傾向スコアを予測できる場合もありますが，そうではない場合もあります．データや分野に依存するのではないでしょうか？）．

では，観測されたデータを用いて，ロジスティック回帰を当てはめて傾向スコアを推定します．

> #------------------------------
> # 傾向スコアの推定
> #------------------------------
> 
> #観測データから傾向スコアの推定を行う．
> #傾向スコアは共変量が与えられたもとでの，d=1になる確率P(d=1|x)．
> #ただ，何もモデルを仮定しないままでは推定できないので，ここではロジスティック回帰を使う
> 
> #割り付け結果を，共変量で回帰する
> res = glm(d~x1+x2+x3+x4,data=as.data.frame(data),family=binomial(link="logit"))
> 
> summary(res)
Call:
glm(formula = d ~ x1 + x2 + x3 + x4, family = binomial(link = "logit"), 
    data = as.data.frame(data))

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.3775  -0.9835   0.4383   0.9341   2.5166  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  0.15718    0.08026   1.958   0.0502 .  
x1          -0.71110    0.08596  -8.272  < 2e-16 ***
x2           0.81237    0.09085   8.941  < 2e-16 ***
x3           0.36507    0.08447   4.322 1.55e-05 ***
x4          -0.33182    0.08230  -4.032 5.53e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1106.39  on 799  degrees of freedom
Residual deviance:  916.61  on 795  degrees of freedom
AIC: 926.61

Number of Fisher Scoring iterations: 4

> #ここでの傾向スコアは，ロジスティック回帰を行った時の確率の推定値（当てはめ値）
> hat.e = res$fitted
> head(hat.e)
        1         2         3         4         5         6 
0.6784036 0.8307323 0.7208275 0.7212536 0.6900221 0.4944932 
[f:id:tomoshige_n:20160619031457p:plain]

これで，傾向スコアはhat.eとして推定されたことになります．ここで，傾向スコアの推定値は次の通りです． $\hat\beta$ はロジスティック回帰で推定された各共変量の回帰係数です．
$\hat{e}_{i}(x_{i}) = \cfrac{\exp(x_{i}'\hat{\beta})}{1+\exp(x_{i}'\hat{\beta})}$ .

また，最初に設定した真のモデルがロジスティック回帰モデルの平均構造でしたから，そこから発生したデータに同じモデルを当てはめていますので，推定された回帰係数は，もちろん最初に設定した真のモデルの回帰係数に近くなります．

ここからは，処置の平均処置効果（処置の因果効果）を求めていきましょう．

マッチングによる推定量の構成

マッチングは，傾向スコアが同じ標本の，処置群の標本平均と，対照群の標本平均を計算して，その差を計算して，計算結果を群ごとに傾向スコアについて平均することによって，因果効果を推定するというのが基本の考え方です．では，なぜ傾向スコアが同じ標本同士を比較するのでしょうか？マッチングの本質はここにあります．

この理由を理解するためには，そもそもの因果効果を解析する際の仮定を思い出す必要があります．それがなんだったかというと「強く無視できる割り付け」というものでした．これは，「割り付けに関係する共変量をすべて固定してしまえば，割り付け変数dと，潜在結果変数の組 $(y_{1},y_{0}$ が独立になるというものです．この結果は，極めて重要な意味を持ちます．次の式をみてください．

$E(y_{1}-y_{0}) = E\left( E(y_{1} - y_{0} | x) \right) = E\left( E(y_{1} | x)\right) - E\left( E(y_{0} | x) \right) = E\left( E(y_{1} | d =1, x)\right) - E\left( E(y_{0} | x, d=0) \right)$

この式の左辺は，因果効果です．最後の等号は強く無視できる割り付けの仮定から成立します．この式の意味は，処置群と対照群で同じ共変量 $x$ を持つ標本の平均を計算した後に，その差を計算するという操作を，すべての $x$ で平均することができれば，平均処置効果は求めることができるという意味です．もっとわかりやすく（厳密さを欠いて言えば），「各共変量の値について，その値を持つ標本を処置群・対照群から選び出し，その標本の結果変数の差を計算するという操作を，すべての共変量の取りうる値について平均化すればいいんだ」と言っています．

これは，つまり，どちらの群に割り付けるのかを決める共変量 $x$ の値が等しい標本ごとの差を計算せよ！と言っています．つまり，その共変量を固定した箇所で見れば，標本はランダムに割り付けられていることになります．なので，これが処置の効果になるんです．しかし，すべての共変量で揃えるのは大変です．．．が，ここからがとても大事な結果です．実は上の式は，強く無視できる割り付けが成立する元では， $x$ のところを $e(x)$ （傾向スコア）に変えても成立するということが示されています（Rosembaum, and Rubin 1983）．すなわち，

$E(y_{1}-y_{0}) = E\left( E(y_{1} | e(x), d =1)\right) - E\left( E(y_{0} | e(x), d=0) \right)$

これはすなわち，傾向スコアが等しいところでの群平均の差を，傾向スコアで足しあわせて平均かしたものが処置の効果になると言っているのと同じことですね！そして，この考え方を利用するのがマッチングなわけです．

ただし，傾向スコアの値がぴったり同じという標本を見つけるのは難しいわけですね．そこで，傾向スコアを適当な離散区間に区切ります．今回は0.05幅に区切った結果を示します．すると，結果は-7.221となり，見事に因果効果-7に近い値が得られましたね！

しかし，ここで区間幅の決め方をどうするのかが恣意的になっていることが，ちょっと厄介ですね．実際，区切りを荒くしてみると結果は変わりますし，細かくしすぎると群ごとの平均が取れなくなるので，結果は不安定になります．これは困ったということで，別の方法が次のIPW推定量です．
（＊興味がある方は，他の区切り方も試してください．）

>#----------------------------------
># 傾向スコアマッチングによる因果効果推定
>#----------------------------------
> interval = seq(0,1,0.05) #machingを行う時の傾向スコアの幅を指定する
> maching_res = rep(0,10)
> 
> for(i in 1:(length(interval)-1)){
+ 	temp0 = y[(d==0)&((interval[i]<hat.e)&(hat.e<interval[i+1]))]
+ 	temp1 = y[(d==1)&((interval[i]<hat.e)&(hat.e<interval[i+1]))]
+ 	maching_res[i] = mean(temp1) - mean(temp0)
+ }
> 
> maching_res
 [1] -5.1711970        NaN -0.3080318 -4.2663056 -8.1488541 -7.9177035 -5.4518013
 [8] -6.9564249 -6.8824864 -7.6636229 -7.1455726 -7.2535686 -7.4277173 -5.7344403
[15] -6.7976385 -8.0770380 -6.7897113 -6.7753327 -8.8878558        NaN
> 
> mean(maching_res,na.rm=T)
[1] -6.536406

IPW推定量へと

最後にIPW推定量についてです．IPW推定量は，その定義から明らかなように，傾向スコアの逆数で重みをつけて因果効果を推定します．これにはどんな意味があるのでしょうか？また，1つの標本を考えます．この標本の共変量は $x$ であるとします．このとき，この標本が処置群に行く確率は80%で，対照群に行くか確率は20％であるとします．さらに，ここでは簡単のため処置群と対照群の人数はそれぞれ1000人であるとします．つまり，100人の同じ共変量 $x$ の人がいれば，期待値の意味で，80人は処置群へ，20人は対照群へ行くわけです．そうすると，共変量すべての標本について和をとるとき，処置群では，共変量 $x$ を持つ人の結果変数は80/1000 = 8%考慮され，対照群では20/1000=2%考慮されます．

何が言いたいかというと，共変量 $x$ を持つ人の1000人中に占める割合が異なってしまうために，両軍でそれぞれ平均を計算すると共変量のせいで差が生まれます．実際，結果変数が $x$ に依存している場合には，当たり前ですが，差が生まれます．

今回の共変量に依存する割り付けというのは次の問題のように捉えてもいいと思います．

今80人の子供がいます．身長は150cmか，140cmの人しかおらず，150cmの人の体重は50kg, 140cmの人の体重は40kgです．そして，50kgの人は確率0.8でクラス1へ，確率0.2でクラス2に入ります．一方で，40kgの人は確率0.8でクラス2へ，確率0.2でクラス1へ入ります．そして，体重が50kgと40kgの人はそれぞれ40人です．このとき，クラス1とクラス2の平均身長を比較すると異なる理由はなぜでしょうか？もしランダムに割り付けられたとすると，クラスの平均身長はどうなるでしょうか？

こう考えると，これは体重という共変量が，割り付け（クラス分け）と，結果変数（身長）に影響をするから，クラス1とクラス2の平均身長が異なるということになります．では，補正するためにはというと，クラス1では40kgの人に1/0.2の重みで身長を評価し，50kgの人を1/0.8で評価します．これは，割り付け確率の逆確率（つまり傾向スコアの逆重み付けです）です．一方で，クラス2では50kgの人に1/0.2の重みで身長を評価し，40kgの人を1/0.8で評価します．すると，クラスの平均身長の差の平均は0になることがわかります．なぜなら，最初の割り付けにおいて，身長150cmの人は0.8の確率でクラス1に入りますから，平均すると40*0.8 = 32人がクラス1に入ることになります．そう考えると，クラス1とクラス2の調整済み平均は
クラス1：(150/0.8×32 + 140/0.2×8)/40 = 145
クラス2：(150/0.2×8 + 140/0.8×32)/40 = 145
となり，割り付けをきちんと戻せば，平均に差がないことがわかります．

これと同じことを考えているのが，傾向スコアによるIPW推定量ということになります．つまり，傾向スコアの逆数による重み付けというのは，この共変量によって影響される人数比の偏りを補正する役目（つまりは，ランダムな割り付けにするための処置！！）があるわけですね．さて，この考え方からIPW推定量は，ここで書かれているようになります．

smrmkt.hatenablog.jp

> #-----------------------------
> # IPW推定量の計算
> #-----------------------------
> 
> #傾向スコアが求まったらIPW推定量を計算します．
> IPW.Treated = IPW(n,hat.e,d,y)
> IPW.Untreated = IPW(n,1-hat.e,1-d,y)
> 
> #因果効果
> IPW.Treated - IPW.Untreated
[1] -6.953027

近い値が得られました！めでたし！（ただ，set.seedを変えると．．．ちょっとざんない結果になることもありますので，ご了承を！ただ，それは推定量のばらつきなので，問題ありません．）

おわりに

ここまで，長い道のりをお付き合いくださりありがとうございました．
ここまでなんと，11,000字もあり，簡単な論文程度の文字数になってしまいました．
8時間もかかってしまいました．外明るくなってきてる気がする．．．
傾向スコアがどういうもので，IPW推定量・因果効果とは何かについて少しでも助けになれば幸いです．

岩波データサイエンスvol.3を片手に持って，
一つ一つのRコードの意味を岩波データサイエンスで確かめながら，
読んでいただけれいれば，さらに嬉しく思う次第です．

野球データの解析について

岩波データサイエンスに寄稿させていただいた野球データの解析では，
上記の議論に加えて，平滑化スプラインなど発展的な手法を用いて解析しています．
この辺りについてはサポートページに書かせていただきましたので参照ください．
ただし，傾向スコアや，IPW推定量の計算はここに書かせていただいた方法と大きくは変わりません．

・サポートページ
【特集】因果推論 ― 実世界のデータから因果を読む - 岩波データサイエンス