UPDATE: 2023-12-10 10:41:17

はじめに

ここでは頻度に関するデータ分析の手法についてまとめていきます。母比率\(p\)に関する話からはじめて、対応ありなしの2つの母比率\(p\)の差の区間推定と検定、2×2分割表に関するリスク、オッズや検定などをまとめていきます。

比率の信頼区間の推定

2項分布に従う確率変数の期待値と分散は\(E(X)=np, V(X)=np(1-p)\)となり、\(n\)が十分大きいとき、\(B(n,p)\)は\(N(\mu = np, \sigma^{2} = np(1-p))\)で2項分布は正規近似できる性質があります。これにより変数\(Z\)への標準化は下記のようになります。

\[ Z = \frac{\bar{X} - \mu}{\sigma} = \frac{\bar{X} - np}{\sqrt{np(1-p)}}=\frac{\frac{\bar{X}}{n} - p}{\sqrt{\frac{p(1-p)}{n}}} = \frac{\hat{p} - p}{\sqrt{\frac{p(1-p)}{n}}} \approx N(0,1) \]

また、標本比率\(\hat{p}\)の母比率、母分散は、

\[ E(\hat{p}) = p,\quad V(\hat{p}) = \frac{p(1-p)}{n},\quad SE(\hat{p}) = \sqrt{\frac{p(1-p)}{n}} \]

となります。区間推定の式の\(Z\)に先程の\(Z\)を代入して整理すると、

\[ P(|Z| \le z_{\frac{\alpha}{2}}) = 1 - \alpha \]

信頼区間は下記のようになるものの、未知のパラメータ\(p\)が含まれているので、

\[ P \left( \hat{p} - z_{\frac{\alpha}{2}}\sqrt{\frac{p(1-p)}{n}} \le p \le \hat{p} + z_{\frac{\alpha}{2}}\sqrt{\frac{p(1-p)}{n}} \right) = 1 - \alpha \]

\(p\)を\(\hat{p}\)で代用し、下記を信頼区間として利用することになります。

\[ P \left( \hat{p} - z_{\frac{\alpha}{2}}\sqrt{\frac{\hat{p}(1-\hat{p})}{n}} \le p \le \hat{p} + z_{\frac{\alpha}{2}}\sqrt{\frac{\hat{p}(1-\hat{p})}{n}} \right) = 1 - \alpha \]

サンプルサイズと信頼区間

この信頼区間の計算式を利用すれば、誤差を\(\varepsilon\)以内に抑えた上で、確率\(1-\alpha\)で標本比率を推定できます。ここで点推定値は\(\hat{p}=0.2\)だとして、誤差を\(\varepsilon=0.02=2\%\)、確率\(1-\alpha=1-0.05=95\%\)で標本比率を推定したいとします。

区間を決めているのは、下記の部分なので、

\[ Z_{\frac{\alpha}{2}}\sqrt{\frac{\hat{p}(1-\hat{p})}{n}} = \varepsilon \]

これを変形し、\(n\)について解くことで、精度をコントロールできます。

\[ n = \hat{p}(1-\hat{p})\left(\frac{Z_{\frac{\alpha}{2}}}{\varepsilon} \right)^{2} \]

実際に計算すると、標本比率\(\hat{p}=0.2=20\%\)という設定で、誤差を\(\varepsilon=0.02=2\%\)、確率\(1-\alpha=1-0.05=95\%\)で標本比率を推定したい場合、

p <- 0.20
alpha <- 0.05
eps <- 0.02
z <- qnorm(alpha/2, mean = 0, sd = 1, lower.tail = FALSE)

# サンプルサイズの計算式
n <- p*(1-p)*(z / eps)^2
sprintf('必要なサンプルサイズn = %.0f人', n)
## [1] "必要なサンプルサイズn = 1537人"

必要なサンプルサイズは1537人となります。さらに精度を上げるために、\(\varepsilon\)、確率\(1-\alpha\)を厳しくすると、サンプルサイズは自然と大きくなります。

alpha <- 0.01
eps <- 0.01
z <- qnorm(alpha/2, mean = 0, sd = 1, lower.tail = FALSE)

# サンプルサイズの計算式
n <- p*(1-p)*(z / eps)^2
sprintf('必要なサンプルサイズn = %.0f人', n)
## [1] "必要なサンプルサイズn = 10616人"

実際のところ、調査する前から標本比率\(\hat{p}\)がわかっている場合はまれなので、標本比率\(\hat{p}=0.50=50\%\) と設定して、サンプルサイズを求めることになります。標本比率\(\hat{p}=0.50=50\%\)とするのは、必要なサンプルサイズが一番大きなる標本比率\(\hat{p}\)が50%だからです。

誤差を\(\varepsilon=0.10=10\%\)、確率\(1-\alpha=1-0.05=95\%\)と固定しておいて、\(n\)が最大となる\(p\)を探してみます。

alpha <- 0.05
eps <- 0.10
z <- qnorm(alpha/2, mean = 0, sd = 1, lower.tail = FALSE)
p <- seq(0,1,0.01)
len <- length(p)
n <- vector(mode = "numeric", length = len)

