GBMを用いた傾向スコア法

R

前回のまとめ記事からtwangパッケージを選んで実際の手順を回してみる。
twangを選んだ理由は、3群以上の群に対する適用が可能だったこと＋モデルがロジスティック回帰に比べて性能の良いGBMだったことの2点。
なお、twangの使い方についてはビネットにまとまっている。
http://cran.r-project.org/web/packages/twang/vignettes/twang.pdf
以下ビネットをかいつまみながら、下記手順に従って簡単にコードを追う。

  1. 傾向スコアの推定及びマッチング
  2. マッチング前後のバランスの確認
  3. 効果の推定

前準備

データについてはlalondeデータを用いる。
説明についてはヘルプを参照のこと。

install.packages("twang")
library(twang)
set.seed(1)
data(lalonde)

傾向スコアの推定及びマッチング

twangパッケージは傾向スコアの推定にGBMを用いる。
推定の際はps関数を使う。
この際パラメータをたくさん指定しているが、shurinkageまでがGBMのモデルのパラメータ、perm.test.iters以降がps独自のパラメータである。
stop.methodで収束条件を指定できる。
収束条件は介入群と非介入群における共変量の分布に差がないことをなにで測るか（標準化バイアスもしくはカルバックライブラー統計量）、
共変量全体の要約統計量を何とするか（平均もしくは最大値）の組み合わせで設定できる。
デフォルトで指定されているのは標準化バイアスの平均値（es.mean）、カルバックライブラー統計量(ks.max)
結果はsummary関数で確認する。
また収束の状況についてはplotでplotsパラメータを1に設定することで図示して確認できる。

# ps estimation and matching
ps.lalonde <- ps(treat~age+educ+black+hispan+nodegree+married+re74+re75,
                 data=lalonde,
                 n.trees=5000,
                 interaction.depth=2,
                 shrinkage=0.01,
                 perm.test.iters=0, #0より大きい数値にするとモンテカルロ法を利用する
                 stop.method=c("es.mean","ks.max"),
                 estimand="ATT",
                 verbose=FALSE
                 )
summary(ps.lalonde)
plot(ps.lalonde, plots=1) # 収束の状況をチェック

マッチング前後のバランスの確認

バランスはbal.tableで一望できる。
plotを用いて箱ひげ図でも確認できる。

(lalonde.balance <- bal.table(ps.lalonde))
plot(ps.lalonde, plots=2) # 2はboxplot

効果の推定

上記で各ケースの重みが求まるのでそれを元データ（lalonde）に付加した後、その重みを用いて効果を推定する。
重みを用いるにあたってsurveyパッケージを使う。

library(survey)
lalonde$w <- get.weights(ps.lalonde, stop.method="es.mean")
design.ps <- svydesign(ids=~1, weights=~w, data=lalonde)
res.glm1 <- svyglm(re78~treat,
               design=design.ps
               )
summary(res.glm1)

res.glm2 <- svyglm(re78~treat+age+educ+black+hispan+nodegree+married+re74+re75,
               design=design.ps
               )
summary(res.glm2)

参考資料

各統計量の意味などについては前回紹介した勉強会資料が詳しい。
http://blue.zero.jp/yokumura/Rhtml/session08.html
またsurveyパッケージの使い方については下記書籍が詳しい。
http://www.amazon.co.jp/dp/425412791X
全体を通したコードは下記の通り。

install.packages("twang")
library(twang)
set.seed(1)
data(lalonde)
# 傾向スコアの推定及びマッチング
ps.lalonde <- ps(treat~age+educ+black+hispan+nodegree+married+re74+re75,
                 data=lalonde,
                 n.trees=5000,
                 interaction.depth=2,
                 shrinkage=0.01,
                 perm.test.iters=0, #0より大きい数値にするとモンテカルロ法を利用する
                 stop.method=c("es.mean","ks.max"),
                 estimand="ATT",
                 verbose=FALSE
                 )
plot(ps.lalonde,plots=1) # 収束の状況をチェック

# マッチング前後のバランスの確認
(lalonde.balance <- bal.table(ps.lalonde))
summary(ps.lalonde)
plot(ps.lalonde, plots=2) # 2はboxplot

# 効果の推定
library(survey)
lalonde$w <- get.weights(ps.lalonde, stop.method="es.mean")
design.ps <- svydesign(ids=~1, weights=~w, data=lalonde)
res.glm1 <- svyglm(re78~treat,
                   design=design.ps
)
summary(res.glm1)
res.glm2 <- svyglm(re78~treat+age+educ+black+hispan+nodegree+married+re74+re75,
               design=design.ps
               )
summary(res.glm2)