観察データでの効果推定（傾向スコア、IPW、DR）

まずは教科書を紹介します。

Rosenbaum先生は傾向スコア（propensity score）を提案した方です。

この教科書に書いていあるのは傾向スコアについてだけで、IPWやDRは書いてありません。

日本語はこちらの星野先生の本。

IPWやDRのことも書いてあり、とても充実した内容になっています（今回とても参考になりました）。

観察研究はランダム化研究と違い、治療（曝露）対象者を選ぶことができません。

よって様々な交絡要因が混入している可能性があり、場合によっては調整する必要があります。

しかし調整できるのは、データとして収集している変数のみで、未測定の交絡要因は調整できません（ランダム化試験は理論上未測定の交絡要因も調整されています）。

交絡要因の調整方法は、次のようなものがあります。

1. 傾向スコア
   • restriction（制限、調整）
   • stratification（層別）
   • matching（マッチング）
2. IPW（Inverse Probability Weighting）
3. DR推定量（ダブリーロバスト推定量、doubly robust）

Rのパッケージはあまり見当たりませんが、意外に簡単にRで推定できます。

データはMatchingパッケージに入っているlalondeを使います。

• 傾向スコア

治療を割り付けられる「傾向」を使って調整する方法です。

まずロジスティック回帰で、ある対象者が治療が割り付けられる確率を推定します。

library(Matching)
data(lalonde)

#------傾向スコアの計算
Logit <- glm(treat ~ age + educ + black + hisp + married + nodegr + re74 + re75 + u74 + u75, 
	family = binomial, data = lalonde)
PS    <- Logit$fitted

次のように傾向スコアを使った調整（restriction）を行います。

#---restriction（制限、調整）
Restriction <- lm(lalonde$re78 ~ lalonde$treat + PS)
summary(Restriction)

これだけですので、傾向スコアでの調整は難しくありません。

ただ、データによっては効果が逆転してしまったり、標準誤差が大きくなってしまったりします。

またPSでは群間差は推定できるのですが、各群の期待値は推定できません。

• IPW

IPWは傾向スコア（PS）を利用して、「その治療群に割り付けられる確率」の逆数で重み付け推定します。

#------IPW（「その治療群に割り付けられる確率」の逆数で重み付けする）
Weight   <- lalonde$treat/PS + (1 - lalonde$treat)/(1 - PS)
IPweight <- lm(lalonde$re78 ~ lalonde$treat, weights = Weight)
summary(IPweight)

Y    <- lalonde$re78
(IPW1 <- sum( (Y/PS)[lalonde$treat == 1])/sum(1/PS[lalonde$treat == 1]))
(IPW0 <- sum( (Y/(1-PS))[lalonde$treat == 0])/sum(1/(1-PS)[lalonde$treat == 0]))

• DR

ダブリーロバスト推定量はIPWをさらに拡張したもので、「2種のモデルのどちらかが間違っていてもバイアスのない推定ができる」推定方法です。

次の2モデルです。

1. 割り付けに関するロジスティックモデル
2. 結果変数に関する線形モデル

#------DR（ダブリーロバスト推定量）
n1      <- nrow(lalonde[lalonde$treat == 1, ])
n0      <- nrow(lalonde[lalonde$treat == 0, ])
Linear1 <- lm(re78 ~ age + educ + black + hisp + married + nodegr + re74 + re75 + u74 + u75, data = lalonde[lalonde$treat == 1, ])
Linear0 <- lm(re78 ~ age + educ + black + hisp + married + nodegr + re74 + re75 + u74 + u75, data = lalonde[lalonde$treat == 0, ])
Y1      <- lalonde$re78[lalonde$treat == 1]
Y0      <- lalonde$re78[lalonde$treat == 0]
PS1     <- PS[lalonde$treat == 1]
PS0     <- PS[lalonde$treat == 0]

(DR1 <- 1/n1 * sum(Y1 + ( (1 - PS1)/PS1) * (Y1 - Linear1$fitted)))
(DR0 <- 1/n0 * sum(Y0/(1 - PS0) + (1 - 1/(1 - PS0))*Linear0$fitted))

これらの推定量にどのような違いが出ているかをプロットします。

#------未調整の各群平均値
(Ttest <- t.test(lalonde$re78[lalonde$treat == 1], lalonde$re78[lalonde$treat == 0]))

#------差にエラーバーを付ける
#---Ttest
Ttestmean <- Ttest$estimate[1] - Ttest$estimate[2]
TtestCI <- Ttest$conf.int[1:2]

#---PS restriction
PsEst    <- summary(Restriction)$coefficients
Psmean   <- PsEst[2, 1]
Psse     <- PsEst[2, 2]
PsRestCI <- c(Psmean - 1.96*Psse, Psmean + 1.96*Psse)

#---IPW
IpwEst   <- summary(IPweight)$coefficients
Ipwmean  <- IpwEst[2, 1]
Ipwse    <- IpwEst[2, 2]
IpwCI    <- c(Ipwmean - 1.96*Ipwse, Ipwmean + 1.96*Ipwse)

#---DR
n      <- nrow(lalonde)
DRadj  <- 1/n^2 * (sum(sqrt((1 - PS1)/PS1)*(Linear1$fitted - mean(Y1))) + 
	sum(sqrt(PS0/(1 - PS0))*(Linear0$fitted - mean(Y0))))^2
DRse   <- sqrt(Ipwse^2 - DRadj)
DRmean <- DR1 - DR0
DRCI   <- c(DRmean - 1.96*DRse, DRmean + 1.96*DRse)


#---プロット
x    <- 1:4
mean <- c(Ttestmean, Psmean, Ipwmean, DRmean)
lcl  <- c(TtestCI[1], PsRestCI[1], IpwCI[1], DRCI[1])
ucl  <- c(TtestCI[2], PsRestCI[2], IpwCI[2], DRCI[2])

plot(x, mean, xaxt="n", ylim = c(0, 3500))
arrows(x, ucl, x, lcl, length=.05, angle=90, code=3) 
axis(side=1, at=1:4, labels=c("Ttest", "PS", "IPW", "DR"), cex=0.7)

このように、今回のデータでは調整してもしなくてもあまり結果は変わっていません（DRの計算は少し自信がありませんが、、）。

しかし、sparseなデータ（疎なデータ）では効果が逆転したり過剰に大きくなったりする事が多いです。

結果を良く吟味して利用すると良いかもしれません。

追記

多変量回帰と傾向スコアでの調整を比較した論文がありました。

http://www.sciencedirect.com/science?_ob=ArticleURL&_udi=B7P72-52W3V5B-2&_user=10&_coverDate=05%2F16%2F2011&_rdoc=1&_fmt=high&_orig=gateway&_origin=gateway&_sort=d&_docanchor=&view=c&_acct=C000050221&_version=1&_urlVersion=0&_userid=10&md5=5d4425d80a3e20b79c70783597dfd51f&searchtype=a

本当に傾向スコアって有用なのか？というタイトルで興味深いです。

追記2

傾向スコアの原論文はこれです。

http://faculty.smu.edu/millimet/classes/eco7377/papers/rosenbaum%20rubin%2083a.pdf