データ解析のための統計モデリング入門 GLMの尤度比検定と検定の非対称性


推定された統計モデルを比較する検定を行う．
本章では**尤度比検定**を扱う．


# 統計学的な検定の枠組み

手順

1. 使用するデータを確定する
2. 目的とデータの構造に対応した適切な統計モデルを設計し，パラメータの最尤推定をする

> パラメータの少ないモデルと多いモデル（単純モデルと複雑モデル）を検定では，
**帰無仮説**と**対立仮説**とよぶ．

### Neyman-Pearsonの検定の枠組み

統計モデルの検定では，「帰無仮説は正しい」という命題を否定できるかどうかを調べる．

1. モデルの当てはまりの良さなどを**検定統計量**として扱う
2. 帰無仮説が「真のモデル」であると仮定する
3. 検定統計量のばらつき（確率分布）を調べて，検定統計量の値が”ありがち範囲”を決める
4. ”ありがちな範囲”の大きさが95％のとき，5％の有意水準を設定した言える
5. 対立仮説のモデルで得られた検定統計量が”ありがちな範囲”からはみ出ていれば，帰無仮説を棄却し，対立仮説を支持する


# 尤度比検定の例題：逸脱度の差を調べる

使用する統計モデル：
 
 - **一定モデル**：種子数の平均$\lambda_i$が定数であり，体サイズ$x_i$に依存しないモデル（傾き$\beta_2=0$；パラメータ数$k=1$）
 - **xモデル**：種子数の平均$\lambda_i$が体サイズ$x_i$に依存するモデル（傾き$\beta_2 \neq 0$；パラメータ数$k=2$）


| モデル | パラメータ数$k$ | $\log L ^ \ast$ | deviance（$-2\log L ^ \ast$）| residual deviance |AIC |
|:------|--------------:|----------------:|----------------------------:|------------------:|--:|
| 一定   |             1 |          -237.6 |                       475.3 |              89.5 |477.3|
| x     |             2 |          -235.4 |                       470.8 |              85.0 |474.8|
| full  |           100 |          -192.9 |                       385.6 |               0.0 |585.8|


尤度比検定では，尤度比の対数に-2をかけた逸脱度の差を用いる．本例題の場合は，

```math
\begin{eqnarray}
\frac{L_1 ^\ast}{L_2 ^\ast} &=& \frac{一定モデルの最大尤度：\exp(-237.6)}{xモデルの最大尤度：\exp(-235.4)} \\
\Delta D_{1,2} &=& -2 * (\log L_1 ^ \ast - \log L_2 ^ \ast) \\
\Delta D_{1,2} &=& D_1 - D_2 \\
\Delta D_{1,2} &=& 4.5
\end{eqnarray}
```

この逸脱度が4.5改善されたことを評価する．

- 帰無仮説：一定モデル（パラメータ数$k=1, \beta_2=0$）
- 対立仮説：xモデル（$k=2,\beta_2 \neq 0$）



| 帰無仮説は | $\Delta D_{1,2}$はめったに無い差（棄却） | $\Delta D_{1,2}$はよくある差（棄却できない） |
|:------|:--------------:|:----------------:|
| 真のモデルである   | **第一種の過誤** |  （問題なし） |
| 真のモデルではない  | （問題なし） |  **第二種の過誤** |


## Neyman-Pearsonの検定の枠組みでは**第一種の過誤の検討のみに専念**
手順としては

1. 帰無仮説である一定モデルが正しいもの（真のモデルである）と仮定する
2. 観測データに一定モデルを当てはめると，$\hat{\beta_1}=2.06$になったから，これが真のモデルと一緒と考える．
3. 真のモデルからデータを作成，そのたびにxモデルを当てはめれば$\Delta D_{1, 2}$を計算でき，$\Delta D_{1,2}$の分布がわかる
4. 一定モデルとxモデルの逸脱度の差が$\Delta D_{1,2} \geq 4.5$となる確率Pが求められる．
5. $\Delta D_{1,2}=4.5$がありえない値とみなされた場合，帰無仮説は棄却され，対立仮説が自動的に採択される．


一定モデルとxモデルの逸脱度の差が$\Delta D_{1,2} \geq 4.5$となる確率$P$を**P値**と呼ぶ．

- P値が「大きい」：$\Delta D_{1,2}=4.5$はよくあること → 帰無仮説棄却できない
- P値が「小さい」：$\Delta D_{1,2}=4.5$はとても珍しいことだな→帰無仮説を棄却しよう，残ったxモデルを正しいとする

