調査観察データにおける因果推論(4) - Rで傾向スコアを出す際の共変量選択基準
目次
はじめに
前回は細かい理屈をすっ飛ばして，全変数を共変量として突っ込んだロジスティック回帰を実行しましたが，当然そんなやり方はほめられたものではないわけです．今回は，この傾向スコア算出に使う共変量をどのような基準で選ぶかについてみてきます．
変数同士の関係性
共変量候補と，独立変数，従属変数の関係性は以下の3通りが考えられます．
- 処置前変数であり，従属変数に先行する変数
- 処置後変数であり，従属変数に先行する変数
- 処置後変数であり，従属変数の後になる変数
結論からいうと，1番目の変数群は使用するべきで，2番目の変数群を使うかどうかは目的に依存する，そして3番目の変数群は使用するべきではない，ということになります．
1番目の変数群については，処置変数にも従属変数にも時系列的に先行しているわけですから，論理的に当該変数が両者に影響していると考えてもなんらおかしいことはありません．lalondeデータについては，「年齢，教育年数，黒人，ヒスパニック，既婚，高卒以上，74年の年収，75年の年収，74年に収入ゼロ，75年に収入ゼロ」の10個の変数がありましたが，そのすべてがこのカテゴリに属すると考えられるので，文句なく共変量として用いていいことになります．
2番目の変数群については，lalondeの例でいうなら「再就職先の企業規模」という変数があったとしたら，このカテゴリに属すると考えられます．再就職は処置（=職業訓練プログラム）の後に生じるものである一方，賃金は就職した結果得られるものですので，従属変数（=78年の年収）より時間的に先行するものになります．ここで「再就職先の企業規模」を共変量に入れてしまうと，職業訓練プログラムを受けた群と受けていない群の間で，再就職先の企業規模を等しくした場合の年収の差分を推定することになります．しかし職業訓練プログラムを受けたことによって，そもそも就職可能な企業規模が大きくなる，という因果が成り立つと考えられるため，これを共変量に含めてしまうと，職業訓練プログラムの効果の一部（=企業規模によるベース賃金の差分）が因果効果から取り除かれてしまうため，結果として処置の因果効果を過小推定してしまうことになります．なのでこの場合は，共変量に含めない方が（論理的に）正しく因果効果を推定できることになります．
そして3番目の変数群は，因果効果の意味を考えれば．時間的に後に生じる変数を共変量として使用するべきではありません．
共変量の妥当性確認
さて，いくつかよさそうな共変量を選択したとして，それがどの程度妥当か，というのを確認する必要があります．これまでの説明では詳細を省いてきましたが，実は傾向スコアによって因果効果の推定を行うためには，「強く無視できる割り当て」という条件が成立している必要があります．
「強く無視できる割り当て」とは
これは「処置群と非処置群のどちらに割り当てられるかは，共変量の値に依存し，従属変数による割り当てへの影響はあくまで「共変量と従属変数の関係」を通じてのみ，間接的に存在している」という仮定になります．逆にいえば，「共変量だけから，処置群と非処置群のどちらに当てはまるかが予測でき，従属変数から処置の有無への直接的な因果関係は存在しない」ということになります．
ですが実際問題として，この「強く無視できる割り当て」条件を直接的に検証することは不可能です．そこで間接的に検証する方法として，以下の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")
共変量そのものを調整する
共変量で調整した結果の処置群非処置群の共変量については，{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
以上，共変量選択に関するあれこれをまとめました．