for (i in 1:len) {
  n[[i]] <- p[[i]]*(1-p[[i]])*(z / eps)^2
}

# nが最大となるpのインデックスを取得して、nが最大に対応するpを求める
max_idx <- which.max(n)
max_n <- n[max_idx]
max_p <- p[max_idx]

plot(p, n, type = "l", xaxt = "n", yaxt = "n", xlab = "Probability", ylab = "Sample Size")
points(x = max_p, y = max_n, pch = 16, col = "red", cex = 1)
abline(v = max_p, lty = 2, col = "red")
abline(h = max_n, lty = 2, col = "red")
axis(side = 1, at = seq(0, 1, 0.1))
axis(side = 2, at = seq(0, 110, 5))

標本比率\(\hat{p}\)が50%としておくと、最もサンプルサイズが大きくなることがわかりました。予め 標本比率がわからない場合は、必要なサンプルサイズが大きくなることを前提に標本比率\(\hat{p}\)を50%としておくことで、必要なサンプルサイズを求めることができます。

このような考えのもと、標本誤差早見表という表が作られたりしているのだと思われます(リサーチ会社の人間ではないので実際はわからない)。この表の見方は、表側に標本比率、表頭にサンプルサイズが記載されており、任意の表側、表頭の条件のもとで、ぶつかり合うセルに誤差が記載されています。

alpha <- 0.05
z <- qnorm(alpha/2, mean = 0, sd = 1, lower.tail = FALSE)

p <- c(0.01, 0.05, 0.07, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, 0.50)
len_p <- length(p)
n <- c(50, 100, 200, 300, 400, 500, 1000, 2000, 3000, 4000, 5000)
len_n <- length(n)
eps_mat <- matrix(0, nrow = len_p, ncol = len_n)
rownames(eps_mat) <- paste0("p=", sprintf("%02s", p*100), "%-", (1-p)*100, "%")
colnames(eps_mat) <- paste0("n=",n)

for (i in 1:len_p) {
  for (j in 1:len_n) {
    eps <- z * sqrt((p[[i]] * (1-p[[i]]))/ n[[j]]) 
    eps_mat[i, j] <- round(eps * 100, 1)
  }
}

# ネットに公開されているものと少し数値がずれるのはz=2として近似しているためかもしれない
knitr::kable(eps_mat, align = "c", caption = "信頼度95%の標本誤差早見表")
信頼度95%の標本誤差早見表
n=50 n=100 n=200 n=300 n=400 n=500 n=1000 n=2000 n=3000 n=4000 n=5000
p=01%-99% 2.8 2.0 1.4 1.1 1.0 0.9 0.6 0.4 0.4 0.3 0.3
p=05%-95% 6.0 4.3 3.0 2.5 2.1 1.9 1.4 1.0 0.8 0.7 0.6
p=07%-93% 7.1 5.0 3.5 2.9 2.5 2.2 1.6 1.1 0.9 0.8 0.7
p=10%-90% 8.3 5.9 4.2 3.4 2.9 2.6 1.9 1.3 1.1 0.9 0.8
p=15%-85% 9.9 7.0 4.9 4.0 3.5 3.1 2.2 1.6 1.3 1.1 1.0
p=20%-80% 11.1 7.8 5.5 4.5 3.9 3.5 2.5 1.8 1.4 1.2 1.1
p=25%-75% 12.0 8.5 6.0 4.9 4.2 3.8 2.7 1.9 1.5 1.3 1.2
p=30%-70% 12.7 9.0 6.4 5.2 4.5 4.0 2.8 2.0 1.6 1.4 1.3
p=35%-65% 13.2 9.3 6.6 5.4 4.7 4.2 3.0 2.1 1.7 1.5 1.3
p=40%-60% 13.6 9.6 6.8 5.5 4.8 4.3 3.0 2.1 1.8 1.5 1.4
p=45%-55% 13.8 9.8 6.9 5.6 4.9 4.4 3.1 2.2 1.8 1.5 1.4
p=50%-50% 13.9 9.8 6.9 5.7 4.9 4.4 3.1 2.2 1.8 1.5 1.4

例えばサンプルサイズ\(n\)が100人で、標本比率\(\hat{p}\)が50％の場合、ぶつかり合うセルの標本誤差は9.8とあるので、標本比率の信頼区間はだいたい40%~60%の区間が母比率\(p\)を含む可能性があり、同じ条件のもとで繰り返し100回同じように標本比率を計算すれば95回は母比率\(p\)を含む区間が得られるということを意味しています。

母比率\(p\)の検定

母比率\(p\)の両側検定をする場合、下記のように仮説を設定します。

