[Stan] 傾向スコアを使ったモデル
[1]統計的因果推論（傾向スコア）の勉強会資料をアプしてみた（web）
[2]傾向スコアを用いた共変量調整による因果効果の推定と臨床医学・疫学・薬学・公衆衛生分野での応用について（pdf file）
[3]傾向スコア解析法による因果効果の推定と調査データの調整について（web）
[4]調査観察データにおける因果推論(3) - Rによる傾向スコア，IPW推定量，二重にロバストな推定量の算出（web）
概略は以下です。
例えば喫煙の体への影響（例えばガンになったかどうか）を知りたいとします。被験者のデータ（性別・年齢・学歴・ストレスがあるか・飲酒するかどうかなど）とガンかどうかのデータがあるとします。この場合、喫煙する人は男性で年齢が高めでストレスがあって飲酒をする人がとても多いため、通常の解析では喫煙自体がその人のガンの発症に影響があるか（因果効果）がうまく推定できない場合があります。このように無作為化試験が難しい/機能していない場合に傾向スコアを使った解析が有用な様子です。
具体的には、被験者の性質からロジスティック回帰分析によって喫煙者群か非喫煙者群のどちらに割り付けられるかの確率を推定します。これを傾向スコア（propensity score）と呼びます。実際には割り付けの確率が出力されるもの（ランダムフォレストとかSVMとか）なら何を使ってもよいと思います。次にこの傾向スコアを使って重みづけを行って解析します。解析の仕方は色々ある様子です。詳しくは上記のリンク先を見てください。
個人的には、傾向スコアを使った解析は、偏りや欠測のあるデータにおいて、どのデータのあてはまりをよくしたいかの重み付けを変える一つの枠組みであると感じました。
以下では傾向スコアを題材にする際によく使われるデータセットを使ってStanで解析します。Rの{Matching}パッケージのlalondeデータです（445人分データ）。これは1976年に実施された米国職業訓練プログラムを受講した群としていない群（treat=1 or 0）で、訓練実施後のre78（1978年の年収）にどの程度違いが出たかというデータになります。他の説明変数として、age（年齢）・educ（教育年数）・black（黒人か）・hisp（ヒスパニックか）・married（既婚か）・nodegr（高卒以上か）・re74（74年の年収）・re75（75年の年収）・u74（74年に収入ゼロか）・u75（75年に収入ゼロか） の10個あります。
まずは散布図行列を見てみましょう。
当然、re74・re75・u74・u75の間には相関がありますね。u74・u75を使うのは疑義があると思いますが今はつっこみません。またre78の分布は明らかに0が多く、Zero-Inflated XXX な分布です。このデータを正規分布を使った線形モデルで行うのは相当突っ込みどころがあると思うのですが…今回は説明を傾向スコアにしぼりたいのでそれもスルーします！詳しくはこちらの記事を参照してください。
このデータに対してRで単純に線形回帰（lm）を行うと以下の結果が得られます。
Call:
lm(formula = re78 ~ age + educ + black + hisp + married + nodegr +
re74 + re75 + u74 + u75 + treat, data = lalonde)
Residuals:
Min 1Q Median 3Q Max
-9612 -4355 -1572 3054 53119
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 256.6702 3521.6821 0.07 0.9419
age 53.5709 45.8063 1.17 0.2428
educ 400.7703 228.8217 1.75 0.0806
black -2037.3328 1173.7163 -1.74 0.0833
hisp 425.8185 1564.5533 0.27 0.7856
married -146.3292 882.2569 -0.17 0.8683
nodegr -15.1794 1005.8851 -0.02 0.9880
re74 0.1234 0.0878 1.40 0.1608
re75 0.0197 0.1503 0.13 0.8955
u74 1380.2849 1187.9168 1.16 0.2459
u75 -1071.2155 1024.8779 -1.05 0.2965
treat 1670.7095 641.1323 2.61 0.0095
Residual standard error: 6520 on 433 degrees of freedom
Multiple R-squared: 0.0582, Adjusted R-squared: 0.0343
F-statistic: 2.43 on 11 and 433 DF, p-value: 0.00597
しかしながら、上記[2]によるとオススメできない、とのこと。該当部の引用は以下の通り。
共分散分析的手法の最大の問題点は，従属変数と共変量の関係をモデル化する必要があるということである．このモデル化では回帰関数を間違って指定すると（例えば二次関数の関係があるのに線形と仮定すると）誤った結果が導かれることなどがよく知られている．また，共分散分析的手法で推定することができる回帰係数自体は因果効果と等しくはならない．
というわけで傾向スコアを使ってみます。僕が思いついた範囲で通常の傾向スコアを使ったフローの不満点としては、傾向スコア自体にも確率的な変動があるはずなのに、そこは固定値として求めて最後の因果効果だけ確率的な変動幅を検討する点です。また傾向スコアを求めるところと傾向スコアを使って重み付けを行ってフィッティングするところは同時に推定するべきではないかと思いました。よってStanでの実装は以下になりました。
・19行目: e[n]が傾向スコアになります。
・20行目: 傾向スコアを使った重み付けです。割り付けとしてはGr=1（職業訓練を受けた）だけど説明変数からはGr=1にはいなさそうな人のデータ（Y）ほど当てはまってほしいので、このような傾向スコアの逆数が重みになっています。Gr=0の場合も同様です。逆数を使うことは理論的な背景もあるようですが、MCMC計算では分母が0や1に近いところでは不安定になりやすくちょっとつらいです。
・26行目: 各被験者のYをGrであてはめているのですが、それぞれの対数尤度の重みを傾向スコアで調整しています。この部分は上記の文献[2]のp.234の「（6）重み付け M推定量」に沿って実装しています。星野本ですとp.79の「3.4 一般的な周辺パラメトリックモデルの推定」の部分になります。
・43-45行目: IPW推定量を計算して因果効果を推定しています。
キックするRコードは以下になります。
・7,11,12行目: 適度にスケーリングしてあります。これをやっておかないとlogistic回帰の時にstanの実行時エラーが頻出します。
・8行目: Grは割り付け（グループ）のインデックスになります。
以下は結果になります。計算時間は1chainあたり30秒ぐらいでした。スケーリングをもとに戻すには beta_y や causal_effect は10000倍すればOKです。causal_effectの推定幅は思った以上に狭いですね…。
| mean | se_mean | sd | X2.5. | X25. | X50. | X75. | X97.5. | n_eff | Rhat | |
|---|---|---|---|---|---|---|---|---|---|---|
| alpha_x | -0.703 | 0.016 | 0.464 | -1.566 | -1.016 | -0.724 | -0.409 | 0.283 | 830 | 1.00 |
| beta_x[1] | 0.012 | 0.000 | 0.007 | 0.000 | 0.008 | 0.012 | 0.017 | 0.025 | 1500 | 1.00 |
| beta_x[2] | 0.035 | 0.001 | 0.033 | -0.030 | 0.013 | 0.036 | 0.057 | 0.099 | 824 | 1.00 |
| beta_x[3] | 0.088 | 0.006 | 0.174 | -0.268 | -0.027 | 0.094 | 0.203 | 0.419 | 846 | 1.00 |
| beta_x[4] | -0.163 | 0.008 | 0.231 | -0.605 | -0.329 | -0.160 | 0.002 | 0.282 | 925 | 1.00 |
| beta_x[5] | 0.109 | 0.003 | 0.127 | -0.133 | 0.021 | 0.107 | 0.192 | 0.354 | 1500 | 1.00 |
| beta_x[6] | -0.404 | 0.004 | 0.133 | -0.679 | -0.492 | -0.398 | -0.317 | -0.143 | 885 | 1.01 |
| beta_x[7] | 0.024 | 0.004 | 0.113 | -0.204 | -0.052 | 0.026 | 0.101 | 0.243 | 955 | 1.00 |
| beta_x[8] | 0.142 | 0.008 | 0.228 | -0.312 | -0.008 | 0.139 | 0.292 | 0.597 | 867 | 1.00 |
| beta_x[9] | 0.675 | 0.005 | 0.151 | 0.377 | 0.576 | 0.673 | 0.776 | 0.971 | 935 | 1.00 |
| beta_x[10] | -0.562 | 0.004 | 0.124 | -0.799 | -0.649 | -0.565 | -0.478 | -0.327 | 837 | 1.00 |
| alpha_y | 0.446 | 0.001 | 0.030 | 0.386 | 0.426 | 0.446 | 0.467 | 0.502 | 1198 | 1.00 |
| beta_y | 0.157 | 0.001 | 0.045 | 0.073 | 0.125 | 0.153 | 0.189 | 0.246 | 1296 | 1.00 |
| s_y | 0.639 | 0.000 | 0.016 | 0.609 | 0.628 | 0.640 | 0.650 | 0.670 | 1500 | 1.00 |
| (eは略) | … | … | … | … | … | … | … | … | … | … |
| ipw0 | 0.448 | 0.000 | 0.003 | 0.441 | 0.446 | 0.448 | 0.450 | 0.455 | 1257 | 1.00 |
| ipw1 | 0.602 | 0.000 | 0.005 | 0.593 | 0.599 | 0.602 | 0.605 | 0.611 | 1500 | 1.00 |
| causal_effect | 0.154 | 0.000 | 0.006 | 0.142 | 0.150 | 0.154 | 0.158 | 0.166 | 1500 | 1.00 |
| lp__ | -1169 | 0.098 | 2.575 | -1175 | -1170 | -1169 | -1167 | -1165 | 688 | 1.00 |
3グループ以上の時はsoftmax()を使って各グループへの所属確率を求めて、各グループの頂点ベクトル（(1,0,0)とか）へのユークリッド距離に応じて重み付けしたら拡張できるように思いました。
- 関連記事