読了:Wang, et al. (2006) 観察データからの因果効果推定に使うIPTW推定量は処理の割付についてのある仮定が破られていると歪むのでその歪みの大きさを推測する方法を考えたぞ

Wang, Y., Petersen, M.L., Bangsberg, D., van der Laan, M.J. (2006) Diagnosing Bias in the Inverse Probability of Treatment Weighted Estimator Resulting from Violation of Experimental Treatment Assignment. Working Papter 211, Division of Biostatistics, University of California, Berkeley.

 仕事の関係でこの1ヶ月近く延々と悩んでいることがあるんだけど、あまりにspecificな問題で、より一般的な問題として捉え直したいもののどう捉えたらいいのかわからず悶々としている。で、なんとジャスト・フィットなタイトルを持つ論文をみつけて大喜びし、アブストラクトは理解不能だったが、勢い込んで読んでみた。
 いや、動機は間違ってなかったと思うんだけど… たしかに私が抱えている問題は、ある種の実験条件の割付の話で、しかし割付は完全には無作為化できておらず、分析にあたって割付確率の逆数でウェイティングしようとしていて、でもそこにはある種のバイアスがあって、それを診断したい、という話なんだけど… 蓋をあけてみたら、求めていたのとはまるきり違う内容で、途方に暮れた。
 意地を張って少しだけ目を通したけど、もうね… 地獄でしたよ…

1. イントロダクション
 観察データから処理・介入の因果効果を推定するというのは重要な問題である。Robinsの周辺構造モデル(MSM)では3つの推定量が提案されている。G-computation推定量、Inverse Probability of Treatment Weighted (IPTW)推定量、そしてDouble-Robust(DR)推定量である。このうちIPTWが一番簡単。
 IPTW推定量の一致性はexperimental treatment assignment(ETA)仮定に依存している。ETA仮定とは、個々の処理・介入はある確率で起きる、被験者の過去は関係ない、という仮定である。ETA仮定が破られていたらIPTW推定量は定義できない。
 現実にはETA仮定は破られている。被験者の共変量のせいである層の処理確率がすごく低くなったり、単なる偶然でそうなったりする。このときIPTW推定量は歪む。しかしそのことについて調べる人は少ない。きっとツールがないからだろう。本研究ではツールを提案します。名付けてETAバイアス診断(DEB)。

 IPTW推定量の一致性は処理メカニズムの推定の一致性にも依存する。しかし現実には、正しくないメカニズムを使わざるを得ないこともある。たとえば

  • 処理を受ける確率がすごく低い人がいるとき。ウェイトが極端になりIPTW推定量の分散が増える。
  • 処理の割付が共変量\(W\)に強く依存しているんだけど、関心ある因果効果にとっては\(W\)はたいした交絡変数じゃないというとき。ウェイト算出に\(W\)を使うとETA仮定が破られてすごいバイアスがかかる、むしろ\(W\)を無視することによって生じる交絡によるバイアスのほうがまし、ということがありうる。[なるほど… 因果効果推定じゃなくて、抽出確率が不均一な標本抽出に伴うウェイティングのときにも、たとえば層別非比例抽出である層にすごいウェイトがかかってしまう、でも実は層別変数は関心ある変数と独立で、むしろウェイティングしないほうがましだわ、という場面がある。似たような話だろうな]

提案するツールはこういう効果の定量化にも役立つぞ。

