O’Gorman, T.W, Woolson, R.F., Jones, M.P. (1994) A Comparison of Two Methods of Estimating a Common Risk Difference in a Stratified Analysis of a Multicenter Clinical Trial. Controlled Clinical Trials, 15, 135-153.
仕事の都合で読んだ。層別された2×2クロス表について層を潰したリスク差を求めるとき、マンテル・ヘンツェルのアプローチだとどうなるか、という話。
MHオッズ比についての解説はその辺の教科書に載っているけど、リスク差についての解説は少ないので、仕方なくめくった。勤務先の本棚にあるAgrestiの厚い本に書いてあったような気がするけど、いま自宅に閉じこもっているもので…
いわく。
二値の目的変数があり(有害事象の生起とか)、\(K\)個の層があって(施設とか)、層\(k\)における処置条件の患者数を\(n_k\), 統制条件の患者数を\(m_k\), あわせて\(N_k = n_k + m_k\)とする。目的変数が1の患者数を\(x_k, y_k, t_k=x_k + y_k\)とする。
[2021/05/30 追記: わかりにくいので表にしておく。MathJaxのarray環境で表組みしているんだけど、日本語のフォントをふつうのフォントにするのはどうしたらいいんだろう…] $$ \begin{array}{c|cc|c} \hline & 目的変数が1 & 目的変数が0 & 計 \\ \hline 処置条件 & x_k & n_k – x_k & n_k \\ 統制条件 & y_k & m_k – y_k & m_k \\ \hline 計& t_k & N_k – t_k & N_k \\ \hline \end{array}$$
こういうときはMantel-Haenszelオッズ比をよく使いますね。MHオッズ比では, \(x_k = y_k = 0\)の層は無視される。MHリスク比も同様である。しかしリスク差の場合、\(x_k = y_k = 0\)であることは、その層の患者数が十分であれば、リスク差が小さいよという情報を提供しているわけで、これを無視したくない。
[読者がMantel-Haenszelオッズ比についてよく知っていることを自明視しているので、ここで立ち止まってメモしておく。添字\(k\)を省略し、\(p = x/n, q = y/m\)として、オッズ比の推定量は$$ \hat{OR} = \frac{p/(1-p)}{q/(1-q)} = \frac{p(1-q)}{q(1-p)} =\frac{p-pq}{q-pq} $$ですわね。代入して$$ \hat{OR} = \frac{ (x/n)-(x/n)(y/m) }{ (y/m)-(x/n)(y/m) } = \frac{ mx-xy }{ ny-xy } = \frac{x(m-y)}{y(n-x)}$$ である。さてMHオッズ比とは、分母も分子も総数で割ってから足し上げた$$ \hat{OR}_{MH} = \frac{\sum_k x_k(m_k-y_k) / N_k}{\sum_k y_k(n_k-x_k)/N_k} $$である。なるほど、\(x_k = y_k = 0\)の層は無視されますね。
なお、MHオッズ比と各層のオッズ比との間には、上の式の分母の総和記号の内側を\(S_k = y_k(n_k-x_k)/N_k\)として、すべての層で\(S_k \neq 0\)であれば$$ \hat{OR}_{MH} = \frac{\sum_k S_k (\hat{OR}_k)}{\sum_k S_k} $$ という関係があるのだそうだ。つまり\(S_k\)で重みづけて平均しているわけね。さらに、\(S_k\)は\(H_0: OR=1\)の下での\(\hat{OR}\)の漸近分散推定量の逆数になっているのだそうだ。へー。以上、佐藤・高木・柳川・柳本(1998)よりメモした]
以下では処置条件のリスクを\(p_k\), 統制条件のリスクを\(q_k\)とする[原文では\(p_{1k}, p_{0k}\)だけどうざいので書き換えてメモする]。リスク差を\(\delta_k = p_k – q_k\)と書く。\(d_k = x_k/n_k – y_k/m_k\)はその不偏推定量である。
たいていの分析ではオッズ比が層間で共通だという仮定を置く。これはリスク差がロジットスケール上で共通だと仮定しているのと同じである。
層間でオッズ比が一定かどうかとリスク差が一定かどうかとはちょっと話がちがう。たとえば、処置条件での有害事象の生起確率が処置のせいで一定 \(p_e\)になっているとしよう(\(p_e\)をShepのrelative differenceという)。処置条件の患者に占める、「ああ、処置を受けなかったら有害事象は起こらなかったろうに処置を受けたばっかりに有害事象のリスクにさらされるとは」さんの割合は\(1-q_k\)である。よって\(p_k = q_k + p_e(1-q_k)\)である。リスク差は\(\delta_k = p_e(1-p_k)\)となる。このとき、\(q_k\)がどの層でも小さければリスク差はほぼ一定である。いっぽうオッズ比は一定にならない。
以下では、\(q_k\)がどの層でも小さく、\(p_e\)がほぼ一定の場合についてのみ考える。このときリスク差は、処置のせいで生じた有害事象の割合の近似として解釈できる。もしリスク差が層間で異なるなら、我々の議論は、標本における層別リスク差の平均を推定しているだけである[解釈は難しいよってことね]。そういうときはランダム効果モデルで、リスク差の等質性を検定したうえで母リスク差を推定した方がよろしい。
では本題。共通リスク差の推定において \(\sum_k W_k d_k / \sum_k W_k\)という推定量を使うとき、\(W_k\)をどうしたらよいか。
推定量その１，WLS。最初に全セルに1/2を足しておいて、$$ W_{WLS,k} = [x_k(n_k – x_k)/n^3_k + y_k(m_k -y_k)/m^3_k]^{-1} $$ \(d_k\)の計算にも全セルに1/2を足した奴を使う。
WSL推定量は、各層のリスク差を、その分散の逆数\(\sum W_{WLS,k}\)で重みづけたものである。不偏推定量にはならない。
推定量その２，Cochran-Mantel-Haenszelウェイト。$$ W_{CMH, k} = n_k m_k / N_k $$とする。推定量の分散は…[略]。
[他の方法もいくつか紹介してあるけど、面倒くさくなったので略]
。。。ここまでで5ページ。この論文はここから、リスク差のWLS推定量とCMH推定量をシミュレーションで比較するんだけど、なんだか疲れてきちゃったので読まなかった。アブストラクトによれば、WLS推定量は大きく歪む、処置群と統制群の人数が違ってて標本サイズが小さいときはどちらの推定量の信頼区間もあてにならない、らしい。
ほんとに知りたかったのは、複雑なデザインの横断調査から得た標本で、標本に確率ウェイトがついちゃっているときにどうやって共通リスク差を求めるか、ということなんだけど… おとなしくランダム効果モデルを組めってことでしょうね。そりゃそうか。