\[ \begin{eqnarray} \left\{ \begin{array}{l} H_{0}: p = p_{0} \\ H_{1}: p \neq p_{0} \end{array} \right. \end{eqnarray} \]

\(Z\)は標準正規分布\(N(0,1)\)の\(Z(\alpha/2)\)点を読み取り、その棄却域を超えていれば、有意水準\(\alpha\)で帰無仮説\(H_{0}\)を棄却します。

検定の際には、離散分布から連続分布への近似を修正する連続修正項\(- \frac{1}{2n}\)を\(Z\)が小さくならないように分子に加えます。

\[ Z = \frac{|\hat{p} - p_{0}| - \frac{1}{2n}}{\sqrt{\frac{\ p_{0}(1- p_{0})}{n}}} \]

例えば、シングルアンサーの質問に500人が回答し、60人(12%)がYesと回答したとします。このとき、他の調査の同じ質問では、標本比率が10%だったとします。今回の調査の12%は有意水準5%で検定し、比率が高いといえるのでしょうか。このような問題設定であれば、上記の式を利用して\(Z\)を計算すると、

\[ Z = \frac{|0.12 - 0.10| - \frac{1}{2 \times 500}}{\sqrt{\frac{\ 0.10(1 - 0.10)}{500}}} \]

n <- 500
p <- 60/500
p0 <- 0.10

# Zの計算式
Z <- (abs(p - p0) - (1 / (2 * n))) / sqrt((p0 * (1 - p0) / n))
sprintf('Z=%.3f', Z)
## [1] "Z=1.416"

これを\(Z(0.05/2)=1.96\)と比較することで、この棄却域を超えていれば\(H_{0}\)を棄却します。

alpha <- 0.05
Z >= qnorm(alpha/2, mean = 0, sd = 1, lower.tail = FALSE)
## [1] FALSE

結論としては、棄却域を超えないので、今回の調査は前回の調査に比べて、有意水準5%で高いとはいえない結果となります。また、p値は帰無仮説が正しいとした標準正規分布での検定統計量\(Z\)よりも外側の確率のことなので、さきほどの検定統計量\(Z\)を利用して計算できます。

# p値の計算
p <- sum(
  pnorm(q = Z, mean = 0, sd = 1, lower.tail = FALSE),
  pnorm(q = -1*Z,mean = 0, sd = 1, lower.tail = TRUE)
)
sprintf('p-value=%.4f', p)
## [1] "p-value=0.1567"

理論を理解するために手計算をしていますが、普通は便利な関数を使うことをおすすめします。ここまで計算してきた結果と同じ結果が得られています。

prop.test(x = 60, n = 500, p = 0.10, alternative = "two.sided", conf.level = 0.95)
## 
##  1-sample proportions test with continuity correction
## 
## data:  60 out of 500, null probability 0.1
## X-squared = 2.0056, df = 1, p-value = 0.1567
## alternative hypothesis: true p is not equal to 0.1
## 95 percent confidence interval:
##  0.09348364 0.15251247
## sample estimates:
##    p 
## 0.12

対応なしの2つの母比率\(p\)の差の区間推定と検定

2つの母比率\(p_{1},p_{2}\)の差を検定するためには、標本比率\(\hat{p}_{1},\hat{p}_{2}\)の差を利用することになります。母比率、母分散、母標準偏差は下記の通りです。

\[ E(\hat{p}_{1} - \hat{p}_{2}) = p_{1} - p_{2},\quad V(\hat{p}_{1} - \hat{p}_{2}) = \frac{p_{1}(1-p_{1})}{n_{1}} + \frac{p_{2}(1-p_{2})}{n_{2}},\quad SE(\hat{p}_{1} - \hat{p}_{2}) = \sqrt{\frac{p_{1}(1-p_{1})}{n_{1}} + \frac{p_{2}(1-p_{2})}{n_{2}}} \] 両側検定をする場合、下記のように仮説を設定します。

\[ \begin{eqnarray} \left\{ \begin{array}{l} H_{0}: p_{1} = p_{2} \\ H_{1}: p_{1} \neq p_{2} \end{array} \right. \end{eqnarray} \]

帰無仮説が正しいとすれば、\(p_{1}=p_{2}=p\)となるので、\(V(\hat{p}_{1} - \hat{p}_{2})\)はまとめることができます。

\[ V(\hat{p}_{1} - \hat{p}_{2}) = p(1-p)\left(\frac{1}{n_{1}} + \frac{1}{n_{2}} \right) \]

また、帰無仮説が正しいのであれば、両群のデータもまとめることができるので、

\[ \hat{p} = \frac{r_{1} + r_{2}}{n_{1} + n_{2}} \]

\[ SE(\hat{p}_{1} - \hat{p}_{2}) = \sqrt{ \hat{p}(1-\hat{p})\left(\frac{1}{n_{1}} + \frac{1}{n_{1}} \right)} \]

この結果、\(Z\)は標準正規分布\(N(0,1)\)に従うことになります。

\[ Z = \frac{(\hat{p}_{1} - \hat{p}_{2}) - E(\hat{p}_{1}-\hat{p}_{2})}{SE(\hat{p}_{1} - \hat{p}_{2})} = \frac{\hat{p}_{1} - \hat{p}_{2}}{\sqrt{\hat{p}(1-\hat{p})\left( \frac{1}{n_{1}} + \frac{1}{n_{2}} \right)}} \sim N(0, 1) \]

検定する際には、母比率の推定の時と同様に連続修正項として、\(\frac{1}{2}(\frac{1}{n_{1}} + \frac{1}{n_{2}})\)で補正します。

\[ Z = \frac{|\hat{p}_{1} - \hat{p}_{2}| - \frac{1}{2}(\frac{1}{n_{1}} + \frac{1}{n_{2}})}{\sqrt{\hat{p}(1-\hat{p})\left( \frac{1}{n_{1}} + \frac{1}{n_{2}} \right)}} \]

\(SE\)がわかっているので\(p_{1} - p_{2}\)の\(1-\alpha\)信頼区間は、

\[ \hat{p}_{1} - \hat{p}_{2} \pm Z(\alpha/2) \sqrt{\frac{\hat{p}_{1}(1-\hat{p}_{1})}{n_{1}} + \frac{\hat{p}_{2}(1-\hat{p}_{2})}{n_{2}}} \] となります。では、母比率の差の検定に移ります。下記のような2つのグループの実験データがある場合にYesと回答した比率にグループ間で差があるかどうかを両側検定します。

a <- 12; b <- 47
c <-  5; d <- 45

nm <- list(c("Group1", "Group2"), c("Yes", "No"))
mat <- matrix(c(a, b, c, d), 
              nrow = 2, byrow = TRUE,
              dimnames = nm)
mat <- addmargins(mat)
mat
##        Yes No Sum
## Group1  12 47  59
## Group2   5 45  50
## Sum     17 92 109

検定統計量\(Z\)は標準正規分布に従うので、ここでは\(Z\)の値を求め、

n1 <- mat[1,3]
p1 <- mat[1,1] / n1

n2 <- mat[2,3]
p2 <- mat[2,1] / n2

n <- mat[3,3]
p <- mat[3,1] / n

Z <- (abs(p1 - p2) - 1/2*(1/n1 + 1/n2)) / sqrt(p*(1-p)*(1/n1 + 1/n2))
Z
## [1] 1.217563

ここでは、Z値ではなくZ値に対応する\(p\)値をもとに仮説検定を行います。

# p値の計算
p <- sum(
  pnorm(q = Z, mean = 0, sd = 1, lower.tail = FALSE),
  pnorm(q = -1*Z,mean = 0, sd = 1, lower.tail = TRUE)
)
sprintf('p-value=%.4f', p)
## [1] "p-value=0.2234"

\(p\)値は0.2234となっているので、帰無仮説は棄却できず、グループ間で標本比率に差があるとはいえません。prop.test()でも同じように比率の差の検定ができるので、結果を手計算と比較してみます。

prop.test(x = c(mat[1,1], mat[2,1]), n = c(n1, n2), alternative = "two.sided", conf.level = 0.99, correct = TRUE)
## 
##  2-sample test for equality of proportions with continuity correction
## 
## data:  c(mat[1, 1], mat[2, 1]) out of c(n1, n2)
## X-squared = 1.4825, df = 1, p-value = 0.2234
## alternative hypothesis: two.sided
## 99 percent confidence interval:
##  -0.08876026  0.29553992
## sample estimates:
##    prop 1    prop 2 
## 0.2033898 0.1000000

prop.test()は\(\chi^{2}\)値を検定統計量として利用しているので、\(Z\)の値とは比較できませんが、\(Z^2\)は自由度1の\(\chi^2\)分布に従うので、\(Z^2=1.217563^2=1.48246\)となって一致します。また、\(p\)値も一致していることがわかります。検定結果も同じで、グループ間で標本比率に差があるとはいえません。

95％信頼区間も計算しておきます。

z <- qnorm(0.05/2, mean = 0, sd = 1, lower.tail = FALSE)
upr <- (p1 - p2) + 1.96 * sqrt((p1*(1-p1) / n1) + (p2*(1-p2) / n2))
lwr <- (p1 - p2) - 1.96 * sqrt((p1*(1-p1) / n1) + (p2*(1-p2) / n2))
c(sprintf('lower=%.4f', lwr), sprintf('upper=%.4f', upr))
## [1] "lower=-0.0288" "upper=0.2355"

この式からも分かる通り、連続修正項を含んでいないので、prop.test(correct = FALSE)とすることで、同じ結果が得られます。

prop.test(x = c(mat[1,1], mat[2,1]), n = c(n1, n2), alternative = "two.sided", conf.level = 0.95, correct = FALSE)$conf.int
## [1] -0.02876091  0.23554057
## attr(,"conf.level")
## [1] 0.95

対応ありの2つの母比率\(p\)の差の区間推定と検定(McNemar検定)

さきほどの2つの母比率\(p\)の差の区間推定と検定では、回答者は「対応なし」という関係を想定していました。ここでは、「対応あり」の場合の2つの母比率\(p\)の差の区間推定と検定を行います。下記のデータでは、1回目も2回目も同じ54人が回答しており、1回目と2回目の回答の違いをテーブルとしてまとめています。

a <- 25; b <- 7
c <-  1; d <- 21

nm <- list(c("Yes1", "No1"), c("Yes2", "No2"))
mat_ <- matrix(c(a, b, c, d), 
              nrow = 2, byrow = TRUE,
              dimnames = nm)
mat <- addmargins(mat_)
mat
##      Yes2 No2 Sum
## Yes1   25   7  32
## No1     1  21  22
## Sum    26  28  54

ここで興味がある標本比率は、\(\hat{p}_{1} = \frac{a + b}{N}, \hat{p}_{2} = \frac{a + c}{N}\)となります。1回目に Yesと回答した比率、2回目にYesと回答した比率の差が興味の対象です。セルの要素は\(a: [1,1]\quad b: [1,2]\quad c: [2,1]\quad d: [2,2]\)です。

\[ \hat{p}_{1} - \hat{p}_{2} = \frac{a + b - (a + c)}{N} = \frac{b - c}{N} \]

これまでと同じように両側検定の仮説を設定します。

\[ \begin{eqnarray} \left\{ \begin{array}{l} H_{0}: p_{1} = p_{2} \\ H_{1}: p_{1} \neq p_{2} \end{array} \right. \end{eqnarray} \]

対応のある比率の差の検定では、帰無仮説のもとで、標準誤差は下記の通りとなり、

\[ SE(\hat{p}_{1} - \hat{p}_{2}) = \frac{\sqrt{b+c}}{N} \]

連続修正項を加えて、近似的に、検定統計量\(Z\)は標準正規分布\(N(0,1)\)に従います。

\[ Z = \frac{(\hat{p}_{1} - \hat{p}_{2}) \pm \frac{1}{N}}{SE(p_{1} - p_{2})} = \frac{(b-c)/N \pm 1/N}{(b-c/N)} = \frac{(b-c)\pm 1}{\sqrt{b+c}} \sim N(0,1) \]

\(Z\)の値を\(Z^2\)とすることで、これは自由度1の\(\chi^2\)分布に従うので、帰無仮説の棄却域を\(\chi_{1}^{2}\)分布で構築します。

\[ Z^{2}=\chi_{1}^{2} = \left(\frac{(b-c)\pm 1}{\sqrt{b+c}}\right)^{2} = \frac{(|b-c|- 1)^2}{b+c} \]

信頼区間は近似的に下記で構成されます。

\[ \hat{p}_{1} - \hat{p}_{2} \pm Z(\alpha/2)SE(\hat{p}_{1} - \hat{p}_{2}) = \frac{b-c}{N} \pm \frac{Z(\alpha/2)}{N}\sqrt{b+c-\frac{(b-c)^2}{N}} \] 今回も手計算で\(\chi_{1}^{2}\)を求めていきます。

\[ Z^{2} = \chi_{1}^{2} = \frac{(|b-c|- 1)^2}{b+c} = \frac{(|7-1|-1)^2}{7+1} = 3.125 \]

この値を利用して、\(\chi_{1}^{2}\)分布から\(p\)値を計算します。

p <- pchisq(q = 3.125, df = 1, lower.tail = FALSE)
sprintf('p-value=%.4f', p)
## [1] "p-value=0.0771"

\(\chi_{1}^{2}\)、\(p\)値のいずれも、mcnemar.test()の結果と一致しました。今回の検定統計量では、5％有意水準の棄却域を超えず、帰無仮説は棄却できず、標本比率の差に違いはないということになります。

# matではなくmat_を使用する
mcnemar.test(x = mat_, correct = TRUE)
## 
##  McNemar's Chi-squared test with continuity correction
## 
## data:  mat_
## McNemar's chi-squared = 3.125, df = 1, p-value = 0.0771

mcnemar.test()は信頼区間を出力しないので、計算結果の確認はできませんが、手計算で信頼区間を計算しておきます。

\[ \frac{b-c}{N} \pm \frac{Z(\alpha/2)}{N}\sqrt{b+c-\frac{(b-c)^2}{N}} = \frac{7-1}{54} \pm \frac{1.96}{54}\sqrt{7+1-\frac{(7-1)^2}{54}} = 0.0128 \sim 0.2094 \]

#a <- 25; b <- 7
#c <-  1; d <- 21
z <- qnorm(0.05/2, mean = 0, sd = 1, lower.tail = FALSE)
upr <- (b-c)/sum(mat_) + z / sum(mat_) * sqrt(b + c - ((b-c)^2/sum(mat_)))
lwr <- (b-c)/sum(mat_) - z / sum(mat_) * sqrt(b + c - ((b-c)^2/sum(mat_)))
c(sprintf('lower=%.4f', lwr), sprintf('upper=%.4f', upr))
## [1] "lower=0.0128" "upper=0.2094"

\(\chi^{2}\)検定

帰無仮説で与えているモデルと観測データの適合度合いを検定したい場合に適合度の検定が利用できます。適合度の検定やこの後に紹介する独立性の検定は、検定統計量に\(\chi^{2}\)値を利用することからまとめて\(\chi^{2}\)検定と呼ばれることもあります。

適合度の検定では、ベルヌイ試行を一般化した多項試行を考えます。つまり、\(k\)個の互いに排反なカテゴリのいずれかの1つに属し、1回の試行は独立しており、\(i(i=1 \sim k)\)番目のカテゴリに入る確率は試行を通じて\(p_{i}\)とします。\(\sum_{i=1}^{k} p_{i} = 1\)です。

このような試行を\(n\)回行い、確率変数である各カテゴリの度数\(n_{i}\)が得られたとします。\(\sum_{i=1}^{k} n_{i} = n\)で、自由度は\(n\)という制約より\(k-1\)です。

多項試行ではあるものの\(n\)回の試行結果で、カテゴリ\(i\)に入るか、入らないかに注目すれば良いから、\(E(n_{i})=np_{i}, \quad V(n_{i})=np_{i}q_{i}\)の二項分布に従う確率変数とも考えられ、\(np_{i}\)は試行を\(n\)回行うときの、確率\(p_{i}\)をもつセル\(i\)の期待度数となります。

帰無仮説は、カテゴリ確率\(p_{i}\)が帰無仮説\(H_{0}\)によって\(p_{i0}\)として特定される場合を考えます。

\[ \begin{eqnarray} \left\{ \begin{array}{l} H_{0}: p_{1} = p_{10}, \dots, p_{k} = p_{k0}\\ H_{1}: 少なくとも1つのp_{i} \neq p_{i0} \end{array} \right. \end{eqnarray} \]

\(n\)が十分大きい時、近似的に\(\chi^2\)値を用いて実施でき、もし\(H_{0}\)が正しければ、期待度数\(np_{i0}\)と観測度数\(n_{i}\)との差は小さくなるように検定統計量\(\chi^2\)値はつくられています。また、観測度数を\(O_{i}\)とし、期待度数を\(E_{i}\)として表記されていることもあります。

\[ \chi^2_{k-1} \sim \sum_{i=1}^{k} \frac{(n_{i} - np_{i0})}{np_{i0}} = \sum_{i=1}^{k} \frac{(O_{i} - E_{i})}{E_{i}} \]

このような検定統計量をとるので、特段何もなければ、基本的には右側の棄却域のみを設定することが多いです。なぜこれが\(\chi^2\)分布に従うのかは、多項分布、多次元正規分布、多次元確率変数変換などを使って証明されます。それはさておき、実際に検定していきましょう。

特定の質問にYesと回答する人に対して、25人が男性、5人が女性の時、男女で偏りがあるかどうかを調べたいとします。期待度数は\(E=15(=30/2)\)となり観測度数となります。

o <- c(25, 5)
e <- sum(o)/length(o)
chi2 <- sum(
  ((o - e) ^ 2 / e)
  )

# 有意水準は5%
chi2 > qchisq(0.95, df = length(o)-1)
## [1] TRUE

\(\chi^2\)値を計算すると、有意水準5%の棄却域は3.84となり、\(\chi^2\)値は13.33より帰無仮説は棄却されます。手計算する必要はなく、chisq.test()を利用しても同じ結果が得られます。

chisq.test(o)
## 
##  Chi-squared test for given probabilities
## 
## data:  o
## X-squared = 13.333, df = 1, p-value = 0.0002607

独立性の検定を見ていきましょう。独立性の検定は適合度の検定から異なる2つのカテゴリを用いた適合度の検定を拡張したものなので、簡単に紹介しておきます。下記のようなk×lのクロス集計表が利用されます。

a <- 29; b <- 80
c <- 10; d <- 94
mat_ <- matrix(c(a, b, c, d), nrow = 2, byrow = TRUE,
               dimnames = list(c("Group1","Group2"), c("Yes", "No")))
mat <- addmargins(mat_)
mat
##        Yes  No Sum
## Group1  29  80 109
## Group2  10  94 104
## Sum     39 174 213

各\((i,j)\)要素の観測度数\(O_{ij}\)、横計\(n_{i}\)、縦計\(m_{j}\)、全度数\(N\)とします。帰無仮説は、2つのカテゴリが独立、または、分布が一様であるということになります。\(H_{0}\)のもとで、期待度数は\(E_{ij} = \frac{n_{i}m_{j}}{N}\)となります。

\[ \chi^2_{(k-1)(l-1)}{} \sim \sum_{i=1}^{k}\sum_{j=1}^{l} \frac{(O_{ij} - E_{ij})}{E_{ij}} \]

適合度の検定と同様にchisq.test()を利用して結果が得られます。

chisq.test(mat)
## 
##  Pearson's Chi-squared test
## 
## data:  mat
## X-squared = 10.271, df = 4, p-value = 0.0361

リスクとオッズに関する推定

リスク

クロス集計表の分析方法として、リスクとオッズを使った方法もあるので、その内容もまとめておきます。下記のようなクロス集計表があった場合に、

group yes no n
\(group1\) \(a(=r_{1})\) \(b\) \(n_{1}\)
\(group2\) \(c(=r_{2})\) \(d\) \(n_{2}\)

リスク差(Risk Difference)、リスク比(Risk Ratio)は下記の通り定義されます。

\[ \begin{eqnarray} \hat{RD} &=& \frac{\hat{p_{1}}}{\hat{p_{2}}} = \frac{\frac{r_{1}}{n_{1}}}{\frac{r_{2}}{n_{2}}} \\ \hat{RR} &=& \hat{p_{1}} - \hat{p_{2}}= \frac{r_{1}}{n_{1}}-\frac{r_{2}}{n_{2}} \end{eqnarray} \]

Epiパッケージは疫学研究に必要なものがまとまっているパッケージなので、テーブルを渡すだけで豊富なアウトプットを提供してる便利なパッケージです。

library(Epi)
a <- 29; b <- 80
c <- 10; d <- 94
mat_ <- matrix(c(a, b, c, d), nrow = 2, byrow = TRUE,
               dimnames = list(c("Group1","Group2"), c("Yes", "No")))
mat <- addmargins(mat_)
mat
##        Yes  No Sum
## Group1  29  80 109
## Group2  10  94 104
## Sum     39 174 213

ここでは、このパッケージを利用して、リスク比に注目していきます。

res_risk <- Epi::twoby2(mat_)
## 2 by 2 table analysis: 
## ------------------------------------------------------ 
## Outcome   : Yes 
## Comparing : Group1 vs. Group2 
## 
##        Yes No    P(Yes) 95% conf. interval
## Group1  29 80    0.2661    0.1916   0.3567
## Group2  10 94    0.0962    0.0525   0.1696
## 
##                                    95% conf. interval
##              Relative Risk: 2.7670    1.4206   5.3893
##          Sample Odds Ratio: 3.4075    1.5649   7.4195
## Conditional MLE Odds Ratio: 3.3885    1.4937   8.2938
##     Probability difference: 0.1699    0.0668   0.2696
## 
##              Exact P-value: 0.0014 
##         Asymptotic P-value: 0.0020 
## ------------------------------------------------------

結果を見ると、グループ1のリスクは0.27で、グループ2のリスクは0.1となっており、この比をとったリスク比は2.77です。

この結果から、グループ2よりもグループ1のほうが2.77倍起こりやすく、95％信頼区間はlower=0.0128, upper=0.2094であることがわかりました。

信頼区間、検定ではリスク比の対数が正規分布に近似できる性質を利用して、近似的に下記の通り定義されます。まず標準誤差は、

\[ SE \left( log\left[\frac{\hat{p}_{1}}{\hat{p}_{2}}\right] \right) = \sqrt{\frac{1-\hat{p}_{1}}{n_{1}\hat{p}_{1}} + \frac{1-\hat{p}_{2}}{n_{2}\hat{p}_{2}}} = \sqrt{\frac{1}{r_{1}} - \frac{1}{n_{1}} + \frac{1}{r_{2}} - \frac{1}{n_{2}}} \]

であり、この標準誤差を利用し、信頼区間は下記の通り構成されます。

\[ \exp \left( log \left[\frac{\hat{p}_{1}}{\hat{p}_{2}} \right] \pm Z(\alpha/2) \sqrt{\frac{1-\hat{p}_{1}}{n_{1}\hat{p}_{1}} + \frac{1-\hat{p}_{2}}{n_{2}\hat{p}_{2}}}\right) = \exp \left( log \left[\frac{\hat{p}_{1}}{\hat{p}_{2}} \right] \pm Z(\alpha/2) \sqrt{\frac{1}{r_{1}} - \frac{1}{n_{1}} + \frac{1}{r_{2}} - \frac{1}{n_{2}}}\right) \] Epi::twoby2()のリスク比の信頼区間はlower=1.4206, upper=5.3893でした。実際に検算してみます。

rr <- res_risk$measures[1,1]
z <- qnorm(0.05/2, mean = 0, sd = 1, lower.tail = FALSE)
se <- z * sqrt(1/mat[1,1] - 1/mat[1,3] + 1/mat[2,1] - 1/mat[2,3])

upr <- exp(log(rr) + se)
lwr <- exp(log(rr) + -1*se)

c(sprintf('lower=%.4f', lwr), sprintf('upper=%.4f', upr))
## [1] "lower=1.4206" "upper=5.3893"

だいたい同じ値が得られています。次に、検定統計量を計算します。検定統計量は比率の差の検定と同じ計算式で計算できます。帰無仮説はリスク比が0(\(RR=0\))です。

\[ Z = \frac{|\hat{p}_{1} - \hat{p}_{2}| - \frac{1}{2}(\frac{1}{n_{1}} + \frac{1}{n_{2}})}{\sqrt{\hat{p}(1-\hat{p})\left( \frac{1}{n_{1}} + \frac{1}{n_{2}} \right)}} \]

risk1 <- mat[1,1]/mat[1,3]; n1 <- mat[1,3]
risk2 <- mat[2,1]/mat[2,3]; n2 <- mat[2,3]
risk  <- mat[3,1]/mat[3,3]; n <- mat[3,3]

z <-  (abs(risk1-risk2) - 1/2*(1/n1 + 1/n2)) / sqrt(risk*(1-risk)*(1/n1 + 1/n2))
z
## [1] 3.027645

統計量はz=3.03より、5％有意水準のz=1.96よりも大きいので、帰無仮説は棄却されます。

オッズ

リスクとオッズは分析データの取得状態によって使い分けられる事が多く、リスクは前向き(コホート)研究(要因→結果)、オッズは後ろ向き(ケースコントロール)研究(結果→要因)で使用されます。一般に、後ろ向き研究の場合、リスク比、リスク差はリスクの発生指標としてあまり良い推定が得られないことが知られています。そのため、後ろ向き研究ではオッズを利用します。下記の動画で詳しく解説されています。

オッズは当たる確率と、当たらない確率の比で定義され、1よりも大きければ、当たりやすいことを表します。例えば、グループ1のオッズは下記の通りです。

group yes no n
\(group1\) \(a(=r_{1})\) \(b\) \(n_{1}\)
\(group2\) \(c(=r_{2})\) \(d\) \(n_{2}\)

\[ odds = \frac{p}{1-p} = \frac{\frac{a}{a+b}}{\frac{b}{a+b}} = \frac{a}{b} \] 同じように、グループ2のオッズは\(\frac{c}{d}\)となるので、オッズの比をとると、

\[ Odds\ Ratio = \frac{\frac{a}{b}}{\frac{c}{d}} = \frac{ad}{bc} \]

ただ、リスクとは異なり「Hogeリスクがx倍起こりやすい」とは言えず、「オッズ比がx倍である」という言い方、解釈になります。

オッズ比はEpi::twoby2()の出力結果に含まれおり、オッズ比の信頼区間はlower=1.5649, upper=7.4195です。

res_risk
## $table
##        Yes No     P(Yes)  95% conf.  interval
## Group1  29 80 0.26605505 0.19161304 0.3566581
## Group2  10 94 0.09615385 0.05251862 0.1695556
## 
## $measures
##                                        95% conf.  interval
##              Relative Risk: 2.7669725 1.42062629 5.3892686
##          Sample Odds Ratio: 3.4075000 1.56493364 7.4195199
## Conditional MLE Odds Ratio: 3.3884770 1.49367209 8.2937671
##     Probability difference: 0.1699012 0.06678624 0.2695782
## 
## $p.value
## [1] 0.002015153 0.001397784

オッズ比の標準誤差は、

\[ SE \left( log\left[\hat{OR}\right] \right) = \sqrt{\frac{1}{a} + \frac{1}{b} + \frac{1}{c} + \frac{1}{d}} \]

であり、この標準誤差を利用し、オッズ比の信頼区間は下記の通り構成されます。オッズ比の対数が正規分布に近似できる性質が利用されています。

\[ \exp \left( log \left[\hat{OR} \right] \pm Z(\alpha/2) \sqrt{\frac{1}{a} + \frac{1}{b} + \frac{1}{c} + \frac{1}{d}}\right) \]

オッズ比の検定は下記の通りで、\(\hat{OR}=\frac{ad}{bc} \ge 1\)のときマイナスを取り、\(\hat{OR}=\frac{ad}{bc} \lt 1\)のときプラスを取ります。この統計量は超幾何分布が下記の式で近似できることを利用しています。

\[ \chi = \frac{\sqrt{N-1}{(ad-bc) \pm \frac{N}{2}}}{\sqrt{n_{1}n_{2}m_{1}m_{2}}} \sim N(0,1) \]

リスクとオッズの補足

基本的には前向き研究ではリスク比、後ろ向きではオッズ比が使われる。前向き研究では、統制、介入群のサンプルサイズをコントロールすることが可能である一方、後ろ向き研究は統制介入ではなく、結果からサンプルを決定することになる。

前向き研究では、この薬を投与することで病気が寛解するかどうか、という時間の流れである一方、後ろ向き研究では、病気が寛解した(してない人)で、薬が投与されていたかどうか、という時間の流れを遡ることになる。この点でリスク比は結果が一貫しない。

下記のテーブルのリスク比は2.7、オッズ比は3.4である。

res_risk
## $table
##        Yes No     P(Yes)  95% conf.  interval
## Group1  29 80 0.26605505 0.19161304 0.3566581
## Group2  10 94 0.09615385 0.05251862 0.1695556
## 
## $measures
##                                        95% conf.  interval
##              Relative Risk: 2.7669725 1.42062629 5.3892686
##          Sample Odds Ratio: 3.4075000 1.56493364 7.4195199
## Conditional MLE Odds Ratio: 3.3884770 1.49367209 8.2937671
##     Probability difference: 0.1699012 0.06678624 0.2695782
## 
## $p.value
## [1] 0.002015153 0.001397784

この調査に対し、病気が寛解した人を都合よく5倍増やしたとすると、下記のテーブルのリスク比は1.8、オッズ比は3.4のまま変化しない。

mat_[,1] <- mat_[,1]*5
Epi::twoby2(mat_)
## 2 by 2 table analysis: 
## ------------------------------------------------------ 
## Outcome   : Yes 
## Comparing : Group1 vs. Group2 
## 
##        Yes No    P(Yes) 95% conf. interval
## Group1 145 80    0.6444    0.5797   0.7043
## Group2  50 94    0.3472    0.2740   0.4284
## 
##                                    95% conf. interval
##              Relative Risk: 1.8560    1.4540   2.3691
##          Sample Odds Ratio: 3.4075    2.1980   5.2825
## Conditional MLE Odds Ratio: 3.3955    2.1487   5.4158
##     Probability difference: 0.2972    0.1938   0.3914
## 
##              Exact P-value: 0.0000 
##         Asymptotic P-value: 0.0000 
## ------------------------------------------------------

このように、サンプルサイズのコントロールの仕方によって、オッズ比は変化しないが、リスク比は変化してしまう。そのため、後ろ向き研究では結果からサンプルサイズが変化するため、オッズ比が優先的に使われることになる。ただ、オッズ比はリスク比のようには解釈できないのため注意が必要。