尤度比検定：Rによる計算例 - ジョンとヨーコのイマジン日記

モチベーション：尤度比検定をベイズ版と普通のやつ（？）で比べてみたい

ネイマン・ピアソンの補題

ネイマン・ピアソンの補題は数理統計の本にはよく出ている（たとえば竹村『現代数理統計学』）.

新装改訂版現代数理統計学

作者:彰通, 竹村
学術図書出版社

Amazon

ここでも一応紹介するが特にオリジナリティのある書きかたではないので, 飛ばして次のセクションの計算例にすすんでも構わない.

まず使用する記号をいくつか定める.

データ（確率変数）をまとめて $\boldsymbol{x}=(x_1, x_2, \ldots , x_n)'$ と書くことにする.

$\boldsymbol{x}$ に対して 0 または 1 の値を取る関数を検定関数と呼ぶことにし, 検定関数 $S(\boldsymbol{x})$ が 1 のとき帰無仮説を棄却することにする.

尤度比検定を考えたいので, 検定関数 $S(\boldsymbol{x})$ を $\ell(\boldsymbol{x})>c$ のとき 1, $\ell(\boldsymbol{x}) \le c$ のとき 0 の値を取るものとし, $\ell(\boldsymbol{x})$ を尤度比

$\displaystyle \ell(\boldsymbol{x}) = f(\boldsymbol{x} | \theta_1)/ f(\boldsymbol{x} | \theta_0)$

とする. $c$ は分析者が決める定数であり, 帰無仮説の下での期待値

$\displaystyle E[S(\boldsymbol{x}) \mid f(\boldsymbol{x}|\theta_0)] = \alpha$

が望む $\alpha$ になるように $c$ を定める.

一方, 対立仮説の下での期待値

$\displaystyle E[S(\boldsymbol{x}) \mid f(\boldsymbol{x}|\theta_1)] = \beta$

は検出力（power; 検定力という表記もたまに見る）と呼ぶ.

ネイマン・ピアソンの補題は帰無仮説（ $f(\boldsymbol{x}|\theta_0)$ ）, 対立仮説（ $f(\boldsymbol{x}|\theta_1)$ ）, 有意水準 $\alpha$ が同じ検定の中では、尤度比検定の検出力が最大であることを示したものである.

さて, ベイズ統計では推定したい未知パラメータ（ $\theta_0$ や $\theta_1$ ）に対して事前分布（ $\phi_0(\theta)$ , $\phi_1(\theta)$ ）を設定する.

この場合, 尤度比を

$\displaystyle \ell(\boldsymbol{x}) = \frac{\int f(\boldsymbol{x}|\theta)\phi_1(\theta) \, d\theta}{\int f(\boldsymbol{x}|\theta_0)\phi_0(\theta) \, d\theta}$

とすれば尤度が積分を使って書かれているだけなのでネイマン・ピアソンの補題の議論が成立する. ここでは全区間に渡る積分の積分区間を省略して書いてしまったがこの記事内でこれ以外の積分が出てくることはないので混乱はしないと思う.

諸注意＆余談：

ネイマン・ピアソンの補題では $\theta_0$ や $\theta_1$ に対して最尤推定量を代入した尤度比を用いた尤度比検定が強力であることは主張されていない.
ここでは「単純仮説」と呼ばれる1個のパラメータで表せるような仮説のみを扱っているが本当はもうちょっと一般化できる
離散型の分布の（データにタイがある）場合は=（イコール）が入るか入らないかのところをもうちょっと気にする必要がある
ネイマン・ピアソンの補題は通常の感覚では「補題」というよりは「定理」だと思う。ネイマンさんだかピアソンさんは最初これをなんとかの定理のための補題として示したので今でも補題と呼ばれるというのをどこかで聞いたんだけど、すべて忘れてしまった。

簡単な例：正規分布

$\boldsymbol{x}$ を独立に同一の分布に従うとし, 2つの正規分布 $\mathcal{N}(\mu_0,1)$ , $\mathcal{N}(\mu_1,1)$ を比べることにする.

尤度比は標本平均を $\bar x$ で書くと,

$\displaystyle \ell(\boldsymbol{x}) = \exp \left( \{\mu_1 - \mu_0\} n \bar x - \frac{n}{2}(\mu_1^2-\mu_0^2)\right)$

である. これを記事内では conditional な尤度比と呼ぶことにしよう.

