2×2の分割表について，確率の差・確率の比・オッズ比

前置き

2×2の分割表の確率の差，比，オッズ比のワルド検定およびワルド信頼区間についてまとめる．統計モデリングの勉強と接続しやすいよう，仮定している確率モデルを明示することを意識した．

モデルと点推定量

説明変数 $x_i$ を所与として，適当な関数 $g(x)$ を用い，次のような回帰型の確率モデルを考えることができる．

$$z_i \sim \mathrm{Bernoulli}(g(\alpha + \beta x_i)). \tag{M}$$

この式は確率変数 $z_i$ が成功確率 $g(\alpha + \beta x_i)$ のベルヌーイ分布に従うことを表す．確率 $g(\alpha + \beta x_i)$ で $z_i=1$（成功），$1-g(\alpha + \beta x_i)$ で $z_i=0$（失敗）の値を取るということである．また，「適当な」のところで逆関数を持つこと，微分可能であることなどを仮定している．

$(x_i, z_i)$ のペアをクロス集計すると，つまり $y_j$ ($j=0,1$) を

$$y_j = \sum_{i:x_i=j} z_i$$

とすると，次のような2×2の分割表が作れる.

	成功	失敗	横合計
$x=0$	$y_0$	$n_{0} - y_0$	$n_0$
$x=1$	$y_1$	$n_{1} - y_1$	$n_1$

このノートではこのような2×2の分割表について， (M) のモデルを通じて，仮説検定と区間推定を考えていく．

(M) のモデルでは $x_i = 0$ のときの $z_i=1$ が出る確率は

$$p_0 = g(\alpha),$$

$x_i = 1$ のときの $z_i=1$ が出る確率は

$$p_1 = g(\alpha + \beta)$$

である.

$p_i$ ($i=0,1$) の最尤推定量 $\hat p_i$ は

$$\hat p_i = y_i/n_i$$

である．最尤推定量は変数変換に対して不変なので，これを変換して推定量を導くことができる．

最初に一番簡単そうな恒等関数，すなわち $g(x)=x$ の場合について考える．このとき，

$$p_0 = \alpha,$$

$$p_1 = \alpha + \beta$$

である．最尤推定量 $\hat p_j$ を用いて，

$$\hat \alpha = \hat{p}_0 = y_{0}/n_{0},$$

を得る．ここからさらに，

$$\begin{aligned} \hat \beta &= (\hat{\alpha} + \hat{\beta}) - \hat{\alpha} =\hat{p}_1 - \hat{p}_0\\ &= y_{1}/n_{1}-y_{0}/n_{0} \tag{a} \end{aligned}$$

が得られる. $\hat \beta$ について考えることは分割表の「成功」の確率の差について考えるのと同じである．

次に， $g(x)=\exp(x)$ のとき，$p_0 = \exp(\alpha)$, $p_1 = \exp(\alpha + \beta)$ より，

$$\log p_0 = \alpha,$$

$$\log p_1 = \alpha + \beta \tag{b}$$

である．最尤推定量 $\hat p_j$ を用いて，

$$\begin{aligned} \hat \alpha &= \log \hat{p}_0 \\ &= \log(y_{0}/n_{0}), \end{aligned}$$

を得る．ここからさらに，

$$\begin{aligned} \hat \beta &= (\hat{\alpha} + \hat{\beta}) - \hat{\alpha} \\ &=\log(\hat{p}_1) - \log(\hat{p}_0)\\ &=\log(y_{1}/n_{1}) - \log(y_{0}/n_{0})\\ &=\log((y_{1}/n_{1})/(y_{0}/n_{0}))\\ \end{aligned}$$

が得られる. $\hat \beta$ について考えることは分割表の「成功」の確率の比の対数について考えるのと同じである．

最後に，$g(x) = 1/(1+\exp(-x))$ の場合を考える．このとき，

$$\begin{aligned} \log (p_0 /(1-p_0)) &= \alpha,\\ \log (p_1 /(1-p_1)) &= \alpha + \beta \end{aligned}$$

である．最尤推定量 $\hat p_j$ を用いて，

$$\begin{aligned} \hat \alpha &= \log (\hat p_0 /(1- \hat{p}_0))\\ &= \log (y_0/ (n_0- y_0)) \end{aligned}$$

および，

