StatModeling Memorandum

StatModeling Memorandum

StanとRとPythonでベイズ統計モデリングします. たまに書評.

傾向スコアを使ったベイズモデル

BUGS

傾向スコアについては以下を参考にして下さい。

  • [1]統計的因果推論（傾向スコア）の勉強会資料をアプしてみた（web）
  • [2]傾向スコアを用いた共変量調整による因果効果の推定と臨床医学・疫学・薬学・公衆衛生分野での応用について（pdf file）
  • [3]傾向スコア解析法による因果効果の推定と調査データの調整について（web）
  • [4]調査観察データにおける因果推論(3) - Rによる傾向スコア，IPW推定量，二重にロバストな推定量の算出（web）

概略は以下です。 例えば喫煙の体への影響（例えばガンになったかどうか）を知りたいとします。被験者のデータ（性別・年齢・学歴・ストレスがあるか・飲酒するかどうかなど）とガンかどうかのデータがあるとします。この場合、喫煙する人は男性で年齢が高めでストレスがあって飲酒をする人がとても多いため、通常の解析では喫煙自体がその人のガンの発症に影響があるか（因果効果）がうまく推定できない場合があります。このように無作為化試験が難しい/機能していない場合に傾向スコアを使った解析が有用な様子です。

具体的には、被験者の性質からロジスティック回帰分析によって喫煙者群か非喫煙者群のどちらに割り付けられるかの確率を推定します。これを傾向スコア（propensity score）と呼びます。実際には割り付けの確率が出力されるもの（ランダムフォレストとかSVMとか）なら何を使ってもよいと思います。次にこの傾向スコアを使って重みづけを行って解析します。解析の仕方は色々ある様子です。詳しくは上記のリンク先を見てください。

個人的には、傾向スコアを使った解析は、偏りや欠測のあるデータにおいて、どのデータのあてはまりをよくしたいかの重み付けを変える一つの枠組みであると感じました。

以下では傾向スコアを題材にする際によく使われるデータセットを使ってStanで解析します。Rの{Matching}パッケージのlalondeデータです（445人分データ）。これは1976年に実施された米国職業訓練プログラムを受講した群としていない群（treat=1 or 0）で、訓練実施後のre78（1978年の年収）にどの程度違いが出たかというデータになります。他の説明変数として、age（年齢）・educ（教育年数）・black（黒人か）・hisp（ヒスパニックか）・married（既婚か）・nodegr（高卒以上か）・re74（74年の年収）・re75（75年の年収）・u74（74年に収入ゼロか）・u75（75年に収入ゼロか） の10個あります。

まずは散布図行列を見てみましょう。

当然、re74・re75・u74・u75の間には相関がありますね。u74・u75を使うのは疑義があると思いますが今はつっこみません。またre78の分布は明らかに0が多く、Zero-Inflated XXX な分布です。このデータを正規分布を使った線形モデルで行うのは相当突っ込みどころがあると思うのですが…今回は説明を傾向スコアにしぼりたいのでそれもスルーします！詳しくはこちらの記事を参照してください。

このデータに対してRで単純に線形回帰（lm）を行うと以下の結果が得られます。

Call:
lm(formula = re78 ~ age + educ + black + hisp + married + nodegr +
    re74 + re75 + u74 + u75 + treat, data = lalonde)

Residuals:
   Min     1Q Median     3Q    Max
 -9612  -4355  -1572   3054  53119

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept)   256.6702  3521.6821    0.07   0.9419
age            53.5709    45.8063    1.17   0.2428
educ          400.7703   228.8217    1.75   0.0806
black       -2037.3328  1173.7163   -1.74   0.0833
hisp          425.8185  1564.5533    0.27   0.7856
married      -146.3292   882.2569   -0.17   0.8683
nodegr        -15.1794  1005.8851   -0.02   0.9880
re74            0.1234     0.0878    1.40   0.1608
re75            0.0197     0.1503    0.13   0.8955
u74          1380.2849  1187.9168    1.16   0.2459
u75         -1071.2155  1024.8779   -1.05   0.2965
treat        1670.7095   641.1323    2.61   0.0095

Residual standard error: 6520 on 433 degrees of freedom
Multiple R-squared:  0.0582,    Adjusted R-squared:  0.0343
F-statistic: 2.43 on 11 and 433 DF,  p-value: 0.00597

