chapter5.Rmd

---
title: chapter5 ― GLMの尤度比検定と検定の非対称性
author: KADOWAKI, Shuhei
date: 2018/12/03
output:
    html_document:
        toc: true
        toc_float: true
        number_sections: false
        theme: cosmo
        code_folding: show
        df_print: paged
---

```{r echo=FALSE}
### Setting the global code chunk options ###
# Args
#   comment='': won't append any string to the start of each line of results
#   fig.align='center': align figures to the center of document
knitr::opts_chunk$set(comment="", fig.align="center")
```


# AICモデル選択 vs. Neyman-Pearson検定

- 尤度比検定(likelihood ratio test): 逸脱度の差に注目し, どのような統計モデルであっても, ネストしているモデルたちを比較できる検証法

- **Neyman-Pearson検定**: 「帰無仮説は正しい」という命題が否定できるかどうかを調べる検定の枠組み

- 共通手順(step.1 ~ step.3)
    1. 解析対象のデータを確定
    2. データを説明できるような統計モデルを設計
    3. ネストした統計モデルたちのパラメータの最尤推定計算
        - 「単純モデル」(AICモデル選択) vs. 「**帰無仮説(null hypothesis)**」(統計学的な検定)
        - 「複雑モデル」(AICモデル選択) vs. 「**対立仮説(alternative hypothesis)**」(統計学的な検定)

- 異なる手順(step.4 ~ step.5)
    - AICモデル選択の手順
        4. モデル選択基準AICの評価
        5. 予測の良いモデルを選ぶ
    - (Neyman-Pearson)検定の手順
        4. 帰無仮説棄却の危険率の評価
        5. 帰無仮説棄却の可否を判断

- Neyman-Pearson検定のstep.4 ~ step.5の詳細
    - step.4 帰無仮説棄却の危険率の評価
        - step.4.1 モデルの当てはまりの良さや逸脱度などを検定統計量(test statistics)に指定する
        - step.4.2 検定統計量の値が取りうる「ありがちな範囲」を定める
                - e.g.) 「ありがちな範囲」の大きさが95%である場合は, 「5%の**有意水準**(signficant level)」を設定したという
        - step.4.3 帰無仮説が「真のモデル」であると仮定して, そのときに検定統計量の理論的なばらつき(確率分布)を調べる
    - step.5 帰無仮説棄却の可否を判断
        - step.5.1 対立仮説のモデルで得られた検定統計量が, この「ありがちな範囲」からはみ出ているかどうかを確認
        - step.5.2 もしはみ出ていれば帰無仮説は棄却され, 対立仮説が支持されたと結論する
            - もしはみ出ていない(帰無仮説が棄却できない)場合は, 帰無仮説・対立仮説のどちらも**正しいとも正しくないとも言えない**と判断を保留する


# step.4.1 - 統計検定量の指定

例としてchapter3の`Nullモデル`と`xモデル`を比較する

- `Nullモデル`
    - 種子数の平均$\lambda_{i}$が定数: $\lambda_{i} = \exp(\beta_{1})$
    - パラメータ数1: $\beta_{1}$
- `xモデル`
    - 種子数の平均$\lambda_{i}$が体サイズ$x_{i}$に依存: $\lambda_{i} = \exp(\beta_{1} + \beta_{2}x_{i})$
    - パラメータ数2: $\beta_{1}$, $\beta_{2}$

ポアソン回帰の結果は以下のようになる

```{r, echo=FALSE, results='asis'}
d <- read.csv('../data/data3.csv')
 
fit.null <- glm(y ~ 1, family = poisson, data = d)
fit.x <- glm(y ~ x, family = poisson, data = d)
 
k.null <- 1
k.x <- 2
k.full <- 100
loglik.null <- logLik(fit.null)[1]
loglik.x <- logLik(fit.x)[1]
loglik.full <- sum(log(dpois(d$y, lambda = d$y)))
 
d.null <- -2 * loglik.null
d.x <- -2 * loglik.x
d.full <- -2 * loglik.full
 
cat(
    'model | $k$ | loglik $\\log L^*$ | D ($-2\\log L^*$) | residual dev | AIC', '\n',
    '--|--|--|--|--|--', '\n',
    'null', '|', k.null, '|', loglik.null, '|', d.null, '|', fit.null$deviance, '|', fit.null$aic, '\n',
    'x', '|', k.x, '|', loglik.x, '|', d.x, '|', fit.x$deviance, '|', fit.x$aic, '\n',
    'full', '|', k.full, '|', loglik.full, '|', d.full, '|', 0, '|', d.full + 2 * k.full, '\n'
)
```