$$\begin{aligned} \hat \beta &= \log (\hat{p}_1 /(1-\hat{p}_1)) - \log (\hat{p}_0/(1-\hat{p}_0)) \\ &= \log (\{y_1/ (n_1- y_1)\} /\{y_0/ (n_0- y_0)\}) \tag{c} \end{aligned}$$

が得られる．

「成功」と「失敗」の比，$y_j/ (n_j- y_j)$ をオッズと呼ぶ．オッズの比 $\{y_1/ (n_1- y_1)\} /\{y_0/ (n_0- y_0)\}$ をオッズ比と呼ぶ．$\hat \beta$ について考えることは分割表のオッズ比の対数について考えるのと同じである．

まとめると次の表のようになる．

指標	モデルの $g(x)$	モデルのリンク関数
確率の差	恒等関数	恒等関数
確率の比	指数関数	対数関数
オッズ比	逆ロジット関数	ロジット関数

ここでリンク関数という言葉を導入した．一般化線形モデル（GLM; generalized linear model）の文脈では上の $g(x)$ ではなく逆関数の $g^{-1}(x)$ に注目し，これをリンク関数と呼ぶ．$g^{-1}(p_j)$ が線形になることを重視した表現である．

次に知りたいのは $\hat \beta$ の分散である．簡単な近似を使うので少し準備する．

準備：デルタ法

このパートではベルヌーイ分布に従う $Z_i$ の標本平均を次のように書く．

$$U_n = \frac{1}{n}\sum_{i=1}^n Z_i$$

いま，ベルヌーイ分布を考えていたので $U_n$ の平均と分散はそれぞれ

$$E[U_n] = \frac{1}{n} \sum_{i=1}^nE[Z_i] = \theta$$

$$V[U_n] = \frac{1}{n^2}\sum_{i=1}^n V[Z_i]= \frac{\theta(1-\theta)}{n}$$

である．中心極限定理より次が成り立つ．

$$\frac{\sqrt{n}(E[U_n]-\theta)}{\sqrt{\theta(1-\theta)}} \ \dot\sim \ \mathcal{N}(0,1)$$

ここで $\dot\sim$ は近似的に従う（左の確率変数列が $n\to \infty$ で右の分布に収束する）を意味する記号として用いた．

適当な関数の $f(x)$ について，$f(U_n)$ の平均と分散を知りたい．そこで $f(U_n)$ の1次近似，すなわち $\theta$ の周りでの1次の項までのテイラー展開，

$$\begin{aligned} f(U_n) &= f(\theta) + f'(U_n)(U_n-\theta) + \cdots \end{aligned}$$

を利用する．

両辺の平均を取ると，

$$E[f(U_n)] \approx f(\theta)$$

分散を取ると，

$$\begin{aligned} V[f(U_n)] &\approx f'(\theta) V[(U_n-\theta)]\\ &=f'(\theta) \frac{\theta(1-\theta)}{n} \end{aligned}$$

である．中心極限定理より次が成り立つ．

