About connecting the dots.

data science related trivial things

調査観察データにおける因果推論(4) - Rで傾向スコアを出す際の共変量選択基準

はじめに

前回は細かい理屈をすっ飛ばして,全変数を共変量として突っ込んだロジスティック回帰を実行しましたが,当然そんなやり方はほめられたものではないわけです.今回は,この傾向スコア算出に使う共変量をどのような基準で選ぶかについてみてきます.

変数同士の関係性

共変量候補と,独立変数,従属変数の関係性は以下の3通りが考えられます.

  1. 処置前変数であり,従属変数に先行する変数
  2. 処置後変数であり,従属変数に先行する変数
  3. 処置後変数であり,従属変数の後になる変数

結論からいうと,1番目の変数群は使用するべきで,2番目の変数群を使うかどうかは目的に依存する,そして3番目の変数群は使用するべきではない,ということになります.

1番目の変数群については,処置変数にも従属変数にも時系列的に先行しているわけですから,論理的に当該変数が両者に影響していると考えてもなんらおかしいことはありません.lalondeデータについては,「年齢,教育年数,黒人,ヒスパニック,既婚,高卒以上,74年の年収,75年の年収,74年に収入ゼロ,75年に収入ゼロ」の10個の変数がありましたが,そのすべてがこのカテゴリに属すると考えられるので,文句なく共変量として用いていいことになります.

2番目の変数群については,lalondeの例でいうなら「再就職先の企業規模」という変数があったとしたら,このカテゴリに属すると考えられます.再就職は処置(=職業訓練プログラム)の後に生じるものである一方,賃金は就職した結果得られるものですので,従属変数(=78年の年収)より時間的に先行するものになります.ここで「再就職先の企業規模」を共変量に入れてしまうと,職業訓練プログラムを受けた群と受けていない群の間で,再就職先の企業規模を等しくした場合の年収の差分を推定することになります.しかし職業訓練プログラムを受けたことによって,そもそも就職可能な企業規模が大きくなる,という因果が成り立つと考えられるため,これを共変量に含めてしまうと,職業訓練プログラムの効果の一部(=企業規模によるベース賃金の差分)が因果効果から取り除かれてしまうため,結果として処置の因果効果を過小推定してしまうことになります.なのでこの場合は,共変量に含めない方が(論理的に)正しく因果効果を推定できることになります.

そして3番目の変数群は,因果効果の意味を考えれば.時間的に後に生じる変数を共変量として使用するべきではありません.

共変量の妥当性確認

さて,いくつかよさそうな共変量を選択したとして,それがどの程度妥当か,というのを確認する必要があります.これまでの説明では詳細を省いてきましたが,実は傾向スコアによって因果効果の推定を行うためには,「強く無視できる割り当て」という条件が成立している必要があります.

「強く無視できる割り当て」とは

これは「処置群と非処置群のどちらに割り当てられるかは,共変量の値に依存し,従属変数による割り当てへの影響はあくまで「共変量と従属変数の関係」を通じてのみ,間接的に存在している」という仮定になります.逆にいえば,「共変量だけから,処置群と非処置群のどちらに当てはまるかが予測でき,従属変数から処置の有無への直接的な因果関係は存在しない」ということになります.

ですが実際問題として,この「強く無視できる割り当て」条件を直接的に検証することは不可能です.そこで間接的に検証する方法として,以下の3つがあげられます.

  1. 共変量によって割り当てが説明されていることを示す
  2. 共変量そのものを調整する
  3. 理論を精緻化する

このうち1番目の条件は,非常にわかりやすく共変量から処置の有無を予測するロジスティック回帰の当てはまりが良ければOK,ということです.この手法について,下でRのサンプルを用いて説明します.2番目については,共変量を条件づけたときに,実際に群間の共変量の差が無くなっているかどうかで検証可能です.同じようにRのサンプルで説明したいと思います.そして3番目については,学術研究ならいざ知らず,実務レイヤーにおいては理論的な予測というものはほぼ存在しないと思われるので,あまり考慮する必要はないでしょう.

共変量による割り当て説明の確認

共変量選択における,ロジスティック回帰モデルの当てはまりの基準としてよく使われるのは,c統計量というものです.これはROC曲線(詳細は奥村先生の説明を参照してください)の下側面積と等しくなるそうです.Rでc統計量を算出するには,{rms}パッケージを使う必要があります.lrm()メソッドでロジスティック回帰を行った結果に,c統計量が格納されます.具体的なコードは以下の通りです.statsの中には様々な統計量が一緒くたに格納されているので,その6番目の値になります.この値が0.8以上あればOKというのが医学系研究におけるデファクトのようです.

ret <- lrm(treat~age+educ+black+hisp+married+nodegr+re74+re75+u74+u75, data=lalonde)
ret$stats
ret$stats[6]
> ret$stats
         Obs    Max Deriv   Model L.R.         d.f.            P            C          Dxy        Gamma        Tau-a           R2        Brier            g           gr           gp 
