豪華客船タイタニック号は,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
つまり,ファーストクラスでは紳士的に女性の生存が優先されたが,クラスが下がるにつれその傾向が減っていることがわかる。
個数を表す $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 になる。
次のようなコマンドで作られる $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"))