## 尤度比とNeyman-Pearson検定における検定統計量
**尤度比**(likelihood ratio)は以下の式で定義される

$$
    \frac{L_{1}^*}{L_{2}^*} = \frac{Nullモデルの最大尤度: \exp(-237.6)}{xモデルの最大尤度: \exp(-235.4)}
$$

尤度比検定では, 尤度比の対数を取り-2をかけて, **逸脱度の差**に変換して統計検定量として扱う

$$
    \Delta D_{1, 2} = -2 (\log L_{1}^* - \log L_{2}^*) = D_{1} - D_{2}
$$

- $D_{1}$: `Nullモデル`の逸脱度  
- $D_{2}$: `xモデル`の逸脱度

上記の`Nullモデル`と`xモデル`の逸脱度の差$\Delta D_{1, 2}$は
```{r, echo=F, results='asis'}
dd <- fit.null$deviance - fit.x$deviance
cat(dd)
```
となっているおり, 尤度比検定ではこの逸脱度の差が「
```{r, echo=F, results='asis'}
cat(dd)
```
ぐらいでは改善されていない」と言ってよいのかどうかを調べる


## 2種類の過誤とNeyman-Pearson検定の非対称性

Neyman-Pearson検定では比較するモデルを次のように分類する

- 帰無仮説: `Nullモデル` ($k = 1$, $\beta_{2} = 0$)
    - 「棄却されるための仮説」であり, 「無に帰される」ときにのみその役割を果たす特殊な統計モデル
- 対立仮説: `xモデル` ($k = 2$, $\beta_{2} \neq 0$)
    - 対立仮設の意味としては「代替仮設」の方が訳語としては正しい(検定によって帰無仮説が「追放」(棄却)されたあと, 現象の説明を代替するために残されたモデルであるから)

このとき, 次のような2種類の過誤が予期される

 帰無仮説は/を | 棄却する  | 棄却しない
--|---|--
 真のモデルである | 第一種の過誤  |  (問題なし)
 真のモデルでない | (問題なし)  |  第二種の過誤

- **第一種の過誤**(type Ⅰ error): (帰無仮説が真のモデルであるのに)逸脱度の差が「ありえない値」だと判断して, 帰無仮説を棄却してしまう
- **第二種の過誤**(type Ⅱ error): (帰無仮説が真のモデルではないのに)逸脱度の差が「あえりえる値」だと判断して, 帰無仮説を棄却しない
- **検定の非対称性**: Neyman-Pearson検定では**第一種の過誤の検討にだけ専念する**こと
    - 「帰無仮説が正しい」という仮説が否定できるかどうかだけに注目することで以下のように検定に必要な計算を簡単にできる
        - step 4.1: まず帰無仮説(`Nullモデル`)が正しいと仮定し, $beta_{1}=2.06$を真のモデルと考える
        - step 4.2: 帰無仮説を棄却するための有意水準を設定する
        - step 4.3: このモデルからデータを何度も生成し, そのたびに帰無仮説(`Nullモデル`)と対立仮説(`xモデル`)をあてはめれば, たくさんの$\Delta D_{1, 2}$を得ることができ, $\Delta_{1, 2}$の分布がわかり, $\Delta D_{1, 2} \geq 4.51$となる確率$P$を評価できる

```{r, include=F}
cat(fit.null$coefficients[1])
cat(dd)
```

# step.4.2 - 有意水準の設定

