疫学では次のような2×2(2行2列)の表を考えることがよくあります:
| 疾病あり | 疾病なし | |
|---|---|---|
| 曝露あり | $a$ | $b$ |
| 曝露なし | $c$ | $d$ |
「曝露」(ばくろ,exposure)とは,何らかの条件にさらされることです(必ずしも害になるものとは限りません)。「疾病」(しっぺい,disease)は病気のことですが,より一般に「結果」(outcome)というほうがいいかもしれません。具体的には
| 肺がんあり | 肺がんなし | |
|---|---|---|
| 喫煙あり | $a$ | $b$ |
| 喫煙なし | $c$ | $d$ |
のようなことを考えればいいでしょう。
このとき,$a/(a+b)$ や $c/(c+d)$ をそれぞれの場合の危険度(risk)といい,その比 $\dfrac{a/(a+b)}{c/(c+d)}$ を相対危険度(相対リスク,relative risk,RR)といいます。
また,$a/b$ や $c/d$ をそれぞれの場合のオッズ(odds)といい,その比 $\dfrac{a/b}{c/d}$ をオッズ比(odds ratio,OR)といいます。
危険度とオッズは,小さい値のときはどちらもほぼ等しくなります。
相対危険度は,行・列のどちらを原因・結果と考えるかによって,答えが違ってきます(数学的には,行列の転置をとると,答えが違ってきます)。これに対して,オッズ比は行と列を入れ替えても変わりません。
オッズ比の対数 $\log \mathrm{OR} = \log a - \log b - \log c + \log d$ は正規分布で近似でき,その分散は,各度数をポアソン分布とすれば $V(\log a) \approx 1/a$ などが成り立つので,$V(\log \mathrm{OR}) \approx 1/a + 1/b + 1/c + 1/d$ と近似できます。相対危険度も対数をとって $V(\log \mathrm{RR}) \approx 1/a - 1/(a+b) + 1/c - 1/(c+d)$ で近似します。
オッズ比や相対危険度の信頼区間が 1 を含むことと,行・列が独立であること($a:b = c:d$)とは,数学的に同じことですが,信頼区間を求める方法や,独立性の検定の方法によって,結果が一致しないことがあります(「検定」とはそんなもので,$p = 0.04$ と $p = 0.06$ の違いは実質科学的なものではなく,単なる線引きの問題です)。
オッズ比と相対危険度には
\[ \mathrm{RR} = \frac{\mathrm{OR}(1 + c/d)}{1 + (c/d)\mathrm{OR}} = \frac{\mathrm{OR}}{1 - c/(c+d) + (c/(c+d))\mathrm{OR}} \]という関係があり,比較的安定に求められる被曝露群のオッズ $c/d$ または危険度 $c/(c+d)$ を使った式で変換できます。
これら以外に,リスク差(risk difference,RD)$a/(a+b) - c/(c+d)$ も使われます。この分散は2項分布近似で $V(\mathrm{RD}) \approx ab/(a+b)^3 + cd/(c+d)^3$ となります。
このページの最後に,ファイ係数,クラメールの $V$,Yule の $Q$ についても付記しました。
対応がある場合の2×2の表についてはMcNemar検定をご覧ください。
次のような表があったとします(この論文の表IIの男性の一部):
| 肺がんあり | 肺がんなし | |
|---|---|---|
| 喫煙あり | 231 | 23036 |
| 喫煙なし | 26 | 10813 |
この表を縦に読んでいくと,次のようになります:
c(231, 26, 23036, 10813)[1] 231 26 23036 10813
これを2行2列の行列の形にするには,次のように行数(nrow),列数(ncol)のどちらかを指定します:
matrix(c(231,26,23036,10813), nrow=2)[,1] [,2][1,] 231 23036[2,] 26 10813
数値を横に(行ごとに)読んでいった場合は,次のように byrow=TRUE というオプションを与えます:
matrix(c(231,23036,26,10813), nrow=2, byrow=TRUE)[,1] [,2][1,] 231 23036[2,] 26 10813
変数に代入しましょう:
x = matrix(c(231,26,23036,10813), nrow=2)x[,1] [,2][1,] 231 23036[2,] 26 10813
行・列に名前を付けるには次のようにします:
rownames(x) = c("喫煙","非喫煙")colnames(x) = c("肺がんあり","肺がんなし")x肺がんあり 肺がんなし喫煙 231 23036非喫煙 26 10813
手計算と照合しやすいように,数値例を簡単なものに変えました。以下では次のデータについて計算しています:
x = matrix(c(12,5,6,12), nrow=2)x[,1] [,2][1,] 12 6[2,] 5 12
オッズ比は 4.8,その対数(1.5686)の分散は 0.5333,正規分布近似の $p$ 値は 0.0317,95%信頼区間は [1.147127, 20.084960] になります:
(x[1,1]/x[1,2]) / (x[2,1]/x[2,2])[1] 4.8log((x[1,1]/x[1,2]) / (x[2,1]/x[2,2]))[1] 1.5686161/x[1,1] + 1/x[1,2] + 1/x[2,1] + 1/x[2,2][1] 0.5333333pnorm(-1.568616 / sqrt(0.5333333)) * 2[1] 0.03172043exp(1.568616 + qnorm(c(0.025,0.975)) * sqrt(0.5333333))[1] 1.147127 20.084960
関数 fisher.test() を使ってもオッズ比とその信頼区間が求められますが,上で定義したオッズ比とは少し違った定義(周辺分布を固定したときの最尤推定量,conditional MLE)を使っています。詳しくはFisherの正確検定をご覧ください。
fisher.test(x)Fisher's Exact Test for Count Datadata: xp-value = 0.04371alternative hypothesis: true odds ratio is not equal to 195 percent confidence interval:0.9465292 25.7201471sample estimates:odds ratio4.568253
オッズ比はさきほどの 4.8 ではなく,4.568253 となりました。これは定義が微妙に違うためです。詳しくはFisherの正確検定をご参照ください。また,95%信頼区間は $[0.95, 25.72]$ で,1 を含んでいるので有意でないように見えますが,$p$ 値は 0.05 未満です。
これはオッズ比などを求めるものではありませんが,行・列の独立性の検定をします。デフォルトでは連続性の補正をします。
chisq.test(x)# 補正ありPearson's Chi-squared test with Yates' continuity correctiondata: xX-squared = 3.4808, df = 1, p-value = 0.06208chisq.test(x, correct=FALSE)# 補正なしPearson's Chi-squared testdata: xX-squared = 4.8577, df = 1, p-value = 0.02752
補正あり・なしで $p$ 値はそれぞれ 0.06208,0.02752 です。
Epiパッケージで2×2の表を扱う関数は twoby2() です。いろいろ出力しますが,オッズ比は通常の定義(unconditional MLE)を使い,信頼区間はWaldの方法(正規分布による近似)で求めています。
install.packages("Epi")# 初めての場合はパッケージをインストールするlibrary(Epi)# ライブラリとして読み込むtwoby2(x)2 by 2 table analysis:------------------------------------------------------Outcome : Col 1Comparing : Row 1 vs. Row 2Col 1 Col 2 P(Col 1) 95% conf. intervalRow 1 12 6 0.6667 0.4288 0.8420Row 2 5 12 0.2941 0.1280 0.541995% conf. intervalRelative Risk: 2.2667 1.0128 5.0730Sample Odds Ratio: 4.8000 1.1471 20.0850Conditional MLE Odds Ratio: 4.5683 0.9465 25.7201Probability difference: 0.3725 0.0427 0.6073Exact P-value: 0.0437Asymptotic P-value: 0.0317------------------------------------------------------
相対リスクは 2.2667 で,95%信頼区間は $[1,01, 5.07]$ です。また,標本オッズ比は 4.8 で,95%信頼区間は $[1.15, 20.09]$ です。その下の Conditional MLE Odds Ratio は fisher.test() の出してくる値です。$p$ 値は Exact のほうは fisher.test() と同じ 0.0437 で,Asymptotic のほうは近似です。数が多い場合は Asymptotic だけの表示になり,場合によっては「0」と表示されることがありますが,その場合は「$p < 0.001$ で有意」と報告すればいいでしょう(丸める前の値を知りたければ twoby2(x)$p.value と打ち込みます)。
epitoolsパッケージには4つの方法が用意されています:
oddsratio.midp() または単に oddsratio():mid-p法(median-unbiased estimation,exact CI)oddsratio.fisher():Fisherの方法(conditional MLE,exact CI)oddsratio.wald():Waldの方法(unconditional MLE,normal approximation)oddsratio.small():正規分布近似+小標本補正(normal approximation with small sample adjustment)上のEpiパッケージの twoby2() のオッズ比と同じ結果を出すのは oddsratio.wald() です。オッズ比とその信頼区間は $measure の Exposed2 欄に出力されます:
install.packages("epitools")# 初めての場合はパッケージをインストールするlibrary(epitools)# ライブラリとして読み込むoddsratio.wald(x)$dataOutcomePredictor Disease1 Disease2 TotalExposed1 12 6 18Exposed2 5 12 17Total 17 18 35$measureodds ratio with 95% C.I.Predictor estimate lower upperExposed1 1.0 NA NAExposed2 4.8 1.147127 20.08496$p.valuetwo-sidedPredictor midp.exact fisher.exact chi.squareExposed1 NA NA NAExposed2 0.03527143 0.04371017 0.02752225$correction[1] FALSEattr(,"method")[1] "Unconditional MLE & normal approximation (Wald) CI"
ほかの方法もまとめて,違うところだけ記します:
oddsratio.wald(x)Exposed2 4.8 1.147127 20.08496oddsratio.midp(x)Exposed2 4.503795 1.105796 21.10137oddsratio.fisher(x)Exposed2 4.568253 0.9465292 25.72015oddsratio.small(x)Exposed2 3.428571 1.099686 17.37077
install.packages("vcd")library(vcd)oddsratio(x, log=FALSE)[1] 4.8
Waldの方法が使われます。単純明快にオッズ比(unconditional)だけ出力されます。信頼区間や $p$ 値は次のようにします:
confint(oddsratio(x, log=FALSE))lwr upr[1,] 1.147127 20.08496summary(oddsratio(x))Log Odds Ratio Std. Error z value Pr(>|z|)[1,] 1.5686 0.7303 2.1479 0.03172 *---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
中澤 港 先生のfmsbパッケージにもepitoolsと同じ名前の oddsratio(),rateratio(),riskratio() という関数があります。
この oddsratio() は,行列の指定と同様,縦順に数値を4つ並べます。
install.packages("fmsb")library(fmsb)oddsratio(12, 5, 6, 12)Disease Nondisease TotalExposed 12 6 18Nonexposed 5 12 17Total 17 18 35Odds ratio estimate and its significance probabilitydata: 12 5 6 12p-value = 0.0298395 percent confidence interval:1.147127 20.084959sample estimates:[1] 4.8
Waldの方法を使っています。詳細は中澤先生の鵯記878をご覧ください。
[2017-03-31追記] fmsb-0.6.0でvcdと同じ結果を返すオプションが付きました:library(fmsb); oddsratio(14, 1, 41981-14, 18327-1, p.calc.by.independence=FALSE)
exact2x2についてはFisherの正確検定にも少し書きました。このパッケージは exact2x2() という関数を定義していますが,特によく使うオプションについては fisher.exact(),blaker.exact(),mcnemar.exact() という名前でもアクセスできます。
install.packages("exact2x2")library(exact2x2)fisher.exact(x)Two-sided Fisher's Exact Test (usual method using minimum likelihood)data: xp-value = 0.04371alternative hypothesis: true odds ratio is not equal to 195 percent confidence interval:1.0905 22.9610sample estimates:odds ratio4.568253blaker.exact(x)Blaker's Exact Testdata: xp-value = 0.04371alternative hypothesis: true odds ratio is not equal to 195 percent confidence interval:1.0905 23.6488sample estimates:odds ratio4.568253
fisher.exact() の $p$ 値とオッズ比の点推定値(conditional)は fisher.test() と同じですが,信頼区間は $p$ 値と同じ方法で求めていますので,片方だけ有意ということがありません。
まとめると次のようになります:
| パッケージ | 関数 | p値 | オッズ比 | 95%信頼区間 |
|---|---|---|---|---|
| - | 手計算 | 0.03172043 | 4.8 | [1.147127, 20.084960] |
| - | fisher.test() | 0.04371 | 4.568253 | [0.9465292, 25.7201471] |
| - | chisq.test() | 0.06208 | - | - |
| - | chisq.test(correct=FALSE) | 0.02752 | - | - |
| Epi | twoby2() | Exact: 0.0437 | Sample: 4.8000 | Sample: [1.1471, 20.0850] |
| Epi | twoby2() | Asymptotic: 0.0317 | Conditional: 4.5683 | Conditional: [0.9465, 25.7201] |
| epitools | oddsratio.wald() | chisq: 0.02752225 | 4.8 | [1.147127, 20.08496] |
| epitools | oddsratio.midp() | 0.03527143 | 4.503795 | [1.105796, 21.10137] |
| epitools | oddsratio.fisher() | 0.04371017 | 4.568253 | [0.9465292, 25.72015] |
| epitools | oddsratio.small() | - | 3.428571 | [1.099686, 17.37077] |
| vcd | oddsratio() | 0.03172 | 4.8 | [1.147127, 20.08496] |
| fmsb | oddsratio() | 0.02983 | 4.8 | [1.147127, 20.084959] |
| exact2x2 | fisher.exact() | 0.04371 | 4.568253 | [1.0905, 22.9610] |
| exact2x2 | blaker.exact() | 0.04371 | 4.568253 | [1.0905, 23.6488] |
行列ではなく生データが次のように与えられている場合の方法です。logistic.display() はepiDisplay(旧epicalc)パッケージの関数です。
x = rep(0:1, c(18,17))y = rep(c(0,1,0,1), c(12,6,5,12))r = glm(y ~ x, family="binomial")logistic.display(r)Logistic regression predicting yOR(95%CI) P(Wald's test) P(LR-test)x: 1 vs 0 4.8 (1.15,20.08) 0.032 0.026Log-likelihood = -21.7558No. of observations = 35AIC value = 47.5116
津田敏秀先生の「有病オッズ比」問題に計算例があります。
相対危険度の信頼区間の求め方のいろいろについては,R help - [R] Confidence interval for relative risk という議論のMichael Deweyによるまとめ http://www.zen103156.zen.co.uk/rr.pdf や,Calculation of Relative Risk Confidence Interval - Cross Validated の議論をご覧ください。
$2 \times 2$ に限らず,表の各行(各列)の独立性を検定するためによく使われるのが「カイ2乗」($\chi^2$)という統計量です(カイ2乗検定参照)。これを0以上1以下に収まるように変換したものが「クラメールの $V$」(Cramér's $V$)です:
\[ V = \sqrt{\frac{\chi^2}{n \cdot \mathrm{min}(\mathrm{nrow}-1,\mathrm{ncol}-1)}} \]ここで $n$ は表の値を全部合計したもの,nrow と ncol は行数,列数です。特に $2 \times 2$ の表については,ファイ係数(phi coefficient,$\phi$)と呼ばれます:
\[ \phi = \sqrt{\frac{\chi^2}{n}} = \frac{|ad-bc|}{\sqrt{(a+b)(c+d)(a+c)(b+d)}} \]この分子の絶対値を外せば,どちらの向きに関連があるかもわかります。以下では絶対値なしのほうを使います。
2通りの定義でさきほどの行列 x のファイ係数を求めてみましょう:
sqrt(chisq.test(x,correct=FALSE)$statistic / sum(x))X-squared0.04054061a = x[1,1]; b = x[1,2]; c = x[2,1]; d = x[2,2](a*d-b*c) / sqrt((a+b)*(c+d)*(a+c)*(b+d))[1] 0.04054061
ライブラリ psych にファイ係数(絶対値なし)を求める関数があります:
install.packages("psych")library(psych)phi(x)[1] 0.04phi(x, digits=8)[1] 0.04054061
ファイ係数と似たものに Yule の $Q$ があります。これはオッズ比(OR)を -1 から 1 までの範囲に変換したものとも考えられます:
\[ Q = \frac{ad-bc}{ad+bc} = \frac{\mathrm{OR} - 1}{\mathrm{OR} + 1} \]これは psych パッケージの Yule() で計算できます:
Yule(x)[1] 0.6131828(a*d-b*c) / (a*d+b*c)[1] 0.6131828
これらの比較は次の論文をご覧ください:
どれを使うかは分野によると思いますが,医学・疫学方面では相対危険度(RR)やオッズ比(OR)(いずれも対数変換したもの)がよく用いられます(メタアナリシス参照)。What is the difference between odds ratio and relative risk? という議論も参考になります。
どの効果量を使うかで迷ったら,オッズ比を使えばいいでしょう。Michael Borenstein ほか Introduction to Meta-Analysis (Wiley, 2009) はオッズ比について次のように書いています:
Many people find this effect size measure less intuitive than the risk ratio, but the odds ratio has statistical properties that often make it the best choice for a meta-analysis.
$2 \times 2$ より大きい表の「効果量」はクラメールの $V$ でいいか,という話がありますが,まずは効果量の考え方に立ち戻ってデータの扱い方を考え直したほうがいいでしょう。
連続量を大小または大中小に分けて $2 \times 2$ や $2 \times 3$ の表にした場合は,確実に元の連続量に立ち戻るほうがいいでしょう。$n$ 段階の順序尺度の量は 1, 2, …, $n$ の得点に直すほうがいいでしょう。