次にこれのベイズ版を考える. ここでは帰無仮説の分布に $\mu_0$ の一点のみで密度を持つ分布, 対立仮説の分布に区間 $[\mu_0,\mu_0+b$ ] の一様分布を採用することにする.

尤度比は, 標準正規分布の累積分布関数を $\Phi (x)$ で表すことにすると,

$\displaystyle \ell(\boldsymbol{x}) = \frac{\Phi(\mu_0+b-\bar{x})-\Phi(\mu_0-\bar{x})}{\exp(n\{\mu_0 - \bar{x}\}/2)}$

と書ける. これを記事内では marginal な尤度比と呼ぶことにしよう. （この conditional/marginal な尤度比という用語はあまり一般的ではないが, 他に適当な言い方を思いつかなかった.）

この例のモデルはやや単純すぎるが実用性がないわけでもない.

例えば「対応のある（paired）」なんとか検定というのは対応のあるサンプルどうしの差をとって1グループの検定に帰着させるものである.

さて, 有意水準が $\alpha$ となるような $c$ は, 本当は解析的に求まるのだが, ここでは乱数を使って楽をする.

つまり, 正規分布 $\mathcal{N}(\mu_0,1)$ に従う乱数で尤度比の帰無分布を経験的に作っておき, その経験分布から p 値を求める.

R で書くとこんな感じだ。

#conditional な尤度比
lr1 <- function(x,mu1,mu0=0){
  (mu1-mu0)*sum(x)-(mu1^2-mu0^2)*length(x)/2
}

#marginal な尤度比
lr2 <- function(x,b,mu0=0){
  xbar <- mean(x)
  log((pnorm(mu0+b-xbar)-pnorm(mu0-xbar))/exp(-(length(x)*(mu0-xbar)^2)/2))
}

pvsimfunc <- function(n,mu,mu1,b,iter=10000){
  z <- matrix(rnorm(10000*n),10000,n)
  null_lr1 <- ecdf(apply(z,1,lr1,mu1=mu1))
  null_lr2 <- ecdf(apply(z,1,lr2,b=b))
  pv1 <- numeric(10000)
  pv2 <- numeric(10000)
  for(i in 1:10000){
    x <- rnorm(n,mu)
    pv1[i] <- 1-null_lr1(lr1(x,mu1))
    pv2[i] <- 1-null_lr2(lr2(x,b))
  }
  return(data.frame(pv1=pv1,pv2=pv2))
}
set.seed(2024)
res_a <- pvsimfunc(10,0,1,1)
res_b1 <- pvsimfunc(10,1,1,1)
res_b2 <- pvsimfunc(10,0.5,1,1)
res_b3 <- pvsimfunc(10,1.5,1,1)
res_b4 <- pvsimfunc(10,-1,1,1)

以降, $\mu_0=0$ , $mu_1=1$ または $b=1$ , サンプルサイズを10として, まずαエラー（帰無仮説が正しいときに帰無仮説を棄却してしまう確率）のシミュレーション：

これはどちらも名目どおりの水準になっている。

次に検出力のシミュレーション, 真の $\mu$ が 1 （対立仮説の設定と同じ）のとき:

$\mu=0.5$ （対立仮説の設定より小さめ）のとき:

$\mu=1.5$ （対立仮説の設定より大きめ）のとき:

$\mu=-1$ （対立仮説の設定と反対向き）のとき:

検出力は conditional のほうが高い.

marginal のほうは式の形（ $\mu_0=0$ なら $\displaystyle \ell(\boldsymbol{x}) = \{\Phi(b-\bar{x})-\Phi(-\bar{x})\}/\{\exp(-n \bar{x}/2)\}$ ）から両側検定になっている.

図のコードはこちら：

library(ggplot2)

p <- ggplot(res_a)+
  stat_ecdf(aes(x=pv1, colour="conditional"))+
  stat_ecdf(aes(x=pv2, colour="marginal"))+
  geom_abline(slope=1,intercept = 0,linetype=2)+
  theme_bw()+
  labs(colour="method", x="p-value")
print(p)
# ggsave("p_alpha.png")

p %+% res_b1 + ggtitle("mu = 1")
# ggsave("p_1.png")

p %+% res_b2 + ggtitle("mu = 0.5")
# ggsave("p_2.png")

p %+% res_b3 + ggtitle("mu = 1.5")
# ggsave("p_3.png")

p %+% res_b4 + ggtitle("mu = -1")
# ggsave("p_4.png")