4.450000e+02 1.802618e-09 1.994457e+01 1.000000e+01 2.978132e-02 6.220166e-01 2.440333e-01 2.457911e-01 1.188177e-01 5.900912e-02 2.322175e-01 4.952981e-01 1.640987e+00 1.166198e-01 
> ret$stats[6]
        C 
0.6220166 

なお,ROCをプロットすることも{Epi}パッケージを使えば簡単にできます.右下のarea under the curveというところに,0.622という値が示されているのがわかると思います.まぁ要するに,ここで使ってきたlalondeデータは,共変量選択の当てはまりとしてはダメなわけですね.

Epi::ROC(form=treat~age+educ+black+hisp+married+nodegr+re74+re75+u74+u75,
         data=lalonde, plot="ROC")

f:id:SAM:20130928195732j:plain

共変量そのものを調整する

共変量で調整した結果の処置群非処置群の共変量については,{Matching}パッケージで確認できます.どの変数についても,マッチングの結果差が減少しているのは見て取れるかと思います.

mout <- Match(Y=lalonde$re78, Tr=lalonde$treat, X=ps, estimand="ATE")
ret <- MatchBalance(treat~age+educ+black+hisp+married+nodegr+re74+re75+u74+u75,
                    data=lalonde, match.out=mout, nboots=1000)
> mout <- Match(Y=lalonde$re78, Tr=lalonde$treat, X=ps, estimand="ATE")
> ret <- MatchBalance(treat~age+educ+black+hisp+married+nodegr+re74+re75+u74+u75,
+                     data=lalonde, match.out=mout, nboots=1000)

***** (V1) age *****
                       Before Matching 	 	 After Matching
mean treatment........     25.816 	 	     25.704 
mean control..........     25.054 	 	     25.393 
std mean diff.........     10.655 	 	     4.4037 

mean raw eQQ diff.....    0.94054 	 	     0.8378 
med  raw eQQ diff.....          1 	 	          1 
max  raw eQQ diff.....          7 	 	          9 

mean eCDF diff........   0.025364 	 	   0.023314 
med  eCDF diff........   0.022193 	 	   0.016925 
max  eCDF diff........   0.065177 	 	    0.06347 

var ratio (Tr/Co).....     1.0278 	 	    0.91249 
T-test p-value........    0.26594 	 	    0.46162 
KS Bootstrap p-value..      0.527 	 	      0.055 
KS Naive p-value......     0.7481 	 	    0.11496 
KS Statistic..........   0.065177 	 	    0.06347 


***** (V2) educ *****
                       Before Matching 	 	 After Matching
mean treatment........     10.346 	 	      10.26 
mean control..........     10.088 	 	      10.27 
std mean diff.........     12.806 	 	   -0.54649 

mean raw eQQ diff.....    0.40541 	 	    0.12271 
med  raw eQQ diff.....          0 	 	          0 
max  raw eQQ diff.....          2 	 	          2 

mean eCDF diff........   0.028698 	 	  0.0087649 
med  eCDF diff........   0.012682 	 	  0.0056417 
max  eCDF diff........    0.12651 	 	   0.053597 

var ratio (Tr/Co).....     1.5513 	 	     1.0235 
T-test p-value........    0.15017 	 	     0.9182 
KS Bootstrap p-value..      0.014 	 	      0.064 
KS Naive p-value......   0.062873 	 	    0.26035 
KS Statistic..........    0.12651 	 	   0.053597 


***** (V3) black *****
                       Before Matching 	 	 After Matching
mean treatment........    0.84324 	 	    0.83521 
mean control..........    0.82692 	 	    0.83933 
std mean diff.........     4.4767 	 	    -1.1092 

mean raw eQQ diff.....   0.016216 	 	  0.0056417 
med  raw eQQ diff.....          0 	 	          0 
max  raw eQQ diff.....          1 	 	          1 

mean eCDF diff........  0.0081601 	 	  0.0028209 
med  eCDF diff........  0.0081601 	 	  0.0028209 
max  eCDF diff........    0.01632 	 	  0.0056417 

var ratio (Tr/Co).....    0.92503 	 	     1.0206 
T-test p-value........    0.64736 	 	    0.84586 


***** (V4) hisp *****
                       Before Matching 	 	 After Matching
mean treatment........   0.059459 	 	    0.07191 
mean control..........    0.10769 	 	   0.083146 
std mean diff.........    -20.341 	 	    -4.3444 

mean raw eQQ diff.....   0.048649 	 	  0.0098731 
med  raw eQQ diff.....          0 	 	          0 
max  raw eQQ diff.....          1 	 	          1 

mean eCDF diff........   0.024116 	 	  0.0049365 
med  eCDF diff........   0.024116 	 	  0.0049365 
max  eCDF diff........   0.048233 	 	  0.0098731 

var ratio (Tr/Co).....    0.58288 	 	    0.87546 
T-test p-value........   0.064043 	 	     0.3532 


***** (V5) married *****
                       Before Matching 	 	 After Matching