2. Point Treatment Settingにおける因果推論レビュー
2.1 因果推論の反事実フレームワーク
 無作為標本のデータ構造を $$ X = ((Y(a), a \in \mathcal{A}), W) \sim F_{X0} $$ とする。\(W\)は共変量、\(\mathcal{A}\)は可能な処理の集合、\(a\)はある処理、\(Y(a)\)はある被験者がもし処理\(a\)を受けたらアウトカムがどうなるかという反事実である。
 被験者が実際に受けた処理を確率変数\(A\)とする。\(X\)には欠損が空いていて、その欠損の場所を示しているのが\(A\)である。\(A\)の条件付き確率分布を \(g_0(a | X) \equiv P(A=a|X)\)と書く。[えええ?処理は共変量\(W\)だけじゃなくて反事実の集合\(Y(a)\)にも条件づけられてるってこと…???]
 我々は\(n\)個のiidな観察\(O_1, \ldots, O_n\)を得る。そのデータ構造を $$ O = (A, Y \equiv Y(A), W) $$ とする[つまりたまたま観察されたアウトカム\(Y(A)\)を\(Y\)と略記してるんですよってことね]。観察された分布\(P_{F_{X0}, g_0}\)を\(P_0\)と書く。

 処理メカニズム\(g_0(a|X)\)について、無作為化仮定(RA)と実験処理割付仮定(ETA)を導入する。
 無作為化仮定(RA)とは、すべての\( a \in \mathcal{A} \)について \( g_0(a | X) = g_0 (a | W) \)という仮定である[あ、なるほど… 処理は反事実の集合には条件づけられてないってことね]。欠損データモデルにおいて欠損がMARだという仮定である。
 実験処理割付仮定(ETA)とは、\(W\)のそれぞれの層において、それぞれのありうる処理が正の確率を持つという仮定である。すなわち $$ \min_{a \in \mathcal{A}} g_0(a | W) > 0, \ \ F_{0W}-a.e.$$ [末尾に書いてある\( F_{0W}-a.e. \)というのがどういうことなのかわからない。検索してみたらa.e.とは数学用語「ほとんどいたるところで」の略語だそうだ。しらんがな!わて文系やがな! 見なかったことにしよう。まああれだ、要するに、ある処理を絶対受けない人はいないっていう仮定ね。そりゃそうね、もしそんなことがあったらどう逆立ちしても共変量の交絡を取り除けない]
 後述する安定化ウェイトというのを使うとき、ETA仮定の\(V\)-specificな弱いバージョンというのを仮定するだけで、以下を満たす条件つき分布 \(g(\cdot | V)\)が存在するといえる: $$ \sup_{a \in \mathcal{A}} \frac{g(a | V)}{g(a | W)} < 0, \ \ \ F_{0W}-a.e. $$ [わからん。全くわからん。ここは心を無にして先に進もう…]

 RAとETAの下で、\(O \)の密度を\(F_X\)のパート[共変量と反事実集合の分布のパートってことね]と\(g\)のパート[処理の条件つき分布のパートってことね]に分解できる。$$ p(O) = p(W) p(Y|A,W)g(A|W)$$ $$ = p(W) p(Y(a)|W)|_{a=A} g(A|X) $$ [なんでこんな難しい書き方するんだろう… 私のような素人へのいやがらせだろうか]
 反事実アウトカム\(Y(a)\)における処理\(a\)の周辺因果効果を\(V\)で調整したものは\(E_{F_{X0}}(Y_a|V)\)と定義される。この効果を周辺構造モデルでモデル化すると $$ E_{F_{X0}}(Y_a|V) = m(a, V|\beta_0)$$ ここで関心あるパラメータは真の因果パラメータ\(\beta_0\)である。[なにを云っているのか小難しくてわからんが、共変量ベクトル\(W\)のかわりになるような\(V\)があるとして、そのもとでの、処理\(a\)のもとでの反事実アウトカム\(Y_a\)の、元のデータ構造の分布\(F_{X0}\)を通じた期待値を推測できれば、処理の因果効果を推定したことになるよね。それを推定するモデル\(m\)を組みますよ。その入力は\(a\)と\(V\)。モデルのパラメータを\(\beta_0\)としますね。ってことかしらん]

