読了: Brumback, He (2011) ウェイトつきデータから共通オッズ比を推定するためにあなたがお使いのMHオッズ比は、標本の層サイズが固定されていると考えると一致推定量でない

Brumback, B., He, Z. (2011) The Mantel-Haenszel estimator adapted for complex survey designs is not dually consistent. Statistics and Probability Letters, 81, 1465-1470.

 仕事の都合で致し方なく読んだ奴について記録しておくけど、正直、面白くも何ともない話である。(著者の先生、すいません…)

 いわく。
 complex survey design [きちんとした定義はわからないのだが、要は個体に確率ウェイトがついちゃっている状況のことであろう]において、共通オッズ比のMantel-Haenszel推定量を求めるにはどうするか。
 Graubard et al.(1989 Biometrics) がその方法を提案していて、SUDAANにも載っている。この推定量は、(1)large-strata limiting モデル(表の数を固定したとしたとき、個々のセルのサイズが上界なく増える) においては一致推定量である。また (2) sparse-data limiting model (表の数が増えるとき、セルサイズが上界を持つ) においても一致推定量である… といわれている。
 本研究は、(1)は正しいけど(2)は正しくないことを示します。

 [先行研究… 省略]

 クラスタ\(i(=1,\ldots,K)\)における個人\(j(=1,\ldots,L_i)\)の曝露を\(X_{ij}\)、アウトカムを\(Y_{ij}\)とする。\(X_i \equiv (X_{i1}, \ldots, x_{i L_i})\)とする。
 共通オッズ比 \(\exp(\beta)\) を仮定し、有限母集団モデル$$ logit (E(Y_{ij}|X_i, b_i)) = X_{ij} \beta + b_i$$を考える。\(b_i\)はクラスタの切片。残差は\(X_i, b_i\)の下で独立とする。
 この有限母集団は超母集団からの無作為2段標本だと考える。

 個人の選択確率の逆数に比例した値を\(W_{ij}\)とする。また、クラスタ\(i\)におけるすべての個体の同時選択確率の逆数に比例した値を\(W_i\)とする。たとえば、あるクラスタが抽出されたとき、そのクラスタにおける個人が互いに独立に抽出されるなら、クラスタ\(i\)が抽出される確率を\(p_i\)、その下で個人\(j\)が抽出される条件付き確率の逆数を\(W_{j|i}\)として、$$ W_{ij} = W_{j|i} / p_i$$ $$ W_i = \frac{\prod_j W_{j|i} }{p_i}$$ である。
 この研究では、クラスタ数は無限で、各クラスタに含まれる個とは有限で比較的小さい場面について考える。たとえば、クラスタがセンサスの調査ブロックの小さな集合で個体がヒトであるような場合である。

 クラスタ\(i\)について2×2クロス表をつくる。ウェイト付きの値を以下とする。
 [そのままメモするとわかりにくいので表にしてみた。略記したが、表中の総和記号はすべて\(j=1\)から\(l_i\)までを回る。\(T^w\)の下添字の付け方が独特なので注意。一桁目がアウトカム, 二桁目が曝露で、1がon, 2がoffだと思う] $$ \begin{array}{c|ll|l} \hline & Y_{ij} = 1 & Y_{ij} = 0 & \\ \hline X_{ij} = 1 & T^w_{11i} = \sum W_{ij} X_{ij} Y_{ij} & T^w_{21i} = \sum W_{ij} X_{ij} (1-Y_ij) & \\ X_{ij} = 0 & T^w_{12i} = \sum W_{ij} (1-X_{ij}) Y_{ij} & T^w_{22i} = \sum W_{ij} (1-X_{ij}) (1-Y_ij) & \\ \hline & & & T^w_i \\ \hline \end{array} $$
 さて、SUDAANに載っているMHオッズ比はこうである。$$ MHOR = \frac{\sum_i^k \frac{T^w_{11i} T^w_{22i}}{T^w_i}}{\sum_i^k \frac{T^w_{12i} T^w_{21i}}{T^w_i}} $$
 large-strata limiting modelでは、抽出されたクラスタ数が母集団におけるクラスタ数に等しい。つまり\(k = K\)である。[ああ、やっと云ってる意味がわかった。漸近的性質を考える際、なにが無限大に近づくとみるのかという話をしているのね。ここでは、クラスタが全数抽出でクラスタサイズが無限大に近づいたらどうなるのかという話をします、ということだ]
 このとき、MHORは\(\exp(\beta)\)の一致推定量である。なぜなら、\(\frac{T^w_{yxi}}{T^w_i}\)はクラスタ\(i\)における対応表の\(y\)行\(x\)列のセルの確率の一致推定量であり、全クラスタの\(T^w_i\)の合計を\(T^w\)として、\(T^w_i / T^w\) は母集団におけるクラスタ\(i\)のprevalenceの一致推定量だからだ。
 ではsparse-data limiting modelではどうか。今度は、\(K\)は無限大、\(k\)は無限大に近づき、\(l_i\)は固定である。[はいはい。今度は、クラスタサイズを変えずにクラスタ数だけ無限に近づいたらどうなるのかという話をしますってことね]
 Breslow(1981 Biometrika) いわく、その場合でも\(W_{ij}\)が一定であれば一致推定量である[個体の抽出確率が均等だってこと?]。いっぽう本論文では、\(W_{ij}\)が一定でないとき、sparse-data limiting modelの下で一般にMHORは\(\exp(\beta\)の一致推定量でないと論じる。

 本論文ではもう一つの推定量を提案する。ウェイティングしてない人数を\(T_{11i}, T_{12i}, \cdots\)と書いて、$$ AMHOR = \frac{ \sum_i^k (\frac{W_i}{N_i}) (\frac{T_{11i} T_{22i}}{T_i}) }{\sum_i^k (\frac{W_i}{N_i}) (\frac{T_{12i} T_{21i}}{T_i})} $$ AMHORはある種の正値条件の下で、sparse-data limiting modelの下で超母集団の\(\exp(\beta)\)の一致推定量になる。証明は付録をみよ。
 では正値条件とはなにかというと…
 クロスセクショナル・デザインの場合、正値条件とは、固定された\(l_i\)について、クラスタ\(i\)から得たすべての可能なjoint sampleが正の確率を持つということである[どの二人をみても、同時に抽出される可能性が存在するってことね]。たとえばクラスタが世帯を含んでいて、世帯あたり一人しか選ばないという場合には、この条件は破られている。このとき、この条件が破られていることによってバイアスが生じるかどうか考える必要が生じる。
 前向き研究のように、標本抽出が\(X_{ij}\)に条件づけられているとき、正値条件とは…[略]。
 後ろ向き研究のように、標本抽出が\(Y_{ij}\)に条件づけられているとき、正値条件とは…[略]。

 AMHORの標準誤差を求めるには…[とかなんとか。略]

 [シミュレーション。読まずにスキップ]

 [数値例。どのクラスタも同じ構造を持っているけど、各クラスタから二人ずつ抽出していくと、MHORの期待値が真の共通オッズ比にならない、という例を示している。なるほど?]

 考察。
 この論文で示したのと似た指摘がconditional logistic回帰についてもなされていて… [略]
 云々。

 。。。目を通しては見たけれど、正直いって論文の主旨にはあまり関心がなくて、ウェイトつきデータでMHオッズ比を求める方法(この論文で言うMHORの式)が知りたかっただけであった。著者の先生方が想定している場面と異なり(具体的にはNational Health Interview Surveyのことを考えておられる模様)、私の場合は、クラスタが全数抽出でクラスタサイズが無限大に近づくと考えたときに共通ORの一致推定量であるだけで御の字である。
 なあんだ、やっぱり単純にウェイト合計を頻度とみればよかったのか。SASのドキュメントやCRANを検索してみたけど、式が見当たらなかったのである。そうか、こういう場合はSUDAANのドキュメントを調べれば良かったんだなあ。