タイタニック号沈没事故,(Cochran-)Mantel-Haenszel検定,Simpsonのパラドックス

Titanic データ

豪華客船タイタニック号は,1912年,イギリスからアメリカに向かう処女航海で,氷山に衝突し,沈没した。詳しくは,ウィキペディアのタイタニック号沈没事故や,Encyclopedia Titanica を参照されたい。

タイタニック号の生存統計(諸説あるがそのうちの一つ)がRの Titanic というデータセットに収められている。このデータ構造を str() 関数で調べてみよう:

str(Titanic)
 'table' num [1:4, 1:2, 1:2, 1:2] 0 0 35 0 0 0 17 0 118 154 ...
 - attr(*, "dimnames")=List of 4
  ..$ Class   : chr [1:4] "1st" "2nd" "3rd" "Crew"
  ..$ Sex     : chr [1:2] "Male" "Female"
  ..$ Age     : chr [1:2] "Child" "Adult"
  ..$ Survived: chr [1:2] "No" "Yes"

これからわかるように,Titanic は4個の添字を持つ人数の表($4 \times 2 \times 2 \times 2$ の表)で,各添字は,乗客クラス(ファーストクラス,セカンドクラス,サードクラス,乗組員),性別(男,女),年齢(子ども,大人),生存(No,Yes)を表す。

sum(Titanic) # 全員
[1] 2201
sum(Titanic[,,,"Yes"]) # 生存者
[1] 711
sum(Titanic[,,,"No"]) # 死亡者
[1] 1490

これと同じことは margin.table() 関数でもできる:

margin.table(Titanic)
[1] 2201
margin.table(Titanic, 4)
Survived
  No  Yes 
1490  711 

margin.table() は次のようなこともできる:

margin.table(Titanic, c(2,4))
        Survived
Sex        No  Yes
  Male   1364  367
  Female  126  344
margin.table(Titanic, c(2,4,1))
, , Class = 1st

        Survived
Sex       No Yes
  Male   118  62
  Female   4 141

, , Class = 2nd

        Survived
Sex       No Yes
  Male   154  25
  Female  13  93

, , Class = 3rd

        Survived
Sex       No Yes
  Male   422  88
  Female 106  90

, , Class = Crew

        Survived
Sex       No Yes
  Male   670 192
  Female   3  20

これを見れば,全体として,男性の多くが死に,女性の多くは助かったことがわかる。しかし,サードクラスでは女性も死者のほうが多い。

男と比べた女の生存者の多さは,例えばオッズ比(odds ratio)という指標で表せる(2×2の表,オッズ比,相対危険度 参照)。オッズ比を計算する関数はいろいろなパッケージで定義されているが,自前で作っても,次のように1行でできる:

oddsratio = function(x) (x[2,2]/x[2,1])/(x[1,2]/x[1,1])

これを使えば,全体のオッズ比は

oddsratio(margin.table(Titanic, c(2,4)))
[1] 10.14697

乗客クラスごとのオッズ比は

for (i in 1:4) cat(i, oddsratio(margin.table(Titanic, c(2,4,1))[,,i]), "\n")
1 67.08871 
2 44.06769 
3 4.071612 
4 23.26389 

またはもっとRらしく

apply(margin.table(Titanic, c(2,4,1)), 3, oddsratio)
      1st       2nd       3rd      Crew 
67.088710 44.067692  4.071612 23.263889 

つまり,ファーストクラスでは紳士的に女性の生存が優先されたが,クラスが下がるにつれその傾向が減っていることがわかる。

(Cochran-)Mantel-Haenszel 検定

個数を表す $2 \times 2$ の表,あるいは一般に $I \times J$ の表で,行と列の独立性を調べるには,カイ2乗検定またはFisherの正確検定を使う。例えばタイタニック号での生存が性別によらないという帰無仮説を検定するには,次のようにすればよい:

chisq.test(margin.table(Titanic, c(2,4)))

	Pearson's Chi-squared test with Yates' continuity correction

data:  margin.table(Titanic, c(2, 4))
X-squared = 454.5, df = 1, p-value < 2.2e-16

あるいは

fisher.test(margin.table(Titanic, c(2,4)))

	Fisher's Exact Test for Count Data

data:  margin.table(Titanic, c(2, 4))
p-value < 2.2e-16
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
  7.97665 12.92916
sample estimates:
odds ratio 
   10.1319 

このオッズ比が oddsratio(margin.table(Titanic, c(2,4))) の結果 10.14697 と若干異なるのは,Fisherの正確検定に書いた(conditional MLE あたり)。

しかし,生存と性別の関係は乗客クラスに依存するので,乗客クラスを無視した解析は正しくない。正しくは,乗客クラスごとに調べた結果を(重み付き)平均しなければならない。

