About connecting the dots.

data science related trivial things

調査観察データにおける因果推論(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:まぁ二重にロバストの方は,周辺標準誤差を出せてはいませんが...