10月9日（木）
傾向スコア分析
傾向スコアについてＲでやっているサイトがあったので、それをほぼもろパクリして自分でやってみた。
install.packages("Matching") library(Matching) data(lalonde) logi <- glm(treat~age+educ+black+hisp+married+nodegr+re74+re75+u74+u75,family=binomial,data=lalonde) Y78 <- lalonde$re78 Tre <- lalonde$treat lm(Y78~Tre) mout1 <- Match(Y=Y78,Tr=Tre,X=logi$fitted) summary(mout1) > lm(Y78~Tre) Call: lm(formula = Y78 ~ Tre) Coefficients: (Intercept) Tre 4555 1794 > summary(mout1) 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
ちなみに
mout1 <- Match(Y=Y78,Tr=Tre,X=logi$fitted)
には感激した。
これで層化した分析をやってくれているようである。（どうやっているのか厳密なことはよく分からないが…）
これによって、「だいたい同じ条件を持っている人」同士で比較ができるので、バイアスの影響を取り除くことができる、ということらしい。
ただし、「アウトカムに影響しないが曝露率には影響する因子」「アウトカムに影響するが曝露率には影響しない因子」が含まれていたり、あるいは重要な交絡因子が含まれていなかったりすると不都合が起こるはずである。
ちょっと直観的にわかりにくい部分があるので、Ｒでシミュレーションして確かめてみる。
まずは理想的な場合においてのデータ作成＆効果の測定
kouraku1,kouraku2によって曝露率e、outcomeが影響を受けているモデル
effect <- 2 kouraku1 <- runif(100,0,1) kouraku2 <- runif(100,0,1) x <- kouraku1+kouraku2 e <- 1/(exp(-x)+1) bakuro <- c() for (i in 1:100) bakuro <- c(bakuro,sample(c(0,1),1,prob=c(1-e[i],e[i]))) outcome <- kouraku1+kouraku2+effect*bakuro data <- data.frame(cbind(kouraku1,kouraku2,bakuro,outcome)) logi <- glm(bakuro~kouraku1+kouraku2,family=binomial,data=data) mout1 <- Match(Y=outcome,Tr=bakuro,X=logi$fitted) summary(mout1) > summary(mout1) Estimate... 2.0176 AI SE...... 0.018846 T-stat..... 107.06 p.val...... < 2.22e-16 Original number of observations.............. 100 Original number of treated obs............... 77 Matched number of observations............... 77 Matched number of observations (unweighted). 77
これは申し分ない結果。
感動的である。ちなみにこのモデルで交絡を無視すると
> lm(outcome~bakuro) Call: lm(formula = outcome ~ bakuro) Coefficients: (Intercept) bakuro 0.6888 2.3023
こんなかんじ。
それぞれ何度も繰り返して箱ひげ図を書いたものが下の通り。
ここで、「アウトカムに影響しないが曝露率には影響する因子」あり（Ａ）、「アウトカムに影響するが曝露率には影響しない因子」あり（Ｂ）、交絡因子が抜けている（Ｃ）のそれぞれの場合についても、同じようにつくって箱ひげ図を描いて並べてみた。
effect <- 2 #理想的 hako1 <- function(){ kouraku1 <- runif(100,0,1) kouraku2 <- runif(100,0,1) x <- kouraku1+kouraku2 e <- 1/(exp(-x)+1) bakuro <- c() for (i in 1:100) bakuro <- c(bakuro,sample(c(0,1),1,prob=c(1-e[i],e[i]))) outcome <- kouraku1+kouraku2+effect*bakuro data <- data.frame(cbind(kouraku1,kouraku2,bakuro,outcome)) logi <- glm(bakuro~kouraku1+kouraku2+insi1,family=binomial,data=data) mout1 <- Match(Y=outcome,Tr=bakuro,X=logi$fitted) return(mout1$est) } #Ａ hako3 <- function(){ kouraku1 <- runif(100,0,1) kouraku2 <- runif(100,0,1) insi1 <- runif(100,0,1) x <- kouraku1+kouraku2+insi1 e <- 1/(exp(-x)+1) bakuro <- c() for (i in 1:100) bakuro <- c(bakuro,sample(c(0,1),1,prob=c(1-e[i],e[i]))) outcome <- kouraku1+kouraku2+effect*bakuro data <- data.frame(cbind(kouraku1,kouraku2,insi1,bakuro,outcome)) logi <- glm(bakuro~kouraku1+kouraku2+insi1,family=binomial,data=data) mout1 <- Match(Y=outcome,Tr=bakuro,X=logi$fitted) return(mout1$est) } #Ｂ hako4 <- function(){ kouraku1 <- runif(100,0,1) kouraku2 <- runif(100,0,1) insi2 <- runif(100,0,1) x <- kouraku1+kouraku2 e <- 1/(exp(-x)+1) bakuro <- c() for (i in 1:100) bakuro <- c(bakuro,sample(c(0,1),1,prob=c(1-e[i],e[i]))) outcome <- kouraku1+kouraku2+insi2+effect*bakuro data <- data.frame(cbind(kouraku1,kouraku2,insi2,bakuro,outcome)) logi <- glm(bakuro~kouraku1+kouraku2+insi2,family=binomial,data=data) mout1 <- Match(Y=outcome,Tr=bakuro,X=logi$fitted) return(mout1$est) } #Ｃ hako5 <- function(){ kouraku1 <- runif(100,0,1) kouraku2 <- runif(100,0,1) x <- kouraku1+kouraku2 e <- 1/(exp(-x)+1) bakuro <- c() for (i in 1:100) bakuro <- c(bakuro,sample(c(0,1),1,prob=c(1-e[i],e[i]))) outcome <- kouraku1+kouraku2+insi2+effect*bakuro data <- data.frame(cbind(kouraku1,kouraku2,bakuro,outcome)) logi <- glm(bakuro~kouraku1,family=binomial,data=data) mout1 <- Match(Y=outcome,Tr=bakuro,X=logi$fitted) return(mout1$est) } x <- c() for (i in 1:1000) x <- c(x,hako1()) A <- c() for (i in 1:1000) A <- c(A,hako3()) B <- c() for (i in 1:1000) B <- c(B,hako4()) C <- c() for (i in 1:1000) C <- c(C,hako5()) boxplot(x,A,B,C,names=c("理想的","Ａ","Ｂ","Ｃ"))
やはり、Ａ，Ｂはいずれも理想的な場合よりも分散の値が大きくなっているし、Ｃには交絡が残っている。
なお、確認として、Ａのデータについてinsi1を除いて傾向スコアを計算してやると、（結果をＡ.とする）
hako3. <- function(){ kouraku1 <- runif(100,0,1) kouraku2 <- runif(100,0,1) insi1 <- runif(100,0,1) x <- kouraku1+kouraku2+insi1 e <- 1/(exp(-x)+1) bakuro <- c() for (i in 1:100) bakuro <- c(bakuro,sample(c(0,1),1,prob=c(1-e[i],e[i]))) outcome <- kouraku1+kouraku2+effect*bakuro data <- data.frame(cbind(kouraku1,kouraku2,insi1,bakuro,outcome)) logi <- glm(bakuro~kouraku1+kouraku2,family=binomial,data=data) mout1 <- Match(Y=outcome,Tr=bakuro,X=logi$fitted) return(mout1$est) } A. <- c() for (i in 1:1000) A. <- c(A.,hako3.()) boxplot(x,A,A.,names=c("理想的","Ａ","Ａ."))
このように、理想的なときと同じ結果が得られる。
まとめると、傾向スコアの算出に用いる因子には、
「曝露率とアウトカムの両方と関係している交絡因子だけを、できるだけ多く加える」ことが重要であることがシミュレーションからも明らかとなった。
逆に言えば、その部分さえ十分正確に出来れば、そこそこ安心して結果を受け入れることができそうに思える。
（回帰モデルは、入れるべき因子の項はそれでいいのか？とか、そのモデル式で本当に正しいのか？とか、気になる点がたくさんあるが…。）
もちろん厳密な値ではないかもしれないけど、かといって大きく外れてさえなければそれでいいかな、というとき（実際にそういうことはよくあるように思われる）に、このようなセミパラメトリックな方法は相当な威力を発揮すると感じた。
以上
