尤度比検定, ワルド検定, スコア検定から定まる信頼区間：二項分布の例

導入

数理統計の本には良く「信頼区間は検定方式の反転によって得られる」というようなことが書かれている.

と, 言うときにいま僕が手元で参照しているのは『現代数理統計学の基礎』（久保川達也, 共立出版）だが, これ以外にも表紙のどこかに「数理統計学」と入っているような本には同様のことが書かれていることが多いと思う.

両側検定,

帰無仮説: $\theta = \theta_0$（1点）
対立仮説: $\theta \neq \theta_0$

の検定統計量 $T(X)$（標本$X$の関数）と棄却域$R(\theta_0)$（帰無仮説のパラメータに依存）から,

$$P_{\theta_0}(T(X) \in R(\theta_0)) = \alpha$$

とできるなら,

$$P_{\theta_0}(T(X) \in R^c(\theta_0)) = 1-\alpha$$

なので（$R^c(\theta_0)$ は受容域）これを $\theta_0$ について解くと, 区間

$$C(X) = \{\theta_0|T(X) \in A(\theta_0)\}$$

が得られ, これは

$$P_{\theta_0}(\theta_0 \in C(X)) = 1-\alpha$$

となる. ちなみにこの $P_{\theta_0}(\theta_0 \in C(X))$ を 被覆確率 (covarage probability)という.

しかし, このあたりの話題は具体例を複数やらないとちょっとわかりにくいかもしれない.

関連：ウィルコクソン検定から得られる信頼区間

ここでは1標本の二項分布という（実用性を失わない範囲でたぶん一番簡単な）モデルを題材に信頼区間を求めてみる.

地道に手計算してみる

準備

尤度に基づく方法を考える.

モデルを

$$X \sim \mathrm{Binomial}(n, \theta)$$

とする.

推定したい未知パラメータは $\theta \in [0,1]$ である.

このときこのモデルの対数尤度関数は

$$\ell(\theta) = X \log(\theta) + (1-X) \log(1-\theta) + \mathrm{Const.}$$

と書ける. $\mathrm{Const.}$ は自由パラメータ $\theta$ に依存しない定数で微分すると消える.

対数尤度関数をパラメータで1階微分した導関数をスコア関数と呼ぶ.

スコア関数 $S(\theta)$ は

$$S(\theta) = \frac{d}{d\theta}\ell(\theta) = X/\theta - (1-X)/(1-\theta).$$

これを0とおいて $\theta$ について解くことにより, 最尤推定量は

$$\hat \theta = X/n$$

と求まる.

さらに,　フィッシャー情報量を最尤推定量のサンプルサイズ $n$ の大きいときの分布（漸近分布）のために考えておく.

フィッシャー情報量 $I(\theta)$ は対数尤度関数の2階微分の期待値の符号反転である.

$$\begin{align*} I(\theta) &= -E\left[\frac{d^2}{d\theta^2}\ell(\theta)\right] \\ &= -(E[-X /\theta^2 - (1-X)/(1-\theta)^2])\\ &= -n(-1/\theta - 1/(1-\theta))\\ &= \frac{n}{\theta(1-\theta)}. \end{align*}$$

尤度比検定

帰無仮説 $\theta=\theta_0$, 対立仮説 $\theta \neq \theta_0$ の尤度比検定では $\theta = \theta_0$ と固定した自由パラメータが0個のモデルと, 自由パラメータが1個のモデルで尤度の比の対数の2倍が自由度 $1-0 = 1$ のカイ二乗分布に従うことを利用する.

検定統計量は

$$T(X) = 2(\ell(\hat \theta) -\ell(\theta_0))$$

で, これが自由度1のカイ二乗分布の上側確率 $\alpha$ を与える点以下の範囲の $\theta_0$ が信頼区間である.

しかし閉じた形では求まらなかった.

Wald 検定

帰無仮説 $\theta=\theta_0$, 対立仮説 $\theta \neq \theta_0$ の Wald 検定は検定統計量

$$T(X) = |\hat \theta - \theta|\sqrt{I(\hat \theta)}$$

が帰無仮説の下で漸近的に標準正規分布にしたがうことを利用する. すなわち, 帰無仮説の下で

$$P_{\theta_0}(T(X) < z_{\alpha/2}) = 1-\alpha$$

ここで $z_{\alpha/2}$ は標準正規分布の上側確率 $\alpha/2$ を与える分位点である.

信頼区間はこれを $\theta_0$ について解くことで,

$$C(X) = [C_{-}, C_{+}]$$

と求まる.

ただし,

$$C_{-} = \hat \theta - z_{\alpha/2}\sqrt{\hat \theta(1 - \hat\theta)/n},$$

$$C_{+} = \hat \theta +z_{\alpha/2}\sqrt{\hat \theta(1- \hat\theta)/n}$$

とした.

スコア検定

帰無仮説 $\theta=\theta_0$, 対立仮説 $\theta \neq p_0$ のスコア検定は検定統計量

$$T(X) = |S(\theta_0)|/\sqrt{I(\theta_0)}$$

が帰無仮説の下で漸近的に標準正規分布にしたがうことを利用する.

この場合について書き下すと,

$$T(X) = \sqrt{n} \, \frac{|\hat \theta - \theta_0|}{\sqrt{\theta_0(1- \theta_0)}}.$$

$T(X) \le z_{\alpha/2}$ となる（受容域に入る）のは,