- **$P$値**($P$ value): 第一種の過誤をおかす確率
    - e.g.)この例題では, $\Delta_{1, 2} \geq 4.51$となる確率
    - $P$値が大きい($\Delta_{1, 2} = 4.51$はよくあること) ⇒ 帰無仮説は棄却できない
    - $P$値が小さい($\Delta_{1, 2} = 4.51$は珍しい) ⇒ 帰無仮説は棄却し, 残った対立仮説を正しいとする
- **有意水準$\alpha$**: $P$値の大小を判断するために, 事前に決めておく量
    - $P \geq \alpha$: 帰無仮説は棄却できない
    - $P < \alpha$: 帰無仮説は棄却できる
    - **自分で好き勝手に決めるほかない**
        - $\alpha = 0.05$, つまり「めったにないことは, 20回のうち1回より少ない発生件数である」とする値がよく使われている
        - データをとる前の段階であらかじめ決めておくのが正しいお作法とされる


# step 4.3 - 統計検定量の確率分布, $P$値を調べる

## パラメトリックブーツストラップ法

- **パラメトリックブーツストラップ法(parametric bootstrap, PB method)**: いかなる面倒な状況でも必ず$P$値を計算する方法
    - 「データをたくさん生成する」のに乱数発生のシミュレーションを用いる
    - 具体的には次のようなステップを**繰り返す**
        - 帰無仮説を真のモデルとして, データを乱数を用いて生成する
        - 帰無仮説, 対立仮説をそのデータに対して当てはめる
        - 逸脱度の差を計算する

この例題データで帰無仮説(`Nullモデル`)で推定された平均種子数$\lambda$は
```{r, echo=F, results='asis'}
cat(exp(fit.null$coefficients))
```
であったので, 真のモデルを「平均7.83の100個のポアソン乱数」と仮定でき, Rのポアソン乱数生成関数`rpois`を用いてPB法を以下のように実装できる.

```{r}
get.dd <- function(data) {
    n.sample <- nrow(data)
    lambda <- 7.83
    d$y.rand <- rpois(n.sample, lambda = lambda)
    fit.rand.null <- glm(d$y.rand ~ 1, data = data, family = poisson)
    fit.rand.x <- glm(d$y.rand ~ x, data = data, family = poisson)
    fit.rand.null$deviance - fit.rand.x$deviance
}
pb <- function(data, n.boostrap) {
    replicate(n.boostrap, get.dd(data))
}
```

$\Delta_{1, 2}$を10000回計算し, 結果を図示すると以下のようになる ($\Delta_{1, 2}$のサインプルサイズは$10^3$ほどでは十分ではなく, 精度を高めるためには$10^4$以上にする方がよい)

```{r, results='hold', fig.show='hold'}
n.boostrap = 10000
dd12 <- pb(data = d, n.boostrap = n.boostrap)
summary(dd12)
hist(dd12, 100)
abline(v = 4.5, lty = 2)
```


得られた$\Delta_{1, 2}$のうちいくつが,
```{r, echo=FALSE, results='asis'}
cat(dd)
```
よりも大きいのかを数えると以下のように$P$値を計算できる
```{r, echo=T, results='asis'}
dd <- fit.null$deviance - fit.x$deviance
p <- sum(dd12 >= dd) / n.boostrap
```
```{r, echo=F, results='asis'}
cat('$$P = ', p, '$$')
```

よって, この尤度比検定の結論として, 「$\Delta_{1, 2}$の$P$値は
```{r, echo=F, results='asis'}
cat(p)
```
だったので, これは有意水準$\alpha = 0.05$よりも小さい」ので有意差があり(significantly different), 「帰無仮説(`Nullモデル`)は棄却され, 対立仮説(`xモデル`)が残るのでこれを採択」と判断する

- 有意差は説明変数が及ぼす効果の大きさだけで決まるわけではなく, サンプルサイズが大きい場合は小さな差でも統計学的には有意となる場合があることに注意
- $P$値は小さければ小さいほど有意というわけではなく, Neyman-Pearson検定の枠組みの下では$P < \alpha$となっているか, あるいはなっていないかだけが問題であることに注意
- 棄却点と棄却域も気にしてみる
    - 棄却点(critical point): $P = \alpha$となるような$\Delta_{1, 2}$
    - 棄却域(critical region, rejection region): 棄却点よりも大きい$\Delta_{1, 2}$の領域

