傾向スコアによる分析

HAD17.10から傾向スコアによる因果効果の推定ができるようになりました。詳細はこちらの記事をご覧ください。

この記事では、傾向スコアを使った分析で本当に共変量の影響を取り除いた因果効果の推定ができるのかを、Rによるシミュレーションで検証してみようと思います。

◆傾向スコアとは

理論的なところは他に良い資料があるのでざーっと飛ばします。

この記事とか、とても丁寧に書かれてて勉強になります。

すごいざっくり書くと、まず無作為割り当てがなされてない2つの群の平均値の差を見ても、因果効果の推定ができないのはご存知だと思います。

たとえばピアノを習うと頭が良くなるという仮説があるとします。しかし、ピアノを習った人と習ってない人の学力を比較しようと思っても、そもそもピアノを習うのは経済力がある家なので、親の学歴、教育への投資、などなどが交絡します。

そこで、親の学歴、教育への投資、収入などなどの共変量を測定して、それを統制すれば大丈夫と思いがちなんですが、そういうわけでもないのです。心理学で統制するといえば共分散分析、重回帰分析（その他線形モデル）がよく使われます。しかし、共変量から目的変数への実際の予測モデルが線形とは限りません。変数をうまく選び出せても、その効果がどういうモデルなのかまでわからないと、正しい統制ができないのです。

そこで傾向スコア分析が登場します。傾向スコアは、共変量で予測された割付の割合です。親の学歴が高いと子どもにピアノを習わせやすいとしたら、ピアノを習った群への割付が、高くなるということです。その確率を使ってうまく群ごとのバランスをとって、共変量の影響をなくしてしまおうというのが傾向スコアを使った因果推論の狙いです。

傾向スコアを使った分析が万能かといえば、もちろんそうではありません。傾向スコアによるバランシングがうまくいく条件として、「強く無視できる割り当て条件」というのがあります。ざっくりいえば、共変量で完全に割当が説明できるという条件です。なんか結構強い仮定じゃない？と思うかもしれませんが、共分散分析もこれと同様の仮定をおき、さらにモデルも線形である必要があります。前者だけの仮定で良いのが傾向スコアを用いた分析の利点です。

傾向スコアの推定は、ロジスティック回帰分析がよく使われます。共変量を説明変数、群分け変数を目的変数としたロジスティック回帰を行い、その予測値が傾向スコアになります。傾向スコアの推定が正しいためには、モデルをうまく指定しないといけないのですが、共変量の2次の項や、交互作用など、いろんなものをぶっこんでやればよいです。変数が増えると推定が不安定になりそうですが、共分散分析のときと違って、傾向スコアという1つの値に集約するので、そういう問題はないようです（ただ、機械学習の方法を使うとうまくいったりいかなかったりするようです。このあたりまだよくわかってません）。

◆逆確率重み付け推定（IPW推定）

傾向スコアを使った因果推論の方法にはいろいろあります。傾向スコアが同じぐらいの人でマッチングをするとか、傾向スコアを共変量にした回帰分析を使うなどです。しかし、それらの方法も、サンプルサイズが小さくなったり、標準誤差の推定がうまくできなかったり、推定量としての性能がよくなかったりします。

そこで逆確率重み付け推定法（Inverse Probability Weighting; IPW）が推奨されています。これは、群ごとに傾向スコアの逆数で重みづけて、その平均値の差を計算する方法です。この推定量は、平均因果効果の「強く無視できる割り当て条件」と「傾向スコアの推定が正しい」が満たされれば、一致推定量になることがわかっています。

ここまで読んできて、「いままで共分散分析使ってきたけど、具体的にどんな問題があるの？本当にその傾向スコアで重みづけたらうまくいくの？」という疑問が出てきた人もいるでしょう。なんぜ、僕もそう思ったうちの一人です。

そこで、以下では、IPW推定量がほんとうに一致性を持つのか、共分散分析とどう違うのかをシミュレーションしてみます。

◆Rでシミュレーション

以下のコードでシミュレーションしてみました。

set.seed(57)
ATE <- 0.5 #平均因果効果の真値
n <- 500
beta0 <- c(1,0.3, 0.3, 0.3) #群分け←共変量の回帰係数
beta1 <- c(1,0.5, 1.7, 0.6) #目的変数(Z=1の群)←共変量の回帰係数
beta2 <- beta1 + c(0, rnorm(1,0,0.5),rnorm(1,0,0.5),rnorm(1,0,0.5)) #Z=0の群←共変量の回帰係数
p <- length(beta0)
sigma <- 1.5

start.time<-proc.time()
dif <- c() #ただの差
ipw <- c() #逆確率重み付け推定量
reg <- c() #共分散分析
for(m in 1:1000){
  #説明変数作成
  X <- as.matrix(array(dim=c(n,p)))
  X[,1] <- rep(1,n)
  for(i in 2:p){
    X[,i] <- rnorm(n,0,1)
  }
  eta <- 1/(1+exp(-X%*%beta0))

  #群分け変数作成
  Z <- c()
  for(i in 1:n){
    Z[i] <- rbinom(1,size=1,eta[i,1])
  }

  #目的変数作成
  Y <- Z*ATE + Z*X%*%beta1 + (1-Z)*X%*%beta2 + rnorm(n,0,sigma)
  
  #傾向スコア
  result.logi <- glm(Z ~ X[,2] + X[,3] + X[,4],family="binomial")
  ps <- result.logi$fitted.value
  
  #重み
  weight <- Z/ps + (1-Z)/(1-ps)
  
  #ただの差
  dif[m] <- mean(Y[Z==1]) - mean(Y[Z==0]) 
  
  #逆確率重み付け推定量
  result.ipw <- lm(Y ~ Z,weights = weight)
  ipw[m] <- result.ipw$coefficients[2]
  
  #共分散分析
  result.reg <- lm(Y ~ X[,2] + X[,3] + X[,4] + Z)
  reg[m] <- result.reg$coefficients[5]
}
end.time<-proc.time()
(end.time-start.time)

