読了：Rudolph et al. (2014) 大規模調査の標本の一部について別のデータがあるとき、そこで推定された平均処理効果を母集団へと一般化する方法

Rudolph, K., Diaz, I., Rosenblum, M., Stuart, E. (2014) Estimating Population Treatment Effects From a Survey Subsample. Americal Journal of Epidemiology, 180(7), 737-748.

　これ仕事の役に立つんじゃないかしらんと思って読んでみた奴。Google様的な引用件数は20。
　自分の仕事に近づけて言うと、えーっと、大規模な消費者調査のデータがあり、そのなかの一部の対象者についてだけ広告接触有無と製品購買有無がわかっているとき、母集団における広告効果を推定したい、というような話である。RCTの結果を一般化するんじゃなくて観察研究の結果を一般化するというのがポイント。

(イントロ)
　この論文の動機は、国の代表標本による調査NSC-Aの標本からさらに抽出した青少年についてバイオマーカーのデータを得たので、そこから、非無作為な処理(居住地)がコルチゾール勾配(インタビュー中のコルチゾールレベルの減少率)に与える効果を推定できないか、という問題にある。つまり、利用可能なデータと国レベルの代表標本をつなぎ合わせて、結果を国全体へと一般化したい。
　RCTの結果を目標母集団に一般化するという先行研究はある[Stuart, et al.(2011 JRSS), Cole & Stuart(2010 Am.J.Epid.)というのが挙げられている。前者は傾向スコアでもって結合する模様]。しかし観察研究の結果を目標母集団に一般化するという研究は少ない。

  • 二重頑健法はある種のモデル誤指定のもとでも一致性を持つので、非無作為処理割付や非無作為選択とか右側打ち切りとか欠損とかを調整するために用いられてきた。しかし推定量の実装は難しいことがある。
  • もっとシンプルなHorvitz-Thompson逆確率重みづけ(IPW)推定量のほうがよく使われているが、有効性と頑健性に問題がある。
  • 最近ではtargeted ML推定(TMLE)というのがあるけど、ウェイト付きサーヴェイデータという文脈でTMLEを使った研究はみあたらない。[へー]

さらに、これらの手法についての議論は生物統計学ではみかけるけど疫学ではみかけない。

方法
　次のシナリオを考える。サーヴェイにおいて個人が選択される。選択確率は既知。処理についての情報と共変量はサーヴェイで選択されたすべての参加者について完全に観察されている[疫学の論文らしく処理treatmentと書いているけど、ただの観察変数]。しかしアウトカムは標本のサブセットについてしかわかっていない。
　ベースライン共変量のベクトルを\(W\)とする。処理の有無を示す二値変数を\(A\)とする。サーヴェイの標本に入ったかどうかを示す二値変数を\(D_1\)とする[原文では\(\Delta_{svy}\)なんだけどうざいので略記する]。サブ標本に入ったかどうかを示す二値変数を\(D_2\)とする[原文では\(\Delta_{sub}\)]。関心あるアウトカムは連続量で、これを\(Y\)とする。潜在アウトカムの用語を使って、個人\(i\)の\(A=0,1\)の下での[潜在]アウトカムをそれぞれ\(Y_{0i}, Y_{1i}\)とする。
　ここで考えるエスティマンドは目標母集団を通じた平均処理効果 PATE \(E(Y_1-Y_0)\)である。他の種類のATEを考えることもできるけど(サーヴェイ標本ATE \(E(Y_1-Y_0|D_1=1) \)とかね)、そっちの議論は付録をみてくれ。
　PATEの3種類の推定量を比較するぞ。Rのコードは付録をみてくれ。[おお、ありがたいね]