しかしながら、上記[2]によるとオススメできない、とのこと。該当部の引用は以下の通り。

共分散分析的手法の最大の問題点は，従属変数と共変量の関係をモデル化する必要があるということである．このモデル化では回帰関数を間違って指定すると（例えば二次関数の関係があるのに線形と仮定すると）誤った結果が導かれることなどがよく知られている．また，共分散分析的手法で推定することができる回帰係数自体は因果効果と等しくはならない．

というわけで傾向スコアを使ってみます。僕が思いついた範囲で通常の傾向スコアを使ったフローの不満点としては、傾向スコア自体にも確率的な変動があるはずなのに、そこは固定値として求めて最後の因果効果だけ確率的な変動幅を検討する点です。そこでBUGSやStanで実装し、ベイズ推定します。

注意点は「傾向スコアを求めるところ」と「傾向スコアを使って重み付けを行ってフィッティングするところ」を同時に推定してはいけない点です。これをやると回帰のフィッティングの影響が傾向スコアまでいってしまう“逆流”が起きてしまいます。そこで、WinBUGSの上級者向けのcut関数とさらに尤度を自分でいじる方法によりその問題を回避しました。

まずWinBUGSで尤度をいじる方法についてです。これはzeros trickと呼ばれます。BUGS Bookの9.5章「New distributions」に記載があります。

model {
   for (i in 1:N) {
      z[i]   <- 0
      z[i]    ~ dpois(phi[i])
      phi[i] <- -log(my.likelihood)
      my.likehood <- ...
   }
   ...
}

0 ~ dpois(phi[i])はシステムの持っている尤度の値に poisson(0|phi[i]) を乗算することを意味します（Stanでは increment_log_prob(log poisson(0|phi[i])) に相当します）。ここでpoisson(0|phi[i])は exp(-phi[i]) なので、phi[i]に -log(my.likehood) をいれておけばmy.likehoodが乗算されることになります。なんとad hocな。使うと可読性が落ちるのでどうしてもの時だけにしましょう。

これを使ったBUGSコードは以下になりました。cut関数を使って“逆流”を防いでいます（cut関数についてはThe BUGS Bookのp.201-203「9.4 Cutting feedback」の章を参照）

model {
   for(n in 1:N){
      Gr[n] ~ dbern(e[n])
      logit(e[n]) <- alpha.x + inprod(beta.x[], X[n,])
      w.star[n] <- equals(1, Gr[n]) * (1/e[n] - 1/(1-e[n])) + 1/(1-e[n])
      w[n] <- cut(w.star[n])

      mu[n] <- alpha.y + beta.y*Gr[n]
      minus.log.Lik[n] <- log(s.y) + 0.5*pow((Y[n] - mu[n])/s.y, 2)
      phi[n] <- w[n]*minus.log.Lik[n] + Const
      z[n] <- 0
      z[n] ~ dpois(phi[n])
   }

   alpha.x ~ dnorm(0, 1.0E-4)
   for(k in 1:K){
      beta.x[k] ~ dnorm(0, 1.0E-4)
   }
   alpha.y ~ dnorm(0, 1.0E-4)
   beta.y ~ dnorm(0, 1.0E-4)
   s.y ~ dunif(0, 1.0E+2)

   for(n in 1:N){
      w0[n] <- equals(0, Gr[n]) * w[n];
      w1[n] <- equals(1, Gr[n]) * w[n];
      w0.y[n] <- w0[n] * Y[n];
      w1.y[n] <- w1[n] * Y[n];
   }
   ipw0 <- sum(w0.y[])/sum(w0[]);
   ipw1 <- sum(w1.y[])/sum(w1[]);
   causal.effect <- ipw1 - ipw0;
}
  • 3行目：e[n]が傾向スコアになります。
  • 5行目: 傾向スコアを使った重み付けです。割り付けとしてはGr=1（職業訓練を受けた）だけど説明変数からはGr=1にはいなさそうな人のデータ（Y）ほど当てはまってほしいので、このように傾向スコアの逆数が重みになっています。Gr=0の場合も同様です。
  • 8～10行目: 各被験者のYをGrであてはめているのですが、それぞれの対数尤度の重みを傾向スコアで調整しています。この部分は上記の文献[2]のp.234の「（6）重み付け M推定量」に沿って実装しています。星野本ですとp.79の「3.4 一般的な周辺パラメトリックモデルの推定」の部分になります。尤度をいじるときには10行目のようにConstを加え、Constにある程度大きな値を渡すことでポアソン分布のパラメータが負になるのを防いでいます。書籍を参照。
  • 23-31行目: IPW推定量を計算して因果効果を推定しています。

