はじめに
介入の割り当て確率が共変量に依存して異なるとき，統制群と介入群の応答変数の平均の差によって介入効果を推定するとバイアスが生じます． もし，共変量の値毎に割り当て確率を決定する関数が分かっていると，推定量のバイアスを取り除くことができます．
ただし，確率を決定する関数が分かっている点がポイントで，テキトーに関連がありそうな共変量でロジスティック回帰をすればいいわけではありません．
以下にテキトーなロジスティック回帰だと失敗する例を示します．
バイアス割り当て実験
分析に使うパッケージの読み込み
library("MatchIt") library("Matching") library("cobalt") library("WeightIt") library("broom")
共変量xによって介入（トリートメント）の割り当て確率が異なる状況を考えます．
xは-5から5までの範囲に値をとる一様分布とします．
潜在的結果y0は，この共変量と相関しています．
真の介入効果はtau=3とします．よって潜在的結果y1はy1=y0+3で定義できます．
n <- 1000;# sample size tau <- 3;# 真の介入効果 x <- runif(n,-5,5)# 共変量 y0 <-5+2*x+ rnorm(n,mean=0,sd=1)#潜在的結果yo, xの影響を受ける y1 <- y0+tau #潜在的結果y1
介入の割り当て確率をロジスティック関数で定義します．
p <- 1/(1+exp(-x))# make probabilities by logistic function d <- rbinom(n,size=1,p)# treatment with bias
pは割り当て確率です．共変量xの値が大きいほど割り当て確率が高くなっています．
dは介入を表すダミー変数です．
dを確率pのベルヌーイ分布で定義するため，xの値が大きい個体ほど介入を受けやすくなります．
以上のコードで割り当てにバイアスがある状況を表現できました．
観測データyを以下のように定義します*1．
y <- c() for(i in 1:n){# 観測値yの生成 if(d[i]==1){y[i] <- y1[i]}else{y[i] <- y0[i]}} data <- data.frame(y,d,p,x)
介入があった個体については潜在的結果y1を，介入がなかった個体についてはy0を応答変数として観測します．
各個体は，y0もしくはy1のどちらかしか観測できないことに注意します（因果推論の根本問題）
共変量xと割り当て確率pの関係を図で確認しておきましょう．
以下，割り当て確率のことを「傾向スコア」と呼びます*2．
xと傾向スコアpの関係
共変量xと介入dの関係も図で確認しておきましょう．
共変量xの値が小さいほど介入の割り当てがなく，値が高いほど介入が割り当てられています．
このようなランダムでない割り当ては，介入効果の推定にバイアスをもたらします．
計算して確かめてみましょう．
> summary(lm(data=data,y~d)) Call: lm(formula = y ~ d, data = data) Residuals: Min 1Q Median 3Q Max -12.3599 -2.7284 0.0375 2.8791 11.4583 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.8242 0.1699 4.852 1.42e-06 *** d 11.8107 0.2384 49.549 < 2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 3.768 on 998 degrees of freedom Multiple R-squared: 0.711, Adjusted R-squared: 0.7107 F-statistic: 2455 on 1 and 998 DF, p-value: < 2.2e-16
介入効果の推定値は11.8107でした．真の効果は3ですから，過大推定になっています．
回帰ではなく，条件付き期待値の差をとっても同じ結果を得ます．
mean(y[d==1])-mean(y[d==0])
傾向スコアマッチング
傾向スコアはすでに計算済みですから，それをつかってマッチングデータを作り，ATEを推定します．
> matched.data <- Match(Y = data$y, Tr=data$d, X =p,estimand="ATE") > summary(matched.data)#Matched number of observations............... 500 Estimate... 3.4254 AI SE...... 0.23556 T-stat..... 14.541 p.val...... < 2.22e-16 Original number of observations.............. 1000 Original number of treated obs............... 508 Matched number of observations............... 1000 Matched number of observations (unweighted). 1152
直感的に言えば，関数Matchは「傾向スコアが同じだけど介入の有無が異なる個体をペアにしたデータ」を作ります．
結果，バイアスが減り，真の効果に近いATE推定値を得ることができました．
念のためcobaltパッケージで共変量の調整結果を可視化します
data2 <- matchit(d~x,data=data,method = "nearest",replace = T) love.plot(data2, thresholds = 0.1)#SMD
傾向スコアによるマッチングで共変量のバランスがとれたことが分かります．
真の傾向スコアが手に入らない場合
上記の実験では，真の傾向スコアが分かっている条件下でマッチングデータを作りました．しかし現実の観察データを使って分析する場合，傾向スコアが事前に分かっていることなどありません．そもそも傾向スコアを決定する共変量の全てがデータとして入手できる状況など，極めて希でしょう．
そこで次に，必要な共変量が手に入らない場合にGLMで傾向スコアを推定してマッチングに利用した場合，どのくらい真の効果量の推定にバイアスが生じるのかを確かめます． まず分析用のデータの生成です．
n <- 3000;# sample size tau <- 3;# 真の介入効果 x1 <- runif(n,-5,5)# 共変量 x2 <- runif(n,-5,5)# 共変量 y0 <-x1+x2#潜在的結果y0, x1,x2の影響を受ける y1 <- y0+tau #潜在的結果y1 p <- 1/(1+exp(-x1-x2))# make probabilities by logistic function # pは真の傾向スコア．共変量x1とx2の値に依存する． d <- rbinom(n,size=1,p)# treatment with bias y <- c() for(i in 1:n){# 観測値yの生成 if(d[i]==1){y[i] <- y1[i]}else{y[i] <- y0[i]}}
真の傾向スコアは共変量x1とx2によって定まると仮定しています．
次に，必要な共変量が入手できない場合とできた場合を比較します．
out.glm <- glm(d~x1,family = "binomial")#傾向スコアを過小定式化モデルで推定 out.glm2 <- glm(d~x1+x2,family = "binomial")#フルモデルで推定 p.predict <- out.glm$fitted.values p.predict2 <- out.glm2$fitted.values
p.predictはx2が足りない場合の傾向スコアの推定値で，p.predict2はx1とx2が正しくGLMに含まれる場合の推定値です．
結果を比較してみます．
まず共変量が足りない場合
matched.data <- Match(Y = y, Tr=d, X =p.predict,estimand="ATE") summary(matched.data) Estimate... 7.2681 AI SE...... 0.10906 T-stat..... 66.645 p.val...... < 2.22e-16 Original number of observations.............. 3000 Original number of treated obs............... 1509 Matched number of observations............... 3000 Matched number of observations (unweighted). 7892
真の効果3に対して推定値は7.2681ですから過大推定しています（GLMのMcFadden R2 は 0.2209495でした）．
このとき共変量のバランスをチェックすると，共変量x1に関しては調整後にバランスはとれています．
dat <- data.frame(d,x1,x2,y) match.out <- matchit(data=dat,d~x1,method = "nearest",replace = T) love.plot(match.out, thresholds = 0.1)
つまり「投入した共変量のバランスがとれたとしても，効果量の推定には失敗する可能性がある」ので注意が必要です．
次に必要な共変量が入手できる場合のマッチング推定の結果です．
matched.data2 <- Match(Y = y, Tr=d, X =p.predict2,estimand="ATE") summary(matched.data2) Estimate... 3.4448 AI SE...... 0.1594 T-stat..... 21.611 p.val...... < 2.22e-16 Original number of observations.............. 3000 Original number of treated obs............... 1509 Matched number of observations............... 3000 Matched number of observations (unweighted). 4921
バイアスが補正されました（GLMのMcFadden R2 は 0.5866154でした）．
実際の観測データを用いた分析では，傾向スコアの正しい定義に必要な共変量は不明です． データに含まれる共変量のバランスがとれたとしても，処置への割り当て確率を計算するモデルに欠落変数があると，処置効果の推定に失敗します． つまり処置に影響しそうな共変量を分析者が（勝手に）考えて傾向スコアを予測しても，うまくいく保証はありません*3．
補足
必要な共変量が揃っている場合の，誤差項の分散の影響も確かめてみました． 複数回マッチングさせるために関数化します．
f1 <- function(s){ n <- 1000;# sample size tau <- 3;# 真の介入効果 x <- runif(n,-5,5)# 共変量 y0 <-x+ rnorm(n,mean=0,sd=3)#潜在的結果yo, xの影響を受ける y1 <- y0+tau #潜在的結果y1 p <- 1/(1+exp(-x-rnorm(n,mean=0,sd=s)))# make probabilities by logistic function # pは真の傾向スコア．共変量xの値が大きいほど割り当て確率が高い d <- rbinom(n,size=1,p)# treatment with bias y <- c() for(i in 1:n){# 観測値yの生成 if(d[i]==1){y[i] <- y1[i]}else{y[i] <- y0[i]}} out.glm <- glm(d~x,family = binomial) Lm <- out.glm$deviance/(-2) L0 <- out.glm$null.deviance/(-2) p.predict <- pout.glm$fitted.values #予測した傾向スコアを使いMatchでマッチさせる matched.data <- Match(Y = y, Tr=d, X =p.predict,estimand="ATE") c(matched.data$est[1], 1-Lm[1]/L0[1])# McFadden R^2 }
7行目p <- 1/(1+exp(-x-rnorm(n,mean=0,sd=s)))の部分で，共変量xから傾向スコアを計算する過程で誤差を導入します．誤差項の標準偏差sが大きくなると，共変量xによって予測できる傾向スコアの精度が低下します．
xと傾向スコアpの関係．誤差項の標準偏差s=3
トリートメントの割り当て確率は，この傾向スコアによって決まるので，割り当てdを共変量xに回帰すると傾向スコアの予測値を得ます．
最後にMatch関数に予測した傾向スコアを与えてマッチングデータを作成します．計算の手続きをまとめ関数f1の引数は，誤差の標準偏差sのみです．アウトプットは「マッチングデータから推定したATE」と「傾向スコアを予測したロジスティック回帰モデルのMcFadden 疑似決定係数R2」です．
計算例を確認してみましょう．
Mac <- c() Est <- c() for(i in 1:50){ Est <- append(Est,f1(s=3)[1]) Mac <- append(Mac,f1(s=3)[2]) } > mean(Mac) [1] 0.2395028 > mean(Est) [1] 2.998437
McFadden疑似決定係数が0.24程度でも，真の効果をうまく推定できています．
Mac <- c() Est <- c() for(i in 1:50){ Est <- append(Est,f1(s=10)[1]) Mac <- append(Mac,f1(s=10)[2]) } > mean(Mac) [1] 0.03893244 > mean(Est) [1] 3.002857
McFadden疑似決定係数が0.039程度でも，真の効果をうまく推定できています．この計算例から，疑似決定係数の大きさはバイアス補正にとって本質的ではないことが分かります．たとえ疑似決定係数が高くても，必要な共変量が足りなければ推定には失敗します．
逆確率による重み付け
傾向スコアを使ったIPWによる調整も試しておきます．ウェイトの計算にはWeightItパッケージを使いました．
weightit関数がdataを引数にとるのでdata.frameを定義しておきます．
データセットの定義
n <- 3000;# sample size tau <- 3;# 真の介入効果 x1 <- runif(n,-5,5)# 共変量 x2 <- runif(n,-5,5)# 共変量 y0 <-x1+x2#潜在的結果y0, x1,x2の影響を受ける y1 <- y0+tau #潜在的結果y1 p <- 1/(1+exp(-x1-x2))# make probabilities by logistic function d <- rbinom(n,size=1,p)# treatment with bias y <- c() for(i in 1:n){# 観測値yの生成 if(d[i]==1){y[i] <- y1[i]}else{y[i] <- y0[i]}} dat <- data.frame(d,x1,x2,y)
過小定式化モデル（共変量x2が足りないモデル）と適切なモデルによって傾向スコアを計算し，それぞれ重みに変換して回帰に使用します．
## 重み付きデータでの効果の推定 傾向スコアを過小定式化モデルで推定 weighting <- weightit(data=dat, d ~ x1,method = "ps",estimand = "ATE") IPW_result <- lm(y ~ d,weights = weighting$weights) %>% tidy() IPW_result ## 重み付きデータでの効果の推定 フルモデルで推定 weighting2 <- weightit(data=dat, d ~ x1+x2,method = "ps",estimand = "ATE") IPW_result2 <- lm(y ~ d,weights = weighting2$weights) %>% tidy() IPW_result2
結果の比較です．
> IPW_result # A tibble: 2 x 5 term estimate std.error statistic p.value <chr> <dbl> <dbl> <dbl> <dbl> 1 (Intercept) -2.16 0.0679 -31.9 4.56e-192 2 d 7.45 0.0959 77.7 0. > IPW_result2 # A tibble: 2 x 5 term estimate std.error statistic p.value <chr> <dbl> <dbl> <dbl> <dbl> 1 (Intercept) -0.269 0.0923 -2.92 3.58e- 3 2 d 4.02 0.134 30.1 6.63e-174
共変量が足りないモデルではうまくバイアスが補正できませんでした． 共変量がそろったモデルではバイアスが（相対的にはうまく）補正できました．ただし傾向スコアが0に近い値がデータセットに多く含まれるせいで，推定結果にはまだバイアスが残っています．