**有意水準**を事前に（自分で）決め，

- $P \geq \alpha$：帰無仮説は棄却できない
- $P < \alpha$：帰無仮説は棄却できる

とする．

帰無仮説：一定モデルが真のモデルである世界における検定統計量である$\Delta D_{1,2}が4.5より大きくなる（第一種の過誤をおかす）確率を算出する．


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

上記３．のデータを生成する過程を乱数によるシュミレーションにもとづいて行う方法．
`fit1`に一定モデルの結果が，`fit2`にxモデルの結果が格納されているとする．


```pycon
>>> model1 = smf.glm('y~1', data=d, family=sm.families.Poisson())
>>> model2 = smf.glm('y~x', data=d, family=sm.families.Poisson())
>>> fit1 = model1.fit()
>>> fit2 = model2.fit()

# fit1$deviance - fit2$deviance
>>> fit1.deviance - fit2.deviance
4.5139410788540459
```

一定モデルによって推定された平均種子数$\lambda=\exp(2.06)=7.85$より，
生成すべきデータは「**平均が7.85である100個のポアソン乱数**」となる．


```pycon
# d$y.rnd <- rpois(100, lambda=mean(d$y))
>>> d.y.rnd = sci.poisson.rvs(d.y.mean(), size=100)
# fit1 <- glm(y.rnd ~ 1, data=d, family=poisson)
>>> model1 = smf.glm('y.rnd ~ 1', data=d, family=sm.families.Poisson())
>>> fit1 = model1.fit()
# fit2 <- glm(y.rnd ~ x, data=d, family=poisson)
>>> model2 = smf.glm('y.rnd ~ x', data=d, family=sm.families.Poisson())
>>> fit2 = model2.fit()
# fit1$deviance - fit2$deviance
>>> fit1.deviance - fit2.deviance
0.63270346107955788
```


平均値が一定のポアソン乱数であるデータに対して，逸脱度の差が1.92となる．

> 「真のモデル」である一定モデルよりも無意味な説明変数をもつxモデルのほうが当てはまりが良い

このステップを1000回ほど繰り返せば，**検定統計量の分布**，この例題で言う「逸脱度の差$\Delta D_{1, 2}$の分布」を予測できる．

```py
def get_dd(d):
    sample_d = len(d)
    mean_y = d.y.mean()
    d.y.rnd = sci.poisson.rvs(mu=mean_y, size=sample_d)
    model1 = smf.glm('y.rnd ~ 1', data=d, family=sm.families.Poisson())
    model2 = smf.glm('y.rnd ~ x', data=d, family=sm.families.Poisson())
    fit1 = model1.fit()
    fit2 = model2.fit()
    return fit1.deviance - fit2.deviance


def pb(d, bootstrap=1000):
    return pandas.Series([get_dd(d) for _ in xrange(bootstrap)])

results = pb(d)
results.describe()

results[results >= 4.5].count()


results[results>= 4.5].count()
# 32

results.quantile(.95)
# 3.6070094508948025

results.hist()


```

上記の結果から，「逸脱度の差が4.5より大きくなる確率」は$32/1000$，すなわち$P=0.032$となり，
また$P=0.05$となる逸脱度の差$\Delta D_{1,2}$を求めると．$\Delta D_{1,2} \leq 3.61$となる．

以上から，
「逸脱度の差4.5のP値は0.032だったので，これは有意水準0.05よりも小さい」ので有意差があり，
「帰無仮説（一定モデル）は棄却され，xモデルが残るのでこれを採択」と判断する．





## カイ二乗分布を使った近似計算法


逸脱度の差$\Delta D_{1,2}$の確率分布は，自由度$k_2 - k_1 = 2 - 1= 1$の$\chi^2$分布で近似できる場合がある．

> どうやら`statsmodels`にはGLMに対してANOVAはないようで，，，

$\chi^2$分布近似は，**サンプルサイズが大きい場合に有効な近似計算**である．

調査した個体数が多くない小標本，もしくはデータのばらつきがポアソン分布ではなく，等分散正規分布の場合には，PB法がよいらしい．


# まとめ

尤度比検定とAICによるモデル選択では，両方とも逸脱度という統計量に注目している．

- AIC：「良い予測をするモデル」の選択が目的
- 尤度比検定：帰無仮説の安全な棄却が目的

目的に合わせて選ぶしかないらしい．