« 読了：「往生要集 日本浄土教の夜明け」「浄土真宗とは何か 親鸞の教えとその系譜」「日本精神史 自然宗教の逆襲」 | メイン | 読了：Austin, Jembere, Chiu (2016) 層別クラスタ抽出標本の2群を傾向スコアでマッチングするとき傾向スコアの算出に標本ウェイトを使うべきかどうか調べてみたけどよくわかんなかった »
2017年8月22日 (火)
McCaffrey, D.F., Griffin, B.A., Almirall, D., Slaughter, M.E., Ramchand, R., Burgette, L.F. (2013) A tutorial on propensity score estimation for multiple treatments using generalized boosted models. Statistics in Medicine, 32, 3388-3414.
題名の通り、処理が3水準以上あるときに、generalized boosted modelを用いて傾向スコア調整するやり方についての長ーいチュートリアル。著者らはRのtwangパッケージの中の人。急遽実戦投入を迫られ、事前の儀式としてめくった。
処理の水準数を$M$ とする。ある人が処理 $t$ を受けた時のpotential outcomeを$Y[t]$ とする。ペアワイズの効果を$D[t', t''] = Y[t'] - Y[t'']$と書く。
因果的効果についての統計量として次の2つを考える。
- ATE(平均処理効果)。母集団全体が処理$t'$を受けた時の平均と、母集団全体が処理$t''$を受けた時の平均の差。すなわち$E(D[t', t'']) = E(Y[t']) - E(Y[t''])$。
- ATT(処理群の平均処理効果)。$t''$と比べた$t'$のATTとは、この研究において$t'$を受けた対象者さんたちの平均と、もしその人たちが$t''$を受けていたらどうなっていたかの平均との差、すなわち$E(Y[t'] | T = t') - E(Y[t''] | T = t')$。
$M=3$の場合、ATEは3つ、ATTは6つあることになる。
ATEとATTのちがいは、対象者間の効果の異質性から生まれる。[←あっ、そうか。効果がhomogeneousだったらどちらでも同じことだわな。なるほど...]
どういうときにどっちのestimandが適切か？ すべての処理が潜在的には母集団全員に適用可能なのであれば、ATEが自然。いっぽう、処理$t'$が現在のターゲットに対して適切なものかどうかに関心があるのならATTが自然。
いよいよ本題。どうやって推定するか。
個人$i$について、観察された処理を$T_i$、観察されたアウトカムを$Y_i$、共変量のベクトルを$\mathbf{X}_i$とする。
ここではIPTW (inverse probability of treatment weighting) 推定量について考えよう。この推定量は2つの想定を置く。どちらもデータからは検証できない想定である。
- 十分なオーバーラップ(positivityともいう)。すべての$\mathbf{X}$と$t$について$0 < pr(T_i = t | \mathbf{X}) < 1$。つまり、誰であれ、どの処理であれ、それを受ける確率は0じゃない、という想定である。
- 未知・未測定の交絡因子はないという想定(交換可能性ともいう)。$T_i[t] = I(T_i = t)$として、$T[t]$と$Y[t]$は$\mathbf{X}$の下で条件付き独立。
さて、$p_t (\mathbf{X}) = pr(T[t] = 1 | \mathbf{X})$を傾向スコアと呼ぶ。
上の2つの想定の下で、
$\displaystyle \hat{\mu}_t = \frac{\sum_i T_i[t] t_i w_i[t]}{\sum_i T_i[t] w_i[t]}$
は$E(Y[t])$の一致推定量になる。ただし$w_i[t]$とは傾向スコアの逆数、すなわち$w_i[t] = 1 / p_t(\mathbf{X}_i)$ね。ここからペアワイズATEが推定できる。
いっぽうペアワイズATTは... [めんどくさいので略]
ここからは、傾向スコアの推定方法。
もっとも一般的なのは多項ロジスティック回帰を使う方法である。しかし、共変量の交互作用項とかをどこまで入れるかの判断が難しくて... [いろいろ書いてあるけどパス]
そこでGBMを使おう。いろいろ比べたらGBMが一番良かったという話もあるぞ[McCaffrey, Ridgeway, & Morral (2004 Psych.Methods), Harder, Stuart, & Anthony (2010 Psych.Methods)というのが挙げられている。なお、この論文中には、GBMとはなんぞやという説明はほとんど出てこない。割り切っておるなあ]。
まずは処理が2水準の場合について。
GBMの反復をどこでストップするか。いくつかの基準がある。
- 標準化バイアス(絶対標準化平均差)。ある共変量について、処理群の加重平均と統制群の加重平均の差の絶対値を、重み付けしないSDで割った値。$k$番目の共変量について、ATEなら
$SB_k = | \bar{X}_{k1} - \bar{X}_{k0} | / \hat{\sigma}_k$
ここで$\hat{\sigma}_k$は2群をプールして求める(ATTの場合は処理群だけを使って求める)。目安として、$SB_k$が0.2とか0.25とかを下回っていることが求められる。 - Kolmogorov-Smirnov統計量。$k$番目の共変量の条件$t$における重み付き経験分布関数を以下のように定義する：
$\displaystyle EDF_{tk}(x) = \frac{\sum_i w_i[t] T_i[t] I(X_{ik} \leq x)}{\sum_i w_i[t] T_i[t]}$
で、2本の経験分布関数の差をとって
$KS_k = sup_x | EDF_{1k}(x) - EDF_{0k}(x) |$
ATTの場合、処理群の重みをすべて1とする。この指標はサンプルサイズに依存するので基準を設けにくいのだが、サンプルサイズが中から大のときはだいたい0.1を下回っていてほしい。
次に、処理が多水準でestimandがATEの場合。
傾向スコア推定にあたってのお勧めの方法は、多項のモデルを組むのではなく、ある水準$t$に注目し、「対象者が$t$に属する確率」$\hat{p_t}(\mathbf{X}_i)$を求めるGBMを組む、というのを全水準について繰り返すこと。当然ながら、 $\hat{p_1}(\mathbf{X}_i) + \hat{p_2}(\mathbf{X}_i) + \cdots$は1にならない。でもそんなのどうでもいい。話のポイントは、各水準と全体との間で共変量をバランスさせることなのだ。[←へええええ！]
反復の停止にあたっても「$t$ vs. 全体」での共変量バランスを監視する。標準化バイアスは
$PSB_{tk} = | \bar{X}_{kt} - \bar{X}_{kp} | / \hat{\sigma}_{kp}$
$ \bar{X}_{kp}, \hat{\sigma}_{kp}$ は全群をプールして重み付けなしで算出する。KS統計量は
$KS_{tk} = sup_x | EDF_{tk}(x) - EDF_{pk}(x) |$
$ EDF_{pk}(x)$は全群をプールして重み付けなしで算出する。
なお、バランスの要約統計量を示す際には、PSB, KSの全群を通した最大値を使うとよい。
処理が多水準で、estimandがATTの場合は... [いまあまり関心ないのでパス]
ところで、doubly robust推定というのもあってだね... ウェイティングしても共変量のインバランスは少しは残るわけで、ウェイティングするだけでなく、さらに共変量を投入した重み付き回帰モデルを組むことがある。これを推奨する人もいるし、確かに処理効果の推定はより正確になるらしいんだけど、変数選択しなきゃいけないというのが決定。著者らのお勧めは、基本はウェイティングのみとし、どうしてもdoubly robust 推定したい場合は「まだインバランスが残っている共変量」を実質科学的な観点から選択すること。RCTでも設計段階で共変量を実質科学的に特定するし、事後的に調整するときにあらためて変数選択なんてしないでしょ、という理屈。
最後に、有効サンプルサイズについて。ウェイティングで分散は拡大する。そのインパクトを捉える保守的な指標として、
$\displaystyle ESS_t = \left( \sum_i T_i[t] w_i \right)^2 / \sum_i T_i[t] w_i^2$
を用いる。ここで$w_i$は、ATEなら$1/\hat{p}_t (\mathbf{X}_i)$ね。有効サンプルサイズが小さくなると云うことは、少数のケースにすごいウェイトがついているということ、つまり オーバーラップが十分でないというシグナルである。
やれやれ、疲れた。以上が前半のメモ。
後半は事例紹介。すごく役に立ちそうだが、必要になったときに慌てて読むってことにしよう。[←自分に甘い]
処理の水準数が3以上のときの傾向スコア調整で、共変量から各水準への所属確率を推定するんだけど、その推定はなにも多項ロジスティック回帰のようなひとつのモデルでやらなくても、水準ごとに別々のモデルでよいし、ある対象者について確率の和が1にならなくても別にいいじゃん... というところが意外であった。そういうもんなんすか－。
論文：データ解析(2015-) - 読了：McCaffrey, et al. (2013) 処理の水準数が多いときの傾向スコア推定 by 一般化ブースト回帰