mean(dif) #ただの差
mean(ipw) #逆確率重み付け推定量
mean(reg) #共分散分析

set.seed(57)

ATE <- 0.5 #平均因果効果の真値

n <- 500

beta0 <- c(1,0.3, 0.3, 0.3) #群分け←共変量の回帰係数

beta1 <- c(1,0.5, 1.7, 0.6) #目的変数(Z=1の群)←共変量の回帰係数

beta2 <- beta1 + c(0, rnorm(1,0,0.5),rnorm(1,0,0.5),rnorm(1,0,0.5)) #Z=0の群←共変量の回帰係数

p <- length(beta0)

sigma <- 1.5

start.time<-proc.time()

dif <- c() #ただの差

ipw <- c() #逆確率重み付け推定量

reg <- c() #共分散分析

for(m in 1:1000){

#説明変数作成

X <- as.matrix(array(dim=c(n,p)))

X[,1] <- rep(1,n)

for(i in 2:p){

X[,i] <- rnorm(n,0,1)

}

eta <- 1/(1+exp(-X%*%beta0))

#群分け変数作成

Z <- c()

for(i in 1:n){

Z[i] <- rbinom(1,size=1,eta[i,1])

}

#目的変数作成

Y <- Z*ATE + Z*X%*%beta1 + (1-Z)*X%*%beta2 + rnorm(n,0,sigma)

#傾向スコア

result.logi <- glm(Z ~ X[,2] + X[,3] + X[,4],family="binomial")

ps <- result.logi$fitted.value

#重み

weight <- Z/ps + (1-Z)/(1-ps)

#ただの差

dif[m] <- mean(Y[Z==1]) - mean(Y[Z==0])

#逆確率重み付け推定量

result.ipw <- lm(Y ~ Z,weights = weight)

ipw[m] <- result.ipw$coefficients[2]

#共分散分析

result.reg <- lm(Y ~ X[,2] + X[,3] + X[,4] + Z)

reg[m] <- result.reg$coefficients[5]

}

end.time<-proc.time()

(end.time-start.time)

mean(dif) #ただの差

mean(ipw) #逆確率重み付け推定量

mean(reg) #共分散分析

[crayon-673f0240a26d9717996751/]

コードの説明は・・・上のをみれば分かる人は分かるでしょう（おい）。

何をやっているかを簡単に説明すると、まず平均因果効果の真値を決めます。上の例では0.5にしています。つぎに、共変量X（3変数）を正規分布から生成して、そこからロジスティック回帰で群分け変数Zを生成します。ここでは、傾向スコアの推定はあっているものとします。また共変量と群分け変数を使って目的変数Yを生成しています。このとき、群ごとで共変量の効果が違っている（交互作用がある）ように生成しています。

あとは、傾向スコアをglm()で推定して、逆確率で重みづけた平均値の差をlm()で計算しています。同様に、何も統制しないただの差と、共分散分析（交互作用は入っていない）で推定された差を計算します。

このシミュレーションを1000回行いました。一致推定量なので本来はサンプルサイズを大きくすればいいのですが、不偏性もあるかを見るために推定量のその平均値を計算してみました。

結果はこちら

>mean(dif) #ただの差
[1] 1.113364
> mean(ipw) #逆確率重み付け推定量
[1] 0.5121877
> mean(reg) #共分散分析
[1] 0.3939989

ただの差は大きくバイアスがあるのがわかります。IPWは真値の0.5に近いです。共分散分析の結果は、過小推定されていますが、それは今回のモデルの設定に依存しているだけです。

仮に、共分散分析も正しいモデルの場合はどうなるかというと、

>mean(dif) #ただの差
[1] 1.294028
> mean(ipw) #逆確率重み付け推定量
[1] 0.4961258
> mean(reg) #共分散分析
[1] 0.4967441

というように、ちゃんと推定できます。

次に、傾向スコアのモデルが間違えているが、アホみたいに交互作用を入れまくった場合（正しいモデルは1次線形だが、全部の交互作用をいれて傾向スコアを推定する）にどうかるかといえば、

>mean(dif) #ただの差
[1] 1.294028
> mean(ipw) #逆確率重み付け推定量
[1] 0.5102679
> mean(reg) #共分散分析
[1] 0.4967441

このように、バイアスはあまりありません。

これらの結果から、傾向スコアの推定自体を複雑にしても、あまり問題になりません（最終的に重みとして使う変数は1つなので自由度が減ったりすることもない）。

今回の記事では紹介しませんが、傾向スコアのモデルが間違えていても、目的変数Yに対するモデルがあっていれば正しい推定ができる、「二重にロバストな推定量」というのもあります。こちらのほうがIPW推定量よりも仮定が少ない（傾向スコアか目的変数どちらかのモデルが正しいければいい）上に、標準誤差が小さくなるようです。

このように、傾向スコアを使った重み付けで、平均因果効果の平均的には良い推定量になっているのがわかりました。

めでたしめでたし。

上の記事は

『星野崇宏 (2009).　調査観察データの統計科学　因果推論・選択バイアス・データ融合　岩波書店』

を参考に書いています。