　その1、IPW推定量。[以下、ウェイトが上添字で\(w^{hoge}\)という風に表記されているんだけど、みにくいのでカッコに囲んで\(w (hoge)\)と書く]
　サーヴェイ選択確率の逆確率ウェイトは $$ w(D_1=1) = \frac{1}{P(D_1=1|W)} $$
　処理確率の逆確率のウェイトは $$ w(A=a|D_1=1) = \frac{I(A=a)}{P(A=a|D_1=1, W)} $$ [ちょ、ちょっと待って… 分子はいったいなんなの？ 処理\(A\)が水準\(a\)であるときに1, そうでないとき0である変数？ だとしたら、反事実的な処理水準についても割付確率は0でないのにウェイトは0になるってこと？ あ、処理の水準が0であるサーヴェイ参加者について実際に使う個体ウェイトは\(w(A=0|D_1=1)\)だけで、反事実的なウェイトは使わないから0でいいじゃんってことかな…]
　サブ標本選択確率の逆確率ウェイトは $$ w(\Delta_{sub}=1|A=a, \Delta_{svy}=1) = \frac{I(\Delta_{sub} = 1)}{P(\Delta_{sub}=1|A=a,\Delta_{svy}=1,W}$$ [えええ… サブ標本選択確率も処理に依存してるってこと？ そうか、郊外に住んでる青少年のバイオマーカーデータは集めにくいというような交絡も考えるのか]
　以上を掛けまして、$$ w(A=a, D_1=1, D_2=1) = \frac{I(A=a, D_2 = 1, D_1 = 1)}{P(A=a, D_2 = 1, D_1 = 1|W)} $$ ウェイトはロジスティック回帰で推定できる (たとえば\(P(A=a|D_1=1,W)\)はサーヴェイ参加者における処理を共変量で説明すればよい)。[おいおい、\(P(D_1=1|W)\)はどうやって推定すんの？ まあどうにかなるのかな、センサスとかで\(W\)の母集団分布がわかってれば]
　というわけで、PATEのIPW推定量は、サブ標本のサイズを\(r\)として $$ \hat{PATE} = \frac{\sum_i^r Y_i w_i(A=1, D_1=1, D_2 = 1)}{\sum_i^r w_i(A=1,D_1=1,D_2=1)} – \frac{\sum_i^r Y_i w_i(A=0, D_1=1, D_2 = 1)}{\sum_i^r w_i(A=0,D_1=1,D_2=1)} $$

　その2, TMLE。Rのtmleパッケージをちょっと修正しました。付録をみろ。おおまかにいえばこうだ。

  1. サブ標本について\(A\)と\(W\)による\(Y\)の予測値\(\hat{Y}^0\)を得る。ここでは線形回帰を使うけど機械学習とかでもよろしい。
  2. すべての個人\(i\)について$$ H_i = A_i w_i(A=1,D_1=1,D_2=1) – (1-A_i) w_i(A=1,D_1=1,D_2=1) $$ とする。[その人の実際の処理水準に即したほうのウェイトを\(H_i\)って書きますよってことね]
  3. \(Y\)を\(H\)に回帰するんだけどその際\(\hat{Y}^0\)をオフセット項にする。推定した係数を\( \hat{\beta} \)とする。G-計算によって、サーヴェイ標本について、\(A=1\)に割り付けられた下での反事実的アウトカムの予測値と、\(A=0\)に割り付けられた下での反事実的アウトカムの予測値との差を求める。この差をサーヴェイ標本を通じ、サーヴェイウェイトをウェイトにして合計する。

　G-計算についてご存じない人のために簡単に説明すると、これは共変量の周辺分布を使う方法で、疫学ではよく死亡率を標準母集団の年齢分布で標準化するじゃないですか、あれの拡張である。[うぐぐぐ… まあいいや、いずれ別の機会に勉強しよう]

　その3, 二重頑健WLS(DRWLS)推定。

  1. サブ標本について\(Y\)を\(A\)と\(W\)に回帰するんだけど、\(w(A=a,D_1=1,D_2=1)\)をウェイトにしてWLS推定する。G-計算でもって、サーヴェイ標本に標準化した反事実アウトカムを予測する。
  2. これを使って下位標本ATEを予測する
  3. これをサーヴェイウェイトをウェイトにして合計する。

[ああ… G-計算というのはよくわからんが、やろうとしていることはなんとなくわかる…]
　ここで二重頑健といっているのは、アウトカムモデルと、IPW推定量のところで説明した\(w_i(A=a,D_1=1,D_2=1)\)のモデル、どっちかが正しければ一致推定量になるからで、その意味ではTMLEも二重頑健である。
　
シミュレーション研究
[力尽きたのでまるごとスキップ]

ケース研究
[もちろんまるごとスキップ。おじさんはね…疲れているんですよ…]

考察
　というわけで、処理の効果に異質性があり、処理の割付が無作為でなく、2段階の抽出プロセスがあって、プラクティカルなposivity violationがある[ウェイトにでかい値があるってことらしい]ときの、PATEの推定量について検討した。
　シミュレーション研究では、モデルが正しいときも、誤指定があるときも、IPW推定量よりDRWLS推定量とTMLEのほうがMSEが低かった。DRWLS推定はバイアスと分散がもっとも小さかった。RのtmleパッケージをつかってうまいことTMLEを求めるコードを書いたから付録を読め。
　本研究の限界: PATEでは処理効果の異質性は表現されない[そらそうだ]。positivity violationがあるとき推定は難しい。敏感性分析をお勧めする。Robins et al.(2007 Stat.Sci.), Petersen et al.(2012 Stat.Methods.Med.Res.)をみろ。
　IPW推定量の成績がしょぼいというのは新しい話じゃないけど、疫学では広く使われている。みなさん、DRWLSやTMLEだって実装は簡単だし、IPWよかMSEが良いのよ。このことが広まってみんなもっと使ってくれるようになるといいな。
　敏感性分析の結果、処理ウェイトと選択ウェイトのうち極端な2%の値をトランケートすると、バイアスも分散も下がるけど、モデル誤指定によるバイアスは増えるという結果になった。[…中略…]うまいトランケートの仕方は今後の課題だ。
　[ほかにもいろいろ書いてあるけど省略…]

　。。。というわけで、おまいらいつまでもIPW推定してんじゃねえよ、DRWLS推定量やTMLEを使え、という論文であった。へえええ、そうなんすか。
　自分の仕事とはなんだか距離があるなあという気がしてきてしまい、途中からぱらぱらめくっただけになってしまった。おいお前、マーケティングリサーチだって観察データからPATEを推測したいこと多いだろうが！ おいそうだよな！？ と問い詰められば、それはまあ、確かに仰せのとおりでございます、って感じなんだけど…. ううむ、しかしなんかこう距離があるなあ… なんでだろう…
　まあとにかく、勉強になりましたです。