$$(\hat \theta - \theta_0)^2 < z_{\alpha/2}^2 \,\frac{\theta_0(1- \theta_0)}{n}$$

のときである.

この不等式は次のようにも書ける.

$$\left(1+\frac{z_{\alpha/2}^2}{n}\right)\theta_0^2 - 2\left(\hat\theta\theta_0 + \frac{z_{\alpha/2}}{n}\right)\theta_0 + (\hat \theta)^2 \le 0$$

なので $\theta_0$ について解くことは2次方程式の解を求めることであり, いくらかの計算の後,

$$C(X) = [C_{-} , C_{+}]$$

ただし,

$$C_{-} = (1+\frac{z_{\alpha/2}}{n})^{-1} \left( \hat \theta - \frac{z_{\alpha/2}^2}{2n} - z_{\alpha/2}\sqrt{\frac{z_{\alpha/2}}{4n^2} + \frac{\hat \theta (1-\hat \theta)}{n}} \right),$$

$$C_{+} = (1+\frac{z_{\alpha/2}}{n})^{-1} \left( \hat \theta + \frac{z_{\alpha/2}^2}{2n} + z_{\alpha/2}\sqrt{\frac{z_{\alpha/2}}{4n^2} + \frac{\hat \theta (1-\hat \theta)}{n}}\right)$$

とした.

Rによる実装例

尤度比検定

尤度比検定に基づく p 値を返す関数は例えば次のように書ける.

LRtest <- function(p,x,n){
  phat <- x/n
  lp <- 2*(dbinom(x, n, phat ,log = TRUE)-
             dbinom(x, n, p ,log = TRUE))
  pchisq(lp, 1, lower.tail = FALSE)
}

信頼区間は閉じた形では求まらなかったがせっかくコンピュータを使うので帰無仮説を動かして総当たり（0.005刻み）で p 値を計算してみよう.

set.seed(123)
n <- 20
x <- rbinom(1, n, 0.5)
print(x)
# [1] 9

z <- seq(0.01, 0.99, by=0.005)
pv_lr <- sapply(z, LRtest, x=x, n=n)
df_lr <- data.frame(p=z, pv=pv_lr, method="LR")

library(ggplot2)
ggplot(df_lr, aes(x=p, y=pv))+
  geom_line()+
  theme_bw(18)+labs(x="param.", y="p-value")

このグラフを例えば 0.05 の高さで切ってやれば95%信頼区間が求まる.

実は数値的に解くのはさほど大変ではない.

こんなふうだ:

library(rootSolve)
sol <- uniroot.all(function(p)LRtest(p,x,n)-0.05, c(0.1,0.8)) #0になる点を求める
print(sol)
#[1] 0.2479710 0.6641676
ggplot(df_lr, aes(x=p, y=pv))+
  geom_line()+
  geom_errorbarh(data = NULL, aes(xmin=sol[1], xmax=sol[2], y=0.05), height=0.03, colour="cornflowerblue")+
  theme_bw(18)+labs(x="param.", y="p-value")

Wald 検定

Wald 検定に基づく p 値を返す関数は例えば次のように書ける.

waldtest <- function(p,x,n){
  phat <- x/n
  se <- sqrt(phat*(1-phat)/n)
  2*pnorm(abs(phat-p)/se, lower.tail = FALSE)
}

pv_w <- sapply(z, waldtest, x=x, n=n)
df_w <- data.frame(p=z, pv=pv_w, method="Wald")

abs で絶対値を取って片側だけ計算しているので2倍している.

スコア検定

R の prop.test はスコア検定に基づく. このことは次のように確かめられる.

set.seed(1234)
res_p <- prop.test(x, n, p = 0.5, correct = FALSE, conf.level = 0.95)

CI_score <- function(x, n, level){
  z <- qnorm(0.5*(1-level), lower.tail = FALSE)
  phat <- x/n
  t_1 <- phat+(z^2)/(2*n)
  t_2 <- z*sqrt(z^2/(4*n^2)+phat*(1-phat)/n)
  c((t_1 - t_2)/(1+(z^2)/n),
    (t_1 + t_2)/(1+(z^2)/n))  
}


print(CI_score(x,n,0.95))
print(res_p$conf.int)

> print(CI_score(x,n,0.95))
[1] 0.2581979 0.6579147
> print(res_p$conf.int)
[1] 0.2581979 0.6579147
attr(,"conf.level")
[1] 0.95

帰無仮説を動かしたときの p 値を 3 つの検定で比較してみよう.

pv_p <- sapply(z, function(p)prop.test(x, n, p = p, correct = FALSE)$p.value)
df_p <- data.frame(p=z, pv=pv_p, method="score")
df_pv <- rbind(df_lr,df_w,df_p)
ggplot(data = df_pv, aes(x=p, y=pv, colour=method, group = method, linetype=method))+
  geom_line()+
  scale_color_brewer(palette = "Set2")+
  theme_bw(16)+
  labs(x="param.", y="p-value", colour="method", linetype="method")
ggsave("proptest.png")

より大きい p 値が出やすい検定がより広い信頼区間を与える.

これだけではどの検定が良い検定かはまだわからない（次の記事に続く予定…）.

しかし p 値が計算できれば検定も区間推定もできることがわかった.

Rのコード全体はこちら：

https://github.com/abikoushi/Zenn_content/blob/28b54a4fa9641c360d8df0fb285b29e8d5b9f78c/R/lr_wald_score.R