棄却点はRでは以下のように調べることができる
```{r}
alpha <- 0.05
critical_point <- quantile(dd12, 1 - alpha)
cat('critical point:', critical_point)
```
つまり, 有意水準5%のNeyman-Pearson検定の枠組みの下では, $\Delta_{1, 2} \leq$
```{r, echo=F, results='asis'}
cat(critical_point)
```
ぐらいまでは「よくある差」とみなす


## $\chi^2$分布を使った近似計算法

PB法は計算量が非常に大きくなるが, 近似計算法を使うとよりお手軽に尤度比検定ができる場合がある  
近似計算法で得られた$P$値とPB法で得られる$P$値は一致しないことに注意

- $\Delta_{1, 2}$の確率分布を自由度$k$の$\chi^2$分布($\chi^2$ distribution)で近似できる場合がある
    - e.g.) 今回は自由度1となる
- $\chi^2$分布近似は**サンプルサイズが大きい場合に有効な近似計算**であり, サンプルサイズによってはあまり正確ではない可能性がある
    - **サンプルサイズが大きくない小標本の下では, PB法を使って逸脱度差分布の差をシミュレーションで生成する方がよい**
        - PB法の検定統計量のサンプルサイズは, 問題によるので試行錯誤して決めるしかない

Rでは以下のように, `anova`関数で引数`test = "Chisq"`とすることで利用可能  
(`Pr(>Chi)`が$P$値に対応)

```{r}
anova(fit.null, fit.x, test = "Chisq")
```

- **等分散正規分布**の場合には, サンプルサイズが小さい場合の検定統計量の確率分布を利用できる
    - 平均差を検定統計量とする場合: $t$分布
    - 分散比を検定統計量とする場合: $F$分布

    
# 帰無仮説を棄却できない場合($P \geq \alpha$となるとき)

- 「帰無仮説が棄却できない(failed to reject)」= 「帰無仮説が正しい」**ではない**
    - 「帰無仮説・対立仮説のどちらも**正しいとも正しくないともいえない**と判断を保留する
- Neyman-Pearson検定には非対称性があるので, $P < \alpha$となった場合と, $P \geq \alpha$となった場合では結論できることが大きく異なる
    - 以下のように第一種の過誤のみを排除できただけ, と解釈するべき

 帰無仮説は/を | 棄却する  | 棄却しない
--|---|--
 真のモデルである | ~~第一種の過誤 ($P$)~~  |  (問題なし)
 真のモデルでない | (問題なし)  |  第二種の過誤 ($P_{2}$)

- 第二種の過誤の確率を$P_{2}$と評価することもある
    - Neyman-Pearson検定の枠組みの下では, $P_{2}$を使って何かを定量的に主張する手続きは用意されていない
    - **検定力**(検出力, power): 「帰無仮説が間違っていたときに棄却できる」確率 $1 - P_{2}$


# モデル選択 vs. 検定 (again)

- モデル選択: 
    - 目的: 「良い予測をするモデルを選ぶ」 
    - 解釈: 「予測の良さは平均対数尤度」と明示し, 平均対数尤度を最大対数尤度とパラメータ数から推定する
- 検定: 
    - 目的: 「帰無仮説の安全な棄却」
    - 解釈: 帰無仮説が棄却されたあとに残された対立仮説が, どのような意味で「良い」モデルなのかは明確でない


モデル選択や検定のいずれを使用するにしろ, その結果だけに注目するのではなく, 推定された統計モデルが, 対象となる現象の挙動をどのように予測しているのかも確認するべきである  
- 推定されたパラメータがどのような値なのか(**効果の大きさ**(effect size))
- 標準誤差などで表される推定の誤差はどれほどなのか
- それらを組み合わせたときの統計モデルの挙動はどうなると予測されるのか
- などなど(特にサンプルサイズとの関連を考える)