岩波書店の確率と情報の科学シリーズの星野崇宏 著「調査観測データの統計科学」を読んでいき、まとめていってます。
時折、サンプルデータを利用して実装検証していきます。
今回はやっと第3章「セミパラメトリック解析」で傾向スコアの因果効果についてです。
傾向スコアとは
傾向スコアとはローゼンバウムとルービンにより提案された、無作為割り当てが不可能な場合の因果効果を推定する手法です。
はじめに、バランシングスコア$b(x)$を考える。
バランシングスコアとは、条件付けることにより共変量と割り当てが独立にあるような共変量の関数である。
x \perp\!\!\!\!\perp z|b(x)
または、次のように関数$g$を利用して表現する。
p(z=1|x)=g(b(x))
バランシングスコアが存在することにより、次式が成り立ち割り当てを共変量で説明するのと同じように説明する。
p(z=1|x,b(x))=p(z=1|x))=p(z=1|b(x))
ここで、強く無視できる割り当て条件が成立すると仮定すると、次のようになる。
p(z=1|y_1,y_0,b(x))=p(z=1|x)
以上から、次式が成立し、バランシングスコアを条件付ければ、割り当て$z$と潜在的な結果変量$y_1,y_0$は独立になる。
p(z|y_1,_y0,b(x))=p(z|b(x))
または、
(y_1,y_0) \perp\!\!\!\!\perp z|b(x)
傾向スコアの定義
第$i$対象者の傾向スコアは次のように定義する。
e_i=p(z_i=1|x_i) (0<e_i<1)
実際には、各対象者の傾向スコアの心理がわからないため、ロジスティック回帰モデル等を利用して推定することになる。
p(z_i=1|x_i)=e_i=\frac{1}{1+exp\{-\alpha^tx_i\}}
ノンパラメトリックなモデルにより推定することもできるが、次元の呪いが存在することから一般的には傾向スコアの推定はパラメトリックに行われることが多い。
傾向スコアを用いた具体的な解析方法
傾向スコアを用いると因果効果が推定できる理由
因果効果は傾向スコアを用いると次のように表すことができる。
E(y_1)-E(y_0)=E(y_1-y_0)=E_e(E(y_1-y_0|e))
傾向スコアが同じ値の各軍ごとに$y_1-y_0$を計算し、傾向スコアの分布で期待値を取ることにより因果効果を求めることができる。
また、バランシングスコアを条件付けから、次のような関係が成り立つ。
E(y_1|e)=E(y_1|e,z=1),E(y_0|e)=E(y_0|e,z=0)
そして、因果効果は次にように表すことができ、「$z=0$の時の$y_1$」や「$z=1$の時の$y_0$」のような観測されていない値について考える必要がなくなる。
よって、観測されている従属変数かと傾向スコアの情報を用いれば因果効果を推定することができる。
###傾向スコアによる具体的な調整法
傾向スコアは二段推階推定法であり、次のようなステップをとる。
(1) 傾向スコアの推定
母数の推定ちを用いて、各対象ごとに$z=1$に割り当てられる予測確率=傾向スコアを算出する。
(2)推定された傾向スコアを用いた推定
具体的に傾向スコアを用いて調整する方法はマッチング、層別、分散分析、IPWが提案されている。
ⅰ.マッチング
２つの群で傾向スコアが等しい対象者をペアにして、その差の平均を持って因果効果を推定する。
実際には、等しい値のペアは存在しない場合の方が多いため、最近傍のペアを用いる場合が多い。
ⅱ.層別判別
傾向スコアをサイズが等しくなるようなサブクラスにわけ、各クラスごとに処理群と対照群の差を計算し、期待値を求めるこで因果効果を推定する。
ⅲ.分散分析
割り当て変数と傾向スコアを説明変数とした線形回帰分析を行う。
傾向スコア解析の利点
・結果変量と今日変量の回帰モデルを仮定する必要なく、モデルを誤設定する可能性が低い
・共変量を傾向スコアの１変数に集約するため、２群において局外要因や共変量の値に重なりがない場合でも使える
・ロジスティック回帰分析のモデルがある程度真のモデルと異なっていてもロバストな結果が得られる
###傾向スコア解析の問題点
・３群以上を比較する場合には、各２群ごとに傾向スコアを推定する必要がある
・因果効果の推定値の標準誤差を推定することができない
・マッチングと層別解析では周辺期待値$E(y_1),E(y_0)$の推定ができない
・マッチングでの、ペアを作成する基準が恣意的になる
・マッチングで、サンプル数が２群で偏りがありペアが作成できなかったり、一方の群のサンプルが無駄になる
・共分散分析のモデルでは、傾向スコアと目的変数が線形な関係にある必要がある
Rで傾向スコアを求める
傾向スコアを用いた因果効果の推定をRを用いて試してみました。
用いたデータはMatchingパッケージに含まれるデータセットlalondeを利用した。
このデータは1976年に実施された米国職業訓練プログラムを受講した群としていない群で、その後の年収を追跡したものです。
今回は1978年の年収への職業訓練プログラムを受講の効果を推定します
共変量は年齢、教育年数74年の年収、75年の年収等の10個があります。
はじめに、回帰分析を行い職業訓練プログラムの効果を推定してみます。
> # 回帰分析で推定した場合
> lm.model <- lm(lalonde$re78 ~ ., data = lalonde)
> lm.model <- step(lm.model)
> summary(lm.model)
Call:
lm(formula = lalonde$re78 ~ educ + black + re74 + treat, data = lalonde)
Residuals:
Min 1Q Median 3Q Max
-10006 -4331 -1751 3226 54249
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.074e+03 1.885e+03 1.100 0.27191
educ 4.044e+02 1.726e+02 2.343 0.01958 *
black -2.191e+03 8.272e+02 -2.648 0.00838 **
re74 1.010e-01 5.747e-02 1.757 0.07954 .
treat 1.727e+03 6.259e+02 2.760 0.00603 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 6490 on 440 degrees of freedom
Multiple R-squared: 0.05098, Adjusted R-squared: 0.04235
F-statistic: 5.909 on 4 and 440 DF, p-value: 0.0001223
介入によるre78への因果効果は、treatの回帰係数から1727である。
職業訓練を受けると1978年の年収が1727ドル上昇することが期待できということになります。
> library(Matching)
> data(lalonde)
> # 傾向スコアを求めるロジスティック回帰モデルを生成
> logit.model <- glm(treat~age+educ+black+hisp+married+nodegr+re74+re75+u74+u75,
+ data=lalonde, family=binomial(link = "logit"))
> ps.score <- logit.model$fitted
> summary(logit.model)
Call:
glm(formula = treat ~ age + educ + black + hisp + married + nodegr +
re74 + re75 + u74 + u75, family = binomial(link = "logit"),
data = lalonde)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.4884 -0.9934 -0.8708 1.2242 1.7403
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.622e+00 1.102e+00 1.472 0.1410
age 8.276e-03 1.452e-02 0.570 0.5687
educ -8.282e-02 7.230e-02 -1.145 0.2520
black -2.216e-01 3.684e-01 -0.601 0.5476
hisp -8.557e-01 5.128e-01 -1.669 0.0952 .
married 1.960e-01 2.784e-01 0.704 0.4813
nodegr -8.981e-01 3.146e-01 -2.855 0.0043 **
re74 -4.466e-05 3.010e-05 -1.483 0.1380
re75 2.924e-05 4.788e-05 0.611 0.5414
u74 -1.927e-01 3.765e-01 -0.512 0.6088
u75 -3.369e-01 3.213e-01 -1.048 0.2945
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 604.20 on 444 degrees of freedom
Residual deviance: 584.26 on 434 degrees of freedom
AIC: 606.26
Number of Fisher Scoring iterations: 4
> # 分散分析による調整法
> resalt.model <- lm(lalonde$re78 ~ lalonde$treat+ps.score)
> summary(resalt.model)
Call:
lm(formula = lalonde$re78 ~ lalonde$treat + ps.score)
Residuals:
Min 1Q Median 3Q Max
-6948 -4532 -1798 2813 54023
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3475.3 1289.5 2.695 0.0073 **
lalonde$treat 1674.5 647.4 2.587 0.0100 *
ps.score 2716.5 3078.0 0.883 0.3779
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 6581 on 442 degrees of freedom
Multiple R-squared: 0.01955, Adjusted R-squared: 0.01511
F-statistic: 4.407 on 2 and 442 DF, p-value: 0.01273
介入によるre78への因果効果は、lalonde$treatの回帰係数から1674.5である。
職業訓練を受けると1978年の年収が1674.5ドル上昇することが期待できということになります。
傾向スコアによる調整をしたことで、回帰分析の場合よりも効果は小さいことがわかります。
続いで、マッチングにより因果効果を推定してみます。
> # マッチング
> resalt.maching <- Match(Y=lalonde$re78, Tr=lalonde$treat, X=logit.model$fitted.values)
> summary(resalt.maching)
Estimate... 2138.6
AI SE...... 797.76
T-stat..... 2.6807
p.val...... 0.0073468
Original number of observations.............. 445
Original number of treated obs............... 185
Matched number of observations............... 185
Matched number of observations (unweighted). 322
因果効果の推定量は2138.6となり、分散分析よりも大きく効果を見積もった。
###IPW推定量
ルービンは傾向スコアによる問題点を解消した、傾向スコアによる重み付け推定法を提案している。
逆確率による重み付けを用いたIPW推定量を呼ばれるものであり、$y_1,y_0$となる周辺期待値はそれぞれ次のようになる。
\hat E(y_1) =\frac{\sum_{i=1}^N\frac{z_iy_i}{e_i}}{\sum_{i=1}^N\frac{z_i}{e_i}} ,
\hat E(y_0) =\frac{\sum_{i=1}^N\frac{(1-z_i)y_i}{e_i}}{\sum_{i=1}^N\frac{(1-z_i)}{e_i}}
この方法を用いると、因果効果だけではなく、各郡の周辺期待値を計算することができる。
RでIPW推定量を求める
先ほどと同じデータを用いて、IPW推定量を推定したいと思います。
> # IPW推定値
> y <- lalonde$re78
> z1 <- lalonde$treat
> (ipw1 <- sum((z1*y)/ps.score) / sum(z1/ps.score))
[1] 6212.993
> (ipw0 <- sum(((1-z1)*y)/(1-ps.score)) / sum((1-z1)/(1-ps.score)))
[1] 4571.409
> (score <- ipw1 - ipw0)
[1] 1641.584
傾向スコアによる推定量より、効果量は少しだけ小さくなった。
傾向スコア解析の拡張
一般化推定方程式におけるIPW推定量
無作為割り当てがされていない２群に対して、「共変量の影響が除去された場合の回帰係数の差」の検討などを可能にする方法として一般化推定方程式におけるIPW推定量が提案されている。
一般化推定方程式とは、一般化線形モデルを階層化したデータに対応し、平構造に関心がある場合に利用する手法である。
通常の一般化推定方程式は、次式を解くことで回帰関数$\mu(\omega_i,\beta)$の母数$\beta$の推定量を得る。
\sum_{i=1}^N S_i(\beta)= \sum_{i=1}^N \frac{\partial\mu(\omega_i,\beta)}{\partial\beta^t}V_i^{-1}(y_i-\mu(\omega_i,\beta))=0
ここで、$V_i$は作業共分散行列である,$A_i$を対角成分とする行列、$\phi$を過分散パラメータとすると次のように表すことができる。
V_i=\phi A_i^{\frac{1}{2}} R_i A_i^{\frac{1}{2}}
ここで$R_i$は作業関数行列と呼ばれ、一定の条件が合致すれば、推定の一致性は保証される。
従属変数$y$が欠測する可能性がある、この時欠測を表す$z$が$w$飲みに依存するなら次のような関係が成立する。
p(z_i|yi,w_i)=p(z_i|w_i)
この時、次の方程式の解は$\beta$の一致推定量である。
\sum_{i=1}^N z_i \frac{\partial\mu(\omega_i,\beta)}{\partial\beta^t}V_i^{-1}(y_i-\mu(\omega_i,\beta))=0
これは、$E(zS(\beta)|\omega)=E(z|\omega)E(S(\beta)|\omega)$の関係が成立するためある。
ロビンスらは回帰モデルの母数推定のための重み付け一般化推定方程式を提案している。
欠測が起こるかどうかについて独立変数と共変量で説明する次のモデルを考える。
p(z_i|x_i,w_i,\alpha)
この時の母数の最尤推定量$\hat \alpha$を用いて計算できる。
各対象者の観測確率の推定量の逆数
\chi_i(\hat \alpha)=\frac{1}{p(z_i=1|x_i,\omega_i,\hat \alpha)}
を重みとして用いた次の一般化推定方程式を解くことで得られる$\beta$の推定量が一致性をもつ。
\sum_{i=1}^N \chi(\hat \alpha) z_i \frac{\partial\mu(\omega_i,\beta)}{\partial\beta^t}V_i^{-1}(y_i-\mu(\omega_i,\beta))=0
一般化傾向スコア
2群の比較を行ってきたが、３群以上が存在する場合でも傾向スコアを計算することができる。
割り当て変数$z$が1からJまでの値をとる場合、次のように考える。
\begin{equation}
z_j =
\begin{cases}
1, & z=j \\
0, & \text{otherwise}
\end{cases}
(j=1,\cdots,J)
\end{equation}
この割り当て変数を用いて計算する傾向スコアを一般化傾向スコアと呼ばれる。
e_j=p(z_j=1|x) (j=1,\cdots ,J)
この時の前提条件は、強く無視できる割り当て条件から、「各群jごとに$z_j$とj番目の$y_j$が、共変量$x$を所与として独立」という、より弱い条件に緩和することができる。
今回はここまでで、3章の続きは次回で。