調査観察データにおける因果推論(3) - Rによる傾向スコア,IPW推定量,二重にロバストな推定量の算出
目次
Rでの実行
さて,前回まででざっと理屈を紹介したところで,実際にRのコードをみていきましょう.使うデータは{Matching}パッケージのlalondeデータです.これは1976年に実施された米国職業訓練プログラムを受講した群としていない群で,訓練実施後の1978年の年収にどの程度違いが出たか,という分析のデータになります.共変量として,年齢,教育年数,黒人,ヒスパニック,既婚,高卒以上,74年の年収,75年の年収,74年に収入ゼロ,75年に収入ゼロ,の10個があります.
傾向スコアの算出
まずは傾向スコアを算出しましょう.共変量をすべて使って,ロジスティック回帰で処置群か否かを予測します.ぶっちゃけあまり効果はなさそうな感じなのですが,まぁそれはおいといて次にいきます.
logit.ps <- glm(treat~age+educ+black+hisp+married+nodegr+re74+re75+u74+u75, family=binomial, data=lalonde) ps <- logit.ps$fitted summary(logit.ps)
Call: glm(formula = treat ~ age + educ + black + hisp + married + nodegr + re74 + re75 + u74 + u75, family = binomial, 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
続いて傾向スコアを共変量として,処置群から1978年の年収を予測する回帰分析を行います.
ps.reg <- lm(lalonde$re78~lalonde$treat+ps) summary(ps.reg)
lalonde$treatの係数が1674.5なので,職業訓練を受けることで,1978年の収入が1674.5ドル上昇したことが示されています.しかしこのやり方だと,訓練を受けた群と受けなかった群それぞれの,年収の平均値を取得することはできません.
\Call: lm(formula = lalonde$re78 ~ lalonde$treat + ps) 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 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
IPW推定量の算出
IPW推定量とDR推定量については,周辺期待値まで含めて求められる,既存のパッケージは存在しません*1.なので,手計算する必要があります.そこで前回書いた,IPW推定量の周辺期待値の式を,そのままベタ書きしてみます.
y <- lalonde$re78 z1 <- lalonde$treat (ipwe1 <- sum((z1*y)/ps)/sum(z1/ps)) (ipwe0 <- sum(((1-z1)*y)/(1-ps))/sum((1-z1)/(1-ps)))
> (ipwe1 <- sum((z1*y)/ps)/sum(z1/ps)) [1] 6212.993 > (ipwe0 <- sum(((1-z1)*y)/(1-ps))/sum((1-z1)/(1-ps))) [1] 4571.409
ということで,ぶじに周辺期待値が求められました.共変量調整をおこなったあとの,職業訓練を受けた場合の期待年収は6212.993ドル,受けなかった場合は4571.409となります.両群の差の,6212.993-4571.409=1641.584ドルが因果効果ということになります.ただしこのやり方だと,各群の期待値の標準誤差を求めることができないので,別のやり方で算出することにします.手順をまとめて関数にしました.
ipwe <- function(data, target, treat, ps) { y <- as.matrix(data[target])[,] z1 <- as.matrix(data[treat])[,] z2 <- 1-z1 z <- cbind(z1, z2) ipw1 <- (z1/ps)/(length(z1)/sum(z1)) ipw2 <- (z2/(1-ps))/(length(z2)/sum(z2)) ipw <- ipw1+ipw2 ipw12.reg <- lm(y~z-1, data=data, weights=ipw) return(ipw12.reg) } ret <- ipwe(lalonde, "re78", "treat", ps)
今度は回帰式の形で結果が得られました.群毎の期待値に加え,標準誤差も得られています.期待値が上の結果と等しくなっていることを確認してください*2.ちなみに因果効果(群間の差分)の標準誤差は,sqrt(487.1^2+409.5^2)=636.3621になります.
Call: lm(formula = y ~ z - 1, weights = ipw) Weighted Residuals: Min 1Q Median 3Q Max -7167 -4455 -1727 3036 53962 Coefficients: Estimate Std. Error t value Pr(>|t|) zz1 6213.0 487.1 12.76 <2e-16 *** zz2 4571.4 409.5 11.16 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 6612 on 443 degrees of freedom Multiple R-squared: 0.3934, Adjusted R-squared: 0.3907 F-statistic: 143.7 on 2 and 443 DF, p-value: < 2.2e-16
DR推定量の算出
こちらも前回書いた,DR推定量の算出式をベタに書き起こすことにします.これも処理が長いので,関数化しました.DR推定量を求めるためには,共変量から目的変数を予測する回帰関数が必要なため,第5引数に条件式をそのまま引き渡しています*3.このあたりはiAnalysis倉橋さんのエントリを参考にしました*4.
dre <- function(data, target, treat, ps, formula) { n <- nrow(data) y <- data[target] data1 <- data[data[treat]==1,] data0 <- data[data[treat]==0,] model1 <- lm(formula=formula, data=data1) model0 <- lm(formula=formula, data=data0) fitted1 <- predict(model1, data) fitted0 <- predict(model0, data) dre1 <- (1/n)*sum(y+((z-ps)/ps)*(y-fitted1)) dre0 <- (1/n)*sum(((1-z)*y)/(1-ps)+(1-(1-z)/(1-ps))*fitted0) return(c(dre1, dre0)) } ret <- dre(lalonde, "re78", "treat", ps, re78~age+educ+black+hisp+married+nodegr+re74+re75+u74+u75)
得られた周辺期待値は以下の通りです*5.因果効果は6157.951-4586.173=1571.778となり,IPW推定量とそんなに変わらない値が得られましたね.
> ret [1] 6157.951 4586.173
傾向スコア関連のパッケージ
このエントリの下調べをしている最中に,いくつか傾向スコア関連のパッケージを見つけました.{CBPS}{CausalGAM}{iWeigReg}あたりですが,この中で一番使い勝手が良さそうなのは{CBPS}なので,これについて,上記のベタ書き計算式の検算もかねて確認してみたいと思います.IPW2*6とDREの点推定値が,それぞれ上記の自作関数の値と一致しているのがわかるかと思います.IPW2の標準誤差が大きく食い違っているのが気になりますが,正直ここのところはよくわかりません... 誰かご存知の方がいたら教えてほしいです.
ちなみにこの{CBPS}パッケージ,プリンストン大の今井先生らが提案している,Covariate Balancing Propensity Scoreという改良版の傾向スコアを算出する関数,CBPS()を含んでいるので,そのうち機会があれば,このスコアと通常の傾向スコアの比較も行ってみたいと思います.
ipwe.diff <- IPW(outcome=re78, treat=treat, data=lalonde, pscore=ps, k = length(coef(logit.ps))) dre.diff <- DR(re78~age+educ+black+hisp+married+nodegr+re74+re75+u74+u75, model="lm", data=lalonde, treat=treat, pscore=ps)
> ipwe.diff Point.Est Std.Err t.Statistic IPW1 1605.505 3053.2231 0.5258396 IPW2 1641.584 1818.6308 0.9026485 IPW3 1641.196 691.0218 2.3750277 > dre.diff Point Est Std. Err t-Statistic 1571.777258 669.239690 2.348601
ということで,傾向スコアを用いたIPW推定量,さらには二重にロバストな推定量の期待値と標準誤差,また周辺期待値と周辺標準誤差の算出を行いました*7.今日はこの辺で.
*1:因果効果とその標準誤差を推定できる{CBPS}パッケージはありますが...
*2:丸め誤差によって小数点2桁以下の部分は不正確ですが...
*3:ちなみに関数の中では,条件式を受け取ってlm()で推定を行っていますが,実際にはglm()を使ってより目的変数にフィットしたリンク関数を用いることも可能です
*4:倉橋さんのDR推定量の算出式に間違いがあったので,それを修正したら,modestな値に収まりました.また下で確認しているように,ほかのパッケージの算出結果とも一致したので,この算出式で正しいものと思います
*5:周辺標準誤差に関しては,算出式が星野本自体に書いていなかったので,よくわかりません.ちゃんと元論文を読めば書いてあるのかもしれませんが...
*6:IPWには式のパターンに応じて,IPW1, IPW2, IPW3の3種類があるそうです.詳しくはこちらの論文を参照してください.http://www.uni-konstanz.de/FuF/wiwi/workingpaperseries/WP_05-Pohlmeier-Seiberlich-Uysal_2013.pdf
*7:まぁ二重にロバストの方は,周辺標準誤差を出せてはいませんが...