キックするRコードは以下です。

library(R2WinBUGS)
library(Matching)

data(lalonde)
d <- lalonde
N <- nrow(d)
Y <- d$re78/10000
Gr <- d$treat
Expvar <- c('age', 'educ', 'black', 'hisp', 'married', 'nodegr', 're74', 're75', 'u74', 'u75')
X <- d[ , Expvar]
X$re74 <- X$re74/10000
X$re75 <- X$re75/10000
K <- ncol(X)
data <- list(
   N = N,
   K = K,
   Y = Y,
   X = as.matrix(X),
   Gr = Gr,
   Const = 1000
)

parameters <- c('alpha.x', 'beta.x', 'alpha.y', 'beta.y', 's.y', 'e', 'causal.effect')

inits <- function(){
   list(
      alpha.x = rnorm(1, 0, 0.1),
      beta.x = rnorm(K, 0, 0.1),
      alpha.y = rnorm(1, 0, 0.1),
      beta.y = rnorm(1, 0, 0.1),
      s.y = 1
   )
}

set.seed(1234)
post.bugs <- bugs(
   model.file = 'model.bugs',
   data = data, inits = inits, parameters.to.save = parameters,
   n.chains = 3, n.iter = 5000, n.burnin = 1000, n.thin = 5,
   debug = T, bugs.seed = 1234,
   working.directory = getwd(),
   bugs.directory = "C:/Program Files/WinBUGS14/"
)

post.mcmc <- as.mcmc(post.bugs$sims.matrix)
post.list <- mcmc.list(
   lapply(seq(length = post.bugs$n.chain),
      function(chain) as.mcmc(post.bugs$sims.array[, chain,])
   )
)

# print(post.bugs)
options(scipen = 10)
df <- data.frame(post.bugs$summary)
write.table(df, file = "post.summary.txt", sep = "\t", quote = F, col.names = NA)

以下は結果になります。計算時間は1chainあたり200秒ぐらいでした。スケーリングをもとに戻すにはbeta.yやcausal.effectは10000倍すればOKです。

meansdX2.5.X25.X50.X75.X97.5.Rhatn.eff
alpha.x1.718 1.078 -0.359 0.987 1.702 2.433 3.837 1.00 810
beta.x[1]0.008 0.015 -0.022 -0.002 0.008 0.019 0.036 1.01 400
beta.x[2]-0.085 0.069 -0.222 -0.134 -0.082 -0.038 0.048 1.01 650
beta.x[3]-0.240 0.371 -0.968 -0.497 -0.235 0.006 0.492 1.00 2400
beta.x[4]-0.911 0.524 -1.935 -1.260 -0.909 -0.548 0.124 1.00 2400
beta.x[5]0.203 0.284 -0.363 0.005 0.205 0.397 0.756 1.00 1900
beta.x[6]-0.932 0.310 -1.542 -1.142 -0.927 -0.720 -0.334 1.00 2400
beta.x[7]-0.478 0.313 -1.126 -0.686 -0.467 -0.263 0.112 1.00 440
beta.x[8]0.291 0.502 -0.723 -0.031 0.294 0.615 1.284 1.00 2400
beta.x[9]-0.183 0.391 -0.938 -0.438 -0.184 0.082 0.568 1.01 170
beta.x[10]-0.378 0.325 -1.017 -0.600 -0.374 -0.157 0.273 1.00 750
alpha.y0.459 0.035 0.389 0.435 0.460 0.483 0.526 1.00 470
beta.y0.166 0.050 0.071 0.133 0.165 0.201 0.267 1.00 750
s.y0.686 0.032 0.635 0.666 0.683 0.703 0.758 1.00 950
(eは略)………………………
causal.effect0.168 0.023 0.131 0.153 0.166 0.180 0.215 1.01 640
deviance890821.792 68.480 890700 890800 890800 890900 891000 1.00 1

追記