添字が3個ある($2 \times 2 \times K$ の)人数のデータを $n_{ijk}$ とする。一つの添字について合計した値を例えば $n_{i.k} = \sum_j n_{ijk}$ のように表す。このとき,

\[ \frac{(\sum_k (n_{11k} - \mu_{11k}))^2}{\sum_k \mathrm{var}(n_{11k})} \]

は自由度1のカイ2乗分布に従う。ここで

\[ \mu_{11k} = n_{1.k} n_{.1k} / n_{..k} \]

である。分散 $\mathrm{var}(n_{11k})$ は

\[ \mathrm{var}(n_{11k}) = n_{1.k} n_{2.k} n_{.1k} n_{.2k} / (n_{..k}^2 (n_{..k} - 1)) \]

で推定できる(Mantel and Haenszel, 1959)。これに先立って Cochran (1954) も同様の提案をしているが,$\mathrm{var}(n_{11k})$ の分母を $n_{..k}^3$ としている。この類の検定を (Cochran-)Mantel-Haenszel((コクラン・)マンテル・ヘンツェル)検定という。

Rには標準で mantelhaen.test() という関数が備わっている。これは $2 \times 2 \times K$ の表に対して,(デフォルトでは連続性の補正をした)Mantel-Haenszel 検定を行う:

mantelhaen.test(margin.table(Titanic, c(2,4,1)))

	Mantel-Haenszel chi-squared test with continuity correction

data:  margin.table(Titanic, c(2, 4, 1))
Mantel-Haenszel X-squared = 360.33, df = 1, p-value < 2.2e-16
alternative hypothesis: true common odds ratio is not equal to 1
95 percent confidence interval:
  8.232629 14.185153
sample estimates:
common odds ratio 
         10.80653 

このデータについては,Mantel-Haenszel で求めた common odds ratio 10.80653 は,oddsratio(margin.table(Titanic, c(2,4))) で求めた 10.14697 とあまり違わない。どんな場合にもそうであろうか,というのが次節の話題である。

ちなみに,mantelhaen.test() はオプション exact=TRUE でFisherの正確検定に対応する厳密な方法を行う。この場合,common odds ratio も conditional MLE になる。

Simpsonのパラドックス

次のようなコマンドで作られる $2 \times 2 \times 2$ の表を考えよう。

x = array(c(21, 2,71,42,
            18,70, 1,32),
          dim=c(2,2,2),
          dimnames=list(Treatment=c("Treated","Untreated"),
                        Survival=c("Alive","Dead"),
                        Sex=c("Male","Female")))

これは,男女について,ある治療をした結果を表す:

x
, , Sex = Male

           Survival
Treatment   Alive Dead
  Treated      21   71
  Untreated     2   42

, , Sex = Female

           Survival
Treatment   Alive Dead
  Treated      18    1
  Untreated    70   32

男女をまとめて1つの表にしてしまおう:

y = margin.table(x, c(1,2))
y
           Survival
Treatment   Alive Dead
  Treated      39   72
  Untreated    72   74
oddsratio(y)
[1] 0.556713

これを見ると,治療しなければ生死はほぼ等しいのに,治療したら死亡のほうが増えてしまうという,まずい治療の例のようである。Fisherの正確検定をしてみる:

fisher.test(y)

	Fisher's Exact Test for Count Data

data:  y
p-value = 0.03049
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 0.3242247 0.9524200
sample estimates:
odds ratio 
 0.5579833 

オッズ比が 1 であるという帰無仮説は棄却される($p = 0.03$)。

しかし,男女別にオッズ比を計算すれば,

apply(x, 3, oddsratio)
    Male   Female 
6.211268 8.228571 

どちらもオッズ比は 1 より大きく,治療したほうがよい。Mantel-Haenszel 検定は

mantelhaen.test(x)

	Mantel-Haenszel chi-squared test with continuity correction

data:  x
Mantel-Haenszel X-squared = 11.222, df = 1, p-value = 0.0008084
alternative hypothesis: true common odds ratio is not equal to 1
95 percent confidence interval:
  2.062906 23.283520
sample estimates:
common odds ratio 
         6.930492 

となり,オッズ比はやはり 1 より大きい($p = 0.0008$)。

このように,層別に分けた結果と,分けない結果とでは,結果が逆になることがある(Simpsonのパラドックス)。この場合は層別の結果を使うべきである。

mantelhaen.test() の中を調べると,common odds ratio 6.930492 は次のようにして求めているようである(136,121はそれぞれ男女の総数):

(21*42/136 + 18*32/121) / (71*2/136 + 1*70/121)
[1] 6.930492

メタアナリシスのように考えることもできる:

library(metafor)
y = matrix(as.numeric(x), ncol=4, byrow=TRUE)
dat = escalc(measure="OR", ai=y[,1], bi=y[,2], ci=y[,3], di=y[,4])
res = rma(yi, vi, data=dat)
forest(res, slab=c("Male","Female"))

文献