$$\frac{\sqrt{n}\{E[f(U_n)]-f(\theta)\}}{\sqrt{\theta(1-\theta)}} \ \dot\sim \ \mathcal{N}(0,f'(\theta))$$

この近似を利用することをデルタ法（delta method）と呼ぶ．

要するに，これを使って $f(U_n)$ の分布を次で近似しようとしている．

$$f(U_n) \ \dot\sim \ \mathcal{N}\left(f(\theta), \ \{f'(\theta)\}^2 \frac{\theta(1-\theta)}{n} \right) .$$

特に $f(x)= \log(x)$ （指数関数の逆関数）のときは，

$$\log(U_n) \ \dot\sim \ \mathcal{N}\left(\log \theta, \ \frac{1-\theta}{n \theta} \right). \tag{D1}$$

また $f(x)= \log(x/(1-x))$ （逆ロジット関数の逆関数）のとき，

$$\log(U_n) \ \dot\sim \ \mathcal{N}\left(\log (\theta/(1-\theta)), \ \frac{1}{n \theta (1-\theta)} \right) .\tag{D2}$$

推定量の平均と分散

ワルド（Wald）検定では検定統計量 $T$ が帰無仮説の下で標準正規分布に従うことを利用する（ 尤度比検定, ワルド検定, スコア検定から定まる信頼区間：二項分布の例 ）．$T$ は次式で与える．

$$T = |\hat \theta - \theta| / se(\hat \theta)$$

ここで $se(\hat \theta)$ は推定量の標準誤差であり，標準誤差は推定量の分散の平方根である．

これよりワルド検定を利用するため，推定量の分散を求めていく．

まず $g(x)$ が恒等関数のとき，(a) の分散は2項分布の分散から直接求まる．

$$Var(\hat \beta) = \sum_{j=0, 1} \frac{p_j(1-p_j)}{n_j}$$

しかし $p_j$ は不明なので $\hat p_j = y_j/n_j$ を代入すると， 推定量の分散についての推定量として次を得る．

$$Var(\hat \beta) \approx \sum_{j=0, 1} \frac{\hat p_j(1-\hat p_j)}{n_j} .$$

次に $g(x)=\exp(x)$ のとき，(b) にデルタ法 (D1) を用いると，

$$Var(\hat \beta) \approx \sum_{j=0,1} \left( \frac{1-p_j}{n_j p_j}\right).$$

である．$\hat p_j = y_j/n_j$ を代入して整理すると次を得る．

$$Var(\hat \beta) \approx \left(\sum_{j=0, 1} \frac{1}{y_j}\right)-\left(\sum_{j=0,1} \frac{1}{n_j}\right).$$

最後に $g(x)=1/(1+\exp(-x))$ のとき，(c) にデルタ法 (D2) を用いて，

$$\begin{aligned} Var(\hat \beta) &\approx \sum_{j=0,1} \left( \frac{1}{n_jp_j}+\frac{1}{n_j(1-p_j)}\right) \\ \end{aligned}$$

である． $\hat p_j = y_j/n_j$ を代入して整理し，次を得る．

$$Var(\hat \beta) \approx \left(\sum_{j=0, 1} \frac{1}{y_j} \right)+ \left(\sum_{j=0,1} \frac{1}{n_j-y_j}\right).$$

ちなみにこの $g(x)$ が逆ロジット関数の場合については Stan Lipovetsky (2015). Analytical closed-form solution for binary logit regression by categorical predictors. Journal of Applied Statistics. でも同様の議論がされている.

Rによる実装例

下記のように実装した．上で求めた式をそのまま打ち込んだ，という感じでそんなに工夫したポイントはない．

CI_wald = function(effect, se, level){
    alpha = (1-level)*0.5
    z = qnorm(1-alpha, 0, 1)
    lower = effect - se*z
    upper = effect + se*z
  return( unname(c(lower, upper)) )
}

pvalue_wald = function(effect0, effect, se){
  pvfun = fucntion(v){pchisq(((effect-v)/se)^2, df=1, lower.tail = FALSE)}
  return( sapply(effect0, pvfun) )
}

#差
RD = function(X){
  x = X[,1]
  n = rowSums(X)
  phat = x/n
  delta = unname(phat[1] - phat[2])
  se = sqrt(sum(phat*(1-phat)/n))
  return( list(estimates = delta, SE = se) )
}

#比
RR = function(X){
  x = X[,1]
  n = rowSums(X)
  phat = x/n
  beta = unname(log(phat[1]) - log(phat[2]))
  se = sqrt(sum(1/x)-sum(1/n))
  return( list(estimates=beta, SE = se) )
}

#オッズ比
OR = function(X,level=0.95, confint=TRUE){
  tau = unname(log(X[1,1]/X[1,2]) - log(X[2,1]/X[2,2]))
  se <- sqrt(sum(1/X))
  return( list(estimates=tau, SE = se) )
}

データを 佐藤俊哉『宇宙怪人しまりす統計より重要なことを学ぶ』（朝倉書店） の4話でしまりす君がやったサブグループ解析から引用させてもらう．

library(dplyr)
library(ggplot2)
library(knitr)

lname = list("ヨクナール" = c("使用","未使用"),
     "回復" = c("回復","未回復"),
     "重症度"= c("重症","軽症"))
dat = array(c(40,20,20,20,
              18,42,2,18), dim=c(2,2,2), dimnames = lname) 

# サブグループ解析（しまりす君が層別解析と間違えてやったもの）
res_rd = apply(dat, 3, RD)
res_rr = apply(dat, 3, RR)
res_or = apply(dat, 3, OR)

tab_rd = t(sapply(res_rd,function(res){c(res$estimates, CI_wald(res$estimates, res$SE, 0.95))}))
tab_rr = t(sapply(res_rr,function(res){c(res$estimates, CI_wald(res$estimates, res$SE, 0.95))}))
tab_or = t(sapply(res_or,function(res){c(res$estimates, CI_wald(res$estimates, res$SE, 0.95))}))

#差(％)
t1 = data.frame(round(tab_rd*100, 1)) %>% 
  rename("点推定値"=X1,
         "下限（95％信頼区間）"=X2,
         "上限（95％信頼区間）"=X3) %>% 
  tibble::rownames_to_column("重症度") 

print(t1)
writeLines(knitr::kable( t1 ),  "kt1.txt")

#比
t2 = data.frame(round(exp(tab_rr), 2)) %>% 
  rename("点推定値"=X1,
         "下限（95％信頼区間）"=X2,
         "上限（95％信頼区間）"=X3) %>% 
  tibble::rownames_to_column("重症度")

print(t2)
writeLines(knitr::kable( t2 ),  "kt2.txt")

#オッズ比
t3 = data.frame(round(exp(tab_or), 2)) %>% 
  rename("点推定値"=X1,
         "下限（95％信頼区間）"=X2,
         "上限（95％信頼区間）"=X3) %>% 
  tibble::rownames_to_column("重症度") 

print(t3)
writeLines(knitr::kable( t3 ),  "kt3.txt")

差（％）：

重症度	点推定値	下限（95％信頼区間）	上限（95％信頼区間）
重症	16.7	-2.9	36.2
軽症	20.0	2.5	37.5

比：

重症度	点推定値	下限（95％信頼区間）	上限（95％信頼区間）
重症	1.33	0.93	1.91
軽症	1.29	1.03	1.60

オッズ比：

重症度	点推定値	下限（95％信頼区間）	上限（95％信頼区間）
重症	2.00	0.88	4.54
軽症	3.86	0.81	18.39

テキストと同じ結果が得られた．

ついで，glm の出力とも比べてみる．

dat_c = dat[,,1]

x = c(1, 0)
fit_iden = glm(dat_c~x, family = binomial("identity"))
s_iden = summary(fit_iden)
kable(s_iden$coefficients)
kable( c(simplify2array(res_rd[[1]]), Pvalue = pvalue_wald(0, res_rd[[1]]$estimates, res_rd[[1]]$SE)) )


fit_log = glm(dat_c~x, family = binomial("log"))
s_log = summary(fit_log)
kable(s_log$coefficients)
kable( c(simplify2array(res_rr[[1]]), Pvalue = pvalue_wald(0, res_rr[[1]]$estimates, res_rr[[1]]$SE)) )

fit_logis = glm(dat_c~x, family = binomial("logit"))
s_logis = summary(fit_logis)
kable(s_logis$coefficients)
kable( c(simplify2array(res_or[[1]]), Pvalue = pvalue_wald(0, res_or[[1]]$estimates, res_or[[1]]$SE)) )

恒等リンクと確率の差：

	Estimate	Std. Error	z value	Pr(>|z|)
(Intercept)	0.5000000	0.0790569	6.324555	0.0000000
x	0.1666667	0.0997682	1.670538	0.0948129

	x
estimates	0.1666667
SE	0.0997682
Pvalue	0.0948129

対数リンクと確率の比：

	Estimate	Std. Error	z value	Pr(>|z|)
(Intercept)	-0.6931472	0.1581139	-4.383848	0.0000117
x	0.2876821	0.1825742	1.575700	0.1150950

	x
estimates	0.2876821
SE	0.1825742
Pvalue	0.1150950

ロジットリンクとオッズ比：

	Estimate	Std. Error	z value	Pr(>|z|)
(Intercept)	0.0000000	0.3162278	0.000000	1.0000000
x	0.6931472	0.4183300	1.656939	0.0975319

	x
estimates	0.6931472
SE	0.4183300
Pvalue	0.0975319

同じ結果が得られる．

あとがき

もともとやりたかったのは実装例のサブグループ解析（ 佐藤俊哉『宇宙怪人しまりす統計より重要なことを学ぶ』（朝倉書店） ）の再現で，しかもこのサブグループ解析は層別解析（層化統合解析）と間違えやすいものの例としてでてくるので，この記事の内容は本題の準備の準備くらいにあたる．

一方で，2×2の分割表についてのこれらの指標（差，比，オッズ比）自体は広く使われるものなので，この記事は『宇宙怪人しまりす統計より重要なことを学ぶ』と独立に読んでも意味のある内容になっていると思う．

また，ここでは勉強のため自作の関数を定義したが，2×2の表，オッズ比，相対危険度 - 奥村晴彦 で複数のパッケージが比較・紹介されている．