オッズとオッズ比
2つの割合と
に対して、オッズ
の比
を考えると、このオッズ比が1より大なら
、1より小なら
である。
オッズ比は、ロジスティック回帰において回帰係数がオッズ比の対数に一致することや、患者対照研究のようにオッズ比のみが推定できる研究方法があるためにしばしば用いられる。が小さい場合には、オッズ比は割合の比に近似するので、リスクが何倍になる、といった文脈で使われる場合もある。
独立性のカイ二乗検定
2つの二項割合、
について、帰無仮説を
とする。このとき、割合の推定量の差を分散のルートで補正した量の2乗をピアソンのカイ二乗統計量
と呼び、データ数が大きい時に
が自由度1のカイ二乗分布で近似されるため、この近似によって棄却域やp値が計算できる。
このとき、より近似を良くするためにイェーツの補正と呼ばれる補正を行う場合がある。Rのchisq.test()ではデフォルトである。
> Ctable <- matrix(c(512, 1385, 567, 2419), nr = 2) > chisq.test(Ctable) Pearson's Chi-squared test with Yates' continuity correction data: Ctable X-squared = 42.679, df = 1, p-value = 6.449e-11
この補正を行ってもデータ数が少ないと近似が不正確になり、警告メッセージが表示される場合がある。
> Ctable <- matrix(c(2, 10, 3, 15), nr = 2) > chisq.test(Ctable) Pearson's Chi-squared test data: Ctable X-squared = 0, df = 1, p-value = 1 警告メッセージ: chisq.test(Ctable) で: カイ自乗近似は不正確かもしれません
フィッシャーの直接確率法
紅茶が先かミルクが先かという話は、「ミルクが先」が美味しいということで決着が付いているらしい*1が、ともかく、サンプル数が少ない場合に超幾何分布モデルによって結果が偶然に起こる確率を直接計算する方法をフィッシャーの直接確率法と呼ぶ。
Rではfisher.test()で行う。なお、信頼区間や推定値はオッズ比として表示される。
> table <- matrix(c(3, 1, 1, 3), nr = 2) > fisher.test(table) Fisher's Exact Test for Count Data data: table p-value = 0.4857 alternative hypothesis: true odds ratio is not equal to 1 95 percent confidence interval: 0.2117329 621.9337505 sample estimates: odds ratio 6.408309
かつてはデータ数が大きくなると計算量が増えるという点が問題であったが、2x2程度の表であれば多少データ数が大きくても実行に支障はない。
> Ctable <- matrix(c(512, 1385, 567, 2419), nr = 2) > fisher.test(Ctable) Fisher's Exact Test for Count Data data: Ctable p-value = 8.63e-11 alternative hypothesis: true odds ratio is not equal to 1 95 percent confidence interval: 1.372396 1.811955 sample estimates: odds ratio 1.576932
関連性の指標
2x2表が独立ではない場合にどの程度の関連性があるかを調べる指標をまとめて計算できる関数assocstats()がパッケージvcdに入っている。また、oddsratio()でオッズ比も計算できる。
> library(MASS) > library(vcd) > Ctable <- xtabs(~ Sex + W.Hnd, data = survey) > assocstats(Ctable) X^2 df P(> X^2) Likelihood Ratio 0.54629 1 0.45984 Pearson 0.54351 1 0.46098 Phi-Coefficient : 0.048 Contingency Coeff.: 0.048 Cramer's V : 0.048 > oddsratio(Ctable) log odds ratios for Sex and W.Hnd [1] -0.3750241 > oddsratio(Ctable, log = FALSE) # デフォルトでは対数表示 odds ratios for Sex and W.Hnd [1] 0.6872727
2xJ表の解析
カテゴリカル変数が2つよりも大きい場合もchisq.test()は実行できる。
> Ctable <- xtabs(~ Treatment + Improved, data = Arthritis) > Ctable Improved Treatment None Some Marked Placebo 29 7 7 Treated 13 7 21 > chisq.test(Ctable) Pearson's Chi-squared test data: Ctable X-squared = 13.055, df = 2, p-value = 0.00146
ただし、このときは単にPlaceboとTreatedの分布が違う、程度のことしか言えない。どの値が帰無仮説からのずれが大きいかを調べる方法として、調整済み残差を用いる方法がある。
> Ctable <- xtabs(~ Treatment + Improved, data = Arthritis) > result <- chisq.test(Ctable) > v <- result$expected * (1 - rowSums(Ctable)/sum(Ctable)) %o% + (1 - colSums(Ctable)/sum(Ctable)) > (Ctable - result$expected) / sqrt(v) Improved Treatment None Some Marked Placebo 3.27419655 -0.09761768 -3.39563632 Treated -3.27419655 0.09761768 3.39563632
この残差が1.96より大きければ帰無仮説からのずれが大きいと考える。
カテゴリカル変数に順序関係のある場合の解析
マンテル検定
カテゴリカル変数間に順序関係が有るとき、これに1, 2, 3...のようなスコアを与える検定手法としてマンテル検定がある。
manteltrend.test <- function(x, score = c(1:dim(x)[2])){ ctotal <- colSums(x) # 列和 rtotal <- rowSums(x) # 行和 n <- sum(x) # データ総数 expected <- (rtotal %o% ctotal) / n # 期待度数 tscore <- (x - expected) %*% score v <- n * (score^2 %*% ctotal) - (score %*% ctotal)^2 statistics <- n^2 * (n-1) * tscore[2]^2 / rtotal[1] / rtotal[2] / v p.value <- pchisq(statistics, 1, lower.tail = FALSE) list(statistics = statistics, p_value = p.value) }
スコアは等間隔である必要は無いが、どのようなスコア間隔を用いるかは大勢に影響を与えない場合が多い。
> manteltrend.test(Ctable, score=0:2) $statistics [,1] [1,] 12.85902 $p_value [,1] [1,] 0.0003358568 > manteltrend.test(Ctable, score=c(1, 2, 5)) $statistics [,1] [1,] 12.60627 $p_value [,1] [1,] 0.0003844559
Wilcoxonの順位和検定
Wilcoxonの順位和検定では、スコアを使わず、各データを順位に変換して検定を行う。Rではwilcox.test()で行えるが、あらかじめデータを1人1人のデータに戻しておく必要がある。
> x <- rep(col(Ctable)[1,], Ctable[1,]) > y <- rep(col(Ctable)[2,], Ctable[2,]) > wilcox.test(x, y) Wilcoxon rank sum test with continuity correction data: x and y W = 517.5, p-value = 0.0003666 alternative hypothesis: true location shift is not equal to 0 警告メッセージ: wilcox.test.default(x, y) で: タイがあるため、正確な p 値を計算することができません
コクラン・アーミテージ検定
薬の投与量に応じて効果のありなしが分かれるようなデータがあったとして、効果ありの割合が投与量
に対して、次のような関係にあったとする。
ここでであるとすると、投与量と効果ありの割合の間に関係がないということになる。このとき、ロジスティック回帰モデルのように
を何らかの関数で変換したとしても、検定方法は変わらない。この検定はRではprop.trend.test()で実行できる。
> gd <- c(2, 3, 9, 15, 18) > n <- rep(20, 5) > prop.trend.test(gd, n, score = 1:5) Chi-squared Test for Trend in Proportions data: gd out of n , using scores: 1 2 3 4 5 X-squared = 38.86, df = 1, p-value = 4.553e-10
I x J表の解析
独立性の検定
chisq.test()がこれまでと同様に使用できる。
> Ctable <- xtabs(Freq ~ Status + Alcohol, data = DanishWelfare) > chisq.test(Ctable) Pearson's Chi-squared test data: Ctable X-squared = 87.972, df = 4, p-value < 2.2e-16
線形関連の検定
行と列が共に順序カテゴリカルデータの場合、これらにスコアを与えて相関係数を求めることができる。相関係数を、標本数を
としたとき、
は、標本数が大きくなると近似的に自由度1のカイ二乗分布に近づく。このことを利用した検定を線形関連の検定と呼ぶ。
Rではcoinパッケージのlbl_test()関数で実行できる。
> library(coin) > Ctable <- apply(jobsatisfaction, 1:2, sum) # 性別を合計 > Ctable Job.Satisfaction Income Very Dissatisfied A Little Satisfied Moderately Satisfied <5000 2 4 13 5000-15000 2 6 22 15000-25000 0 1 15 >25000 0 3 13 Job.Satisfaction Income Very Satisfied <5000 3 5000-15000 4 15000-25000 8 >25000 8 > lbl_test(as.table(Ctable)) Asymptotic Linear-by-Linear Association Test data: Job.Satisfaction (ordered) by Income (<5000 < 5000-15000 < 15000-25000 < >25000) Z = 2.7623, p-value = 0.005739 alternative hypothesis: two.sided
スコアは指定しなければ自動的に1, 2, 3, ...と設定されるが、指定したければlistで与える。
lbl_test(as.table(Ctable), scores = list(Job.Satisfaction = 1:4, Income = c(2.5, 10, 20, 30)))
なお、I x J表においてI=2の場合はマンテル検定と同じ結果が得られる。