mean treatment........    0.18919 	 	    0.14719 
mean control..........    0.15385 	 	     0.1603 
std mean diff.........     8.9995 	 	    -3.6957 

mean raw eQQ diff.....   0.037838 	 	   0.019746 
med  raw eQQ diff.....          0 	 	          0 
max  raw eQQ diff.....          1 	 	          1 

mean eCDF diff........   0.017672 	 	  0.0098731 
med  eCDF diff........   0.017672 	 	  0.0098731 
max  eCDF diff........   0.035343 	 	   0.019746 

var ratio (Tr/Co).....     1.1802 	 	    0.93256 
T-test p-value........    0.33425 	 	    0.57205 


***** (V6) nodegr *****
                       Before Matching 	 	 After Matching
mean treatment........    0.70811 	 	    0.78352 
mean control..........    0.83462 	 	     0.7618 
std mean diff.........    -27.751 	 	     5.2686 

mean raw eQQ diff.....    0.12432 	 	  0.0056417 
med  raw eQQ diff.....          0 	 	          0 
max  raw eQQ diff.....          1 	 	          1 

mean eCDF diff........   0.063254 	 	  0.0028209 
med  eCDF diff........   0.063254 	 	  0.0028209 
max  eCDF diff........    0.12651 	 	  0.0056417 

var ratio (Tr/Co).....     1.4998 	 	    0.93472 
T-test p-value........  0.0020368 	 	     0.1463 


***** (V7) re74 *****
                       Before Matching 	 	 After Matching
mean treatment........     2095.6 	 	     1847.3 
mean control..........       2107 	 	     2047.4 
std mean diff.........   -0.23437 	 	    -4.1582 

mean raw eQQ diff.....     487.98 	 	     428.11 
med  raw eQQ diff.....          0 	 	          0 
max  raw eQQ diff.....       8413 	 	       8413 

mean eCDF diff........   0.019223 	 	   0.022367 
med  eCDF diff........     0.0158 	 	   0.023977 
max  eCDF diff........   0.047089 	 	   0.046544 

var ratio (Tr/Co).....     0.7381 	 	    0.79328 
T-test p-value........    0.98186 	 	    0.55065 
KS Bootstrap p-value..      0.561 	 	      0.043 
KS Naive p-value......    0.97023 	 	    0.42621 
KS Statistic..........   0.047089 	 	   0.046544 


***** (V8) re75 *****
                       Before Matching 	 	 After Matching
mean treatment........     1532.1 	 	     1200.7 
mean control..........     1266.9 	 	     1267.6 
std mean diff.........     8.2363 	 	    -2.4343 

mean raw eQQ diff.....     367.61 	 	     204.63 
med  raw eQQ diff.....          0 	 	          0 
max  raw eQQ diff.....     2110.2 	 	     5055.8 

mean eCDF diff........   0.050834 	 	   0.010505 
med  eCDF diff........   0.061954 	 	  0.0084626 
max  eCDF diff........    0.10748 	 	    0.03385 

var ratio (Tr/Co).....     1.0763 	 	    0.77337 
T-test p-value........    0.38527 	 	    0.72357 
KS Bootstrap p-value..      0.051 	 	      0.302 
KS Naive p-value......    0.16449 	 	    0.81133 
KS Statistic..........    0.10748 	 	    0.03385 


***** (V9) u74 *****
                       Before Matching 	 	 After Matching
mean treatment........    0.70811 	 	    0.74007 
mean control..........       0.75 	 	    0.74367 
std mean diff.........    -9.1895 	 	   -0.81886 

mean raw eQQ diff.....   0.037838 	 	   0.035261 
med  raw eQQ diff.....          0 	 	          0 
max  raw eQQ diff.....          1 	 	          1 

mean eCDF diff........   0.020946 	 	    0.01763 
med  eCDF diff........   0.020946 	 	    0.01763 
max  eCDF diff........   0.041892 	 	   0.035261 

var ratio (Tr/Co).....     1.1041 	 	     1.0091 
T-test p-value........    0.33033 	 	    0.88645 


***** (V10) u75 *****
                       Before Matching 	 	 After Matching
mean treatment........        0.6 	 	    0.66929 
mean control..........    0.68462 	 	    0.68075 
std mean diff.........    -17.225 	 	    -2.4333 

mean raw eQQ diff.....   0.081081 	 	   0.016925 
med  raw eQQ diff.....          0 	 	          0 
max  raw eQQ diff.....          1 	 	          1 

mean eCDF diff........   0.042308 	 	  0.0084626 
med  eCDF diff........   0.042308 	 	  0.0084626 
max  eCDF diff........   0.084615 	 	   0.016925 

var ratio (Tr/Co).....     1.1133 	 	     1.0185 
T-test p-value........   0.068031 	 	    0.63289 


Before Matching Minimum p.value: 0.0020368 
Variable Name(s): nodegr  Number(s): 6 

After Matching Minimum p.value: 0.043 
Variable Name(s): re74  Number(s): 7 

以上,共変量選択に関するあれこれをまとめました.