2.2 因果パラメータの推定
 さて、因果パラメータ\(\beta\)の3つの推定量をご紹介しましょう。

  • G-computation推定量。回帰モデル\(E(Y|A,W)=Q(A,W)\)に基づく。ステップ1, \(Y\)を\(A,W\)に回帰してモデル\(\hat{Q}(A,W)\)をあてはめ、すべての被験者、すべての\(a \in \mathcal{A}\)について \(\hat{Y}_{a,W} = \hat{Q}(a, W)\)を求める。ステップ2, \(\hat{Y}_{a,W}\)を\(a,V\)に回帰する。[ちょ、ちょっと待って… \(V\)って共変量ベクトルの代わりになる、傾向スコアみたいなもんのことかと思ってたんだけど、違うの!? 勉強不足でわからん…]
  • IPTW推定量。\(g(A|W)\)を処理モデル、\(h\)を\(A,V\)の関数とする。\(\epsilon(\beta) = Y – m(A,V|\beta)\)とする。このとき、IPTW推定方程式は$$ D_h(O|\beta, g) = \frac{h(A,V)\epsilon(\beta)}{g(A|W)}$$ 処理モデルが正しく、かつETA仮定が成り立っていれば、これは\(\beta_0\)に対して不偏である。[わ・か・ら・ん… こいつを最小化する\(\beta\)が\(\beta_0\)の不偏推定値になるってこと? この式の意味は、観察されたアウトカム\(Y\)のそのモデル上の予測値\(m(A,V|\beta)\)からの残差に、謎の値\(h(A,V)\)を掛けて、その処理を割り付けられる確率\( g(A|W) \)で割った値を、最小化するような\(\beta\)を求めなさいってこと? というか\( h(A,V) \)ってなによ]
     IPTW推定量は\(Y\)を\(A, V\)に重みつき回帰すれば簡単に求められる。

    • 標準的には\(h(A,V) = \frac{d}{d\beta}m(A,V|\beta)\)。つまり重みは\(w_i = \frac{1}{g(A_i|W_i)}\)。[なにがなんだかわからんが、嗚呼、やっと一瞬だけ理解できる話になった。割付確率の逆数をウェイトにしてWLS推定しなさいってことね]
    • 「安定化ウェイト」というのもある。\(h(A,V) = g(A|V) \frac{d}{d\beta} m(A,V|\beta)\)。つまり重みは\(w_i = \frac{g(A_i|V_i)}{g(A_i|W_i)} \)。[なにがなんだかもうさっぱり]
  • double robust 推定量。[一行たりともわからないので省略。DRって、傾向スコアをウェイトとしても共変量としても使うことだと思ってたんだけどな… とにかく、まるっきりからっきし知識不足であることがわかった]

 [くじけそうなので、とにかくここまでのところで自分にも理解できたことをメモしておこう。処理がアウトカムにもたらす因果効果を推定するために、アウトカムを共変量と処理に回帰する、その際にその人がその処理に割り付けられる確率の逆数をウェイトにしてWLS推定する、というのがIPTW推定量。その前提になっているのは、処理への割付は共変量には依存しているけどアウトカムには依存してないというRA仮定、そして誰だってどの処理にだって割付られる確率はゼロじゃないというETA仮定。ってことでいいんですか???]

3. ETA仮定違犯によるIPTW推定量のバイアス
 IPTW推定量にバイアスが生じる原因を挙げよう。

  • 周辺構造モデル自体の誤指定。ただし、ここでは\(\beta\)を単純に、真の因果パラメータを\(m(a, V|\beta)\)にprojectしたものと捉える。[つまり\(m\)の誤指定についてはここでは考えないことでしょうか]
  • 処理メカニズム\(g(A,W)\)の誤指定。しかしここではこれも考えない。
  • ETA仮定を破っていること。これが本研究の焦点。

 IPTW推定量のバイアスを、真の因果パラメータと、IPTW推定量を有限標本に適用したときの期待値(データ生成分布は真とする)との差として定義する。処理メカニズムが正しく指定されているなら、バイアスは主にETA違犯によるバイアスである。ほかにも、真の割付メカニズム(\(g_{P_0}(A|W)\))を有限標本で推定している(\(g_{P_n}(A|W)\))ことによるバイアスがある。
 \(g_{P_n}(A|W)\)が正しいという仮定の下でのIPTW推定量のバイアスをETA.Biasと呼ぼう。観察データの真の分布を\(P_0\), 経験分布を\(P_n\), 関心ある因果パラメータを\(\Psi(P_0)\), IPTW推定量を\(P_0 \rightarrow \hat{\Psi}_{IPTW}(P_n)\)、として、ETA.Biasは、\(\hat{\Psi}_{IPTW}\)を\(P_n\)に適用したときのバイアス、つまり $$ ETA.Bias(\hat{\Psi}_{IPTW}, P_0, n) = E_{P_0} \hat{\Psi}_{IPTW}(P_n) – \Psi(P_0) $$ である。
 これを次のように分解できる:$$ ETA.Bias(\hat{\Psi}_{IPTW}, P_0, n) = (E_{P_0} \hat{\Psi}_{IPTW}(P_n) – \hat{\Psi}_{IPTW}(P_0)) + (\hat{\Psi}_{IPTW}(P_0) – \Psi(P_0)) $$
 右辺のふたつめ[つまり、IPTW推定値と真値の差ね]は、ETAの理論的違犯によるバイアスである。ある被験者がある処理を受ける確率がゼロの時にはいつでも生じる。たとえば、\(A\)が催奇形の薬物で、妊婦には絶対に与えられない、というような場合に生じる。
 右辺のひとつめは[IPTW推定量の期待値とIPTW推定値の差ね]、\(P_n\)が有限標本から得られていることに由来するバイアスである。\(g(a|W)\)は理論的には非ゼロだけど、標本サイズが減るにつれ、ある層の被験者がみんな、たまたまある処理を受けてない、ということが起きる。標本にヒスパニックが5人いて、たまたまみんな処理を受けてない、とか。

3.1 ETA仮定違犯によるバイアスのブートストラップ推定
 観察データの真の分布\(P_0\)はわからないので、ETA.Biasの推定にあたっては\(\hat{P}_0\)を使うしかない。そこでブートストラップを使う。ブートストラップ標本の経験分布を\(P^\#_n\)として、\(E_{\hat{P}_0} \hat{\Psi}_{IPTW}(P^\#_n) – \Psi(\hat{P}_0) \)を診断ツールとする。
 ブートストラップの手順は次のとおり。

  1. アソシエーションのモデル\( Q_{\hat{P}_0}(A, W) \)を推定して\(E_0(Y|A,W)\)を得る。処理モデル\(g_{\hat{P}_0}(A|W)\)を推定して(これはここでは既知とする)、\(g_0(A|W)\)を得る。そして観察データから\( P(W) \)を得る。[言い方がいちいち難しくて理解できない… あとの数値例でクリアになるだろうか…]
  2. \( \hat{P}_0 \)からブートストラップ標本 \( P^\#_n \)を生成する。で、それぞれの標本について\(\hat{\Psi}_{IPTW}(P_n^\#) \)を求める。これを平均したのを、\( P_0 \)から得た因果パラメータと比べる。

[正直、全くついて行けない。なにをやろうとしているのかさえよくわからなくなってきた。たぶん、「処理への割付がランダム化されていないとき、割付確率をウェイトにしてアウトカムを要因と共変量に回帰する。でも原理的には、ある共変量ベクトルの実現値の下で割付確率がゼロになっていたら交絡は排除できないし、実際、標本において運悪くそういうことが起きていることも多い。それを調べるために、その標本から得た因果効果の推定値と、その標本からさらにブートストラップ法で推定した因果効果の推定値の期待値を比べる」ということだと思うんだけど…]

4. シミュレーション研究
[正直、3節までが難しすぎて心が折れました。8頁まるごと、目もくれずにすっ飛ばす]

5. データ分析
 HIV患者に対する抗レトロウィルス薬は罹患率と死亡率を劇的に下げるんだけど、患者の生涯にわたる複雑な治療が必要で、なんとかしてアドヒアランスを維持しないといけない。そこで二つの介入について検討する: (1)pill box organizer の使用[調べたところ、日別に小分けになってる薬入れのことらしい。へー], (2) a once daily 抗レトロウィルス治療レジメン [一日に何度も薬を飲むんじゃなくて一回で済ませるようにしたら患者は薬を飲み続けるか?って話かな]。[処理変数が2つあるわけだけど、別々に分析している。以下のメモでは(1)のほうだけメモする。面倒なので薬入れと略記する]
 データは[…略]
 ある観察は、ある患者がある月に、薬入れを使ったか(\(A\), 二値)、その月のアドヒアランス(\(Y\), 連続)、前月の29個の共変量(\(W\))からなる。
 \( m(a|\beta) = \beta_0 + \beta_1 a \)と仮定する。関心ある因果効果は\(E(Y_1-Y_0) = \beta_1\)である。[あああ… やっと実感が湧いてきた… こういう単純な例にしてくれないと、わかんなくなるのよ、ほんとに]
 \(g(A|W), Q(A,W)\)をどうやって推定したかというと…[たぶん本筋と関係ないと思うので略。要するに、薬入れの使用が無作為割付されてないので、ある人が薬入れを使う確率のモデルを組み、また薬入れ使用有無と共変量でアドヒアランスを説明するモデルを組んだ、ということであろう]
 因果効果のIPTW推定量を求めた。ウェイトは安定化ウェイト\(g(A)/g(A|W)\)を用いた。[共変量のもとでの割付の条件付き確率の逆数に、割付の周辺確率をかけたものをウェイトとする、ってことか。薬入れを使いそうにない人が使ってくれたら大きい重みを付けるけど、そもそも薬入れを使っている人が多かったら使っている人には大きな重みをつける、ってこと? ふうん…]
 結果。薬入れ使用の因果効果のG-computation推定値、DR推定値、IPTW推定値を出せる。で、提案手法によってそのETA.Biasを推定できる。薬入れではそうでもなかったけど、一日一回処方のほうではIPTW推定量が大きめになり、それはETA.Biasのせいだということがわかりました。

6. 考察
[略]

 … 正直なところ、仕事の役にはまるで立たなかったが、まああれだ、世の中にはこんな問題があるのねということがわかったので、よしとしよう。
 周辺構造モデルってちゃんと勉強した方が良さそうだ。いつか仕事の役に立つ日が来るかもしれない。この論文は書き方が難しくてお手上げだったけど、簡単な具体例とともに解説していただければ、ある程度は理解できるんじゃないかという気がしてきた(根拠のない楽観)。