Mantel-Haenszel検定をRで実施する - 統計学入門一歩先へ

はじめに

今回は、以前紹介したMantel-Haenszel検定（MH検定）をRで実施する方法を整理していく。
mstour.hatenablog.com

前も書いたようにCochran-Mantel-Haenszel（CMH検定）と呼ばれることが多いが、Cochranの検定とMantel-Haenszelの検定とは周辺和の固定について厳密には異なるので、本記事ではMH検定という言い方で通すことにする。

以前の紹介記事では2×2表の場合しか取り上げなかったため、今回の事例も2×2表のみとする。今後一般のR×C表の場合の解説を書いた際にはRでの実行例も紹介していきたい。

Rでの実行例

SASという統計ソフトでのMH検定の説明に用いられている例を使うことにする。
https://documentation.sas.com/?cdcId=pgmsascdc&cdcVersion=9.4_3.5&docsetId=procstat&docsetTarget=procstat_freq_examples07.htm&locale=en

f:id:mstour:20201008210722j:plain
この例は、片頭痛の治療薬の効果をプラセボと比較した臨床試験の結果である。ここでは性別で層別した場合の結果を考える（余談だが、治療効果に大きく影響を与えることが確実と思われる因子については、治療の割付においても層別化を行うことが多い）。結果は次の分割表にまとめられる。

女性

女性(female)	改善(Better)	不変(Same)	合計(Sum)
治療薬(Active)	16	11	27
プラセボ(Placebo)	5	20	25
合計(Sum)	21	31	52

男性

男性(male)	改善(Better)	不変(Same)	合計(Sum)
治療薬(Active)	12	16	28
プラセボ(Placebo)	7	19	26
合計(Sum)	19	35	54

まずはデータセットを作成しよう。

st1 = data.frame(Gender = rep("female", 52),
                 Treatment = c(rep("Active", 27), rep("Placebo", 25)),
                 Response = c(rep("Better", 16), rep("Same", 11),
                              rep("Better", 5), rep("Same", 20)))
st2 = data.frame(Gender = rep("male", 54),
                 Treatment = c(rep("Active", 28), rep("Placebo", 26)),
                 Response = c(rep("Better", 12), rep("Same", 16),
                              rep("Better", 7), rep("Same", 19)))
df = rbind(st1, st2)

分割表は例えば以下のようにして出力する。

addmargins(table(df[df$Gender == "female", c("Treatment", "Response")]))
addmargins(table(df[df$Gender == "male", c("Treatment", "Response")]))

出力結果は次のようになった（上は女性、下は男性）。

         Response
Treatment Better Same Sum
  Active      16   11  27
  Placebo      5   20  25
  Sum         21   31  52

         Response
Treatment Better Same Sum
  Active      12   16  28
  Placebo      7   19  26
  Sum         19   35  54

性別(Gender)で層別化したMH検定は次のようにして行う。関数mantelhaen.testは基本機能であるstatsパッケージの関数なので、特にパッケージのインストール・読み込みは必要ない。

mantelhaen.test(x = df$Treatment, y = df$Response, z = df$Gender)

xとyに関連性を見たい変数を、zに層別化する変数を入力する。デフォルトの設定での結果は以下のようになる。
f:id:mstour:20201008210744j:plain
デフォルトでは、「連続修正」をした検定の結果が出力される。連続修正をしない場合にはcorrect = Fと指定する必要がある。MH検定はサンプルサイズがある程度大きい場合の漸近的な性質を利用しているため、サンプルサイズが小さい場合には連続修正をした方が近似がよくなると言われている。一方、サンプルサイズが小さい場合の対応として、exact = Tと指定することによって正確な検定を行うこともできる。
（なお、正確なMH検定がどのような方法なのかはよくわからなかった）

連続修正を行わない場合は以下のようになる。連続修正を行った時と比べて、検定統計量がかなり大きくなりp値が小さくなっているのがわかる。
f:id:mstour:20201008210758j:plain
出力の下側には、各層での共通オッズ比（MH検定は、各層でオッズ比が同じであることを仮定している）のMantel-Haenszel推定量と、その信頼区間が表示されている。信頼区間の計算方法はR documentationに明示されていなかったが、SASで使用されているRobins, Breslow, and Greenland (1986)の方法と結果が一致するので、同様なのだと思う。