読者です 読者をやめる 読者になる 読者になる

2×2×k表の共通オッズ比の検定

共通オッズ比の検定の帰無仮説は、層の数を $k$ として第 $i$ 層のオッズ比を $\theta_i$ とすると、
\[
H_0: \psi_1=\psi_2=\ldots=\psi_k=\psi=1
\] \[
H_1: \mbox{not } H_0
\] 以下である。
これは、オッズ比均質性の検定よりは厳しい帰無仮説になってまいす。
この共通オッズ比 $\psi$ は、今の時代なら最尤推定量 $\hat{\psi}$ を計算機で解くことによって推定してしまう事もできます (解析的には解けない)。
昔はこれも困難だったので、有名な Mantel & Haenszel (1959) の近似推定量がよく用いられていたようで、今でも頻用されていると思います (紙さえあれば計算できるほど簡単だし)。
これらのベースになっているのが、二つの二項分布の積のモデルです。
このモデルにおいて、各層のパラメータをオッズ比に変換し、局外母数に対する十分統計量を与えたもとでの条件付き分布を求めると、非心超幾何分布の積モデルが得られます、非常に綺麗です。
非心超幾何分布の積モデルにおいて、各層のオッズ比が共通でその値が $1$ になる、という仮説が上で挙げた帰無仮説に対応しています。
この場合モデルは、超幾何分布の積モデルになります。

とりあえず今回は手始めに、SAS の logistic procedure でいろいろ見てみます。
共通オッズ比の検定は、$y$ がイベント、$x$ が群変数、$z$ が層別変数としたロジスティックモデル
\[
\mathrm{logit}\psi_i=\alpha+\beta x_i+\gamma z_i
\] において、$\beta=0$ の検定に帰着します。
$z$ が与えられたという条件の下で、$\beta=0$ ならばオッズ比は $\exp(\beta=0)=1$ になるからです。

data hor1;
  call streaminit(041309);
  do i = 1 to 200;
    do x = 0 to 1;
    do z = 0 to 1;
      /* odds ratio */
      /* x: OR = exp(1) = 2.7182 */
      /* x * z: OR = exp(2) = 7.3891 */
      /* logOR = exp(b), logit1 - logit0 = b */
      /* 上より設定するlogitにはマイナスをつける */
      logit = -(-1 * x + 1 * z + 2 * x * z);
      p = 1 / (1 + exp(-logit));
      y = rand('Bernoulli', p);
      output;
    end; end;
  end;
run;

/* オッズ比が均一ではない事を考慮しない場合 */
proc logistic data = hor1;
  class x z / param = ref ref = first;
  model y = x z / clodds = both;
  exact x / estimate = both;
run;

proc freq data = hor1;
  tables z * x * y / nopercent cmh;
  exact comor;
run;

/* オッズ比が均一ではない事を考慮した場合 */
proc logistic data = hor1;
  class x z / param = ref ref = first;
  model y = x z x * z / clodds = both;
  exact x / estimate = both;
run;

/* 共通オッズ比=1を設定した場合 (帰無仮説) */
data hor2;
  call streaminit(041309);
  do i = 1 to 200;
    do x = 0 to 1;
    do z = 0 to 1;
      /* OR: exp(0) = 1 */
      logit = -(0 * x + 1 * z);
      p = 1 / (1 + exp(-logit));
      y = rand('Bernoulli', p);
      output;
    end; end;
  end;
run;

proc logistic data = hor2;
  class x z / param = ref ref = first;
  model y = x z / clodds = both;
  exact x / estimate = both;
run;

proc freq data = hor2;
  tables z * x * y / nopercent cmh;
  exact comor;
run;

/* 共通オッズ比=exp(1)を設定した場合 (対立仮説) */
data hor3;
  call streaminit(041309);
  do i = 1 to 200;
    do x = 0 to 1;
    do z = 0 to 1;
      /* OR = exp(1) = 2.7182 */
      logit = -(1 * x + 1 * z);
      p = 1 / (1 + exp(-logit));
      y = rand('Bernoulli', p);
      output;
    end; end;
  end;
run;

proc logistic data = hor3;
  class x z / param = ref ref = first;
  model y = x z / clodds = both;
  exact x / estimate = both;
run;

proc freq data = hor3;
  tables z * x * y / nopercent cmh;
  exact comor;
run;

オッズ比が均一ではない事を考慮した場合としない場合の二種類、共通オッズ比 $=1$ と $\exp(1)$ の計4パターンを示してある。

オッズ比が等しくないことを考慮しない場合に Mantel-Haenszel 検定を行った結果を見ると、検定では有意にならない、つまり共通オッズ比 $=1$ を棄却する事ができない。
検定としては妥当だが、検出力のロスがあるような状況にあたる。
各層のオッズ比をかなり高めに設定しているので、検定の結果を鵜呑みにすると危ない状況になる。
まあそういう判断をすることはあり得ないと思うが。

オッズ比が等しくないことを考慮した場合、$\delta$ を見ればオッズ比均質性にあたるが、今回は $\beta$ を見る。
この場合はちゃんと $\beta=1$ を棄却できる、だが、もはや上で挙げた意味での共通オッズ比の検定にはならない。

共通オッズ比 $=1$を設定した場合は帰無仮説なので、これを棄却できる確率は有意水準に等しい。

共通オッズ比 $=\exp(1)$ の場合は対立仮説なので棄却できる。

乱数なので同じ結果が得られないかもしれないが、大体こんな感じの結果になると思います。

現実的には、オッズ比が等しくないことを考慮しない場合のロスがどの程度で、どんな場合にやばくて、どうしたら誤った判断をしないで済むのかを、解析する前に頭にたたき込んでおくと良い気がします。

ついでに、Logistic Procedure の Exact Statement では Two-sided = 2 * One-sided になっている模様? (9.1.3)。
$\Pr(\Pr(\mathrm{Test Statistics}) > \Pr(t))$ か、$\Pr(|\mathrm{Test Statistics}-\mathrm{E}_0(T)| \geq |t-\mathrm{E}_0(T)|$ を使った方がいいと思う。
ものすごいついでに、Exact Confidence Limits for the Common Odds Ratio の計算も面白いので見てみるといいかも。

追記

2007/05/03に書いたものを手直し?したものです。