Pythonで実装しながら緑本を学ぶ (第5章 GLMの尤度比検定と検定の非対称性)

データ解析のための統計モデリング入門(通称、緑本)を読み進めています。
述べられている理論を整理しつつ、Rでの実装をPythonに置き換えた際のポイントなども深掘りしていきます。

今回は第5章です。実装は以下で公開しています。

introduction_to_machine_learning_with_python/chapter5.ipynb at master · ohke/introduction_to_machine_learning_with_python · GitHub

5 GLMの尤度比検定と検定の非対称性

検定は、推定された統計モデルを比較する方法のひとつ。どのようなモデルにも適用可能な検定として、尤度比検定を導入する。

5.1 統計学的な検定のわくぐみ

代表的な統計学的な検定の枠組みとして、Neyman-Pearsonの検定があります。この検定は以下の手順に従います。
1.と2.は第4章まで行ってきたモデル作りですが、3.以降は「帰無仮説が正しい」という命題が否定できるかどうかを調べます。

データを確定する
- 最後までそのデータを、常に全て、使う
統計モデルを作り、パラメータの最尤推定を行う
- (比較用の)パラメータの少ないモデルを帰無仮説、帰無仮説をネストしたパラメータの多いモデルを対立仮説と呼ぶ
帰無仮説が真のモデルと仮定して、検定統計量の値が取りうる範囲を定める
- 検定統計量は任意に設定する(今回は尤度比を使う)
- 取りうる値の範囲が95%の場合、5%の有意水準を設定したという
対立仮説のモデルで得られた検定統計量がこの範囲からはみ出ているか確認する
- はみ出ていれば、帰無仮説が棄却(対立仮説が支持)されたという

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

これまで使ってきた種子数データを使い、帰無仮説(一定モデル)が棄却できるかどうかを調べる。

一定モデル: $\lambda_i$ が定数のモデルで、 $\lambda_i=\exp(\beta_1)$
xモデル: $\lambda_i$ が定数と体サイズxで定まるモデルで、 $\lambda_i=\exp(\beta_1+\beta_2x)$

最尤推定(ステップ2.)まで行った結果が、以下表となります。xモデルのほうが逸脱度は改善されていますが、ネストしたモデルの場合、パラメータの多いほうが常に逸脱度は小さくなります。

モデル	最大対数尤度	逸脱度	Interceptの最尤推定値	xの最尤推定値
一定モデル	-237.64	475.28	2.06	-
xモデル	-235.39	470.78	1.29	0.076

検定統計量の計算のために、尤度比を導入する。尤度比は(帰無仮説の最大尤度÷対立仮説の最大尤度)で、値が大きいほど対立仮説の方が当てはまりが良いと言えます。

$\frac{\mbox{帰無仮説の最大尤度}}{\mbox{対立仮説の最大尤度}}=\frac{\mbox{一定モデルの最大尤度:}L_1^{*}}{\mbox{xモデルの最大尤度:}L_2^{*}}=\frac{\exp(-237.6)}{\exp(-235.4)}$

さらに式変形して、この値 $\Delta D_{1,2}$ を検定統計量とする。 $\Delta D_{1,2}$ は逸脱度の差となります。

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

表より、逸脱度の差は4.5でした。この差が、パラメータxを追加したことで改善されたのか、それともたまたま良いあてはまりが得られたのか(=改善されていないのか)を調べていきます。

5.3 2種類の過誤と統計学的な検定の非対称性

Neyman-Pearsonの検定では、2種類の誤り(過誤)が起こりうるが、第一種の過誤のみを防ぐことに注意する(第二種の過誤は検討されません)。検定の非対称性と呼ばれます。

帰無仮説が真のモデルであるにも関わらず、棄却される(第一種の過誤、type I error)
帰無仮説が真のモデルでないにも関わらず、棄却されない(第二種の過誤、type II error)

第一種の過誤のみを防げばいいので、観測データについて帰無仮説モデル上で起こる確率が低い場合に棄却するだけで良くなり、シンプルになります。

帰無仮説が正しいものと仮定して、真のモデルを作る
- ここでは一定モデルを真のモデルとする
この真のモデルからデータを生成し、帰無仮説モデルと対立仮説モデルに当てはめて、たくさんのを得る
- 例えば、100個の $\Delta_{1,2}$ を計算する
が起こる確率Pを評価する
- Pの計測方法は、この後2つ紹介する

5.4 帰無仮説を棄却するための有意水準

確率Pの値によって帰無仮説を棄却できるかどうかを決定します。
Pが大きければ帰無仮説モデルでもよくあることなので棄却できませんし、Pが小さければ帰無仮説モデルでは珍しいので棄却できる(対立仮説モデルを支持する)ことになります。

この棄却できるかどうかのPの閾値αは、有意水準と呼ばれます。多くの場合、有意水準には0.05が使われます(この値自体は決め打ちで良く、特に根拠は無いようです)。

P値の計測方法を2つを紹介します。

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

パラメトリックブートストラップ法(PB法)はシミュレーションに基づいた方法で、帰無仮説モデルから乱数でデータを大量生成し、検定統計量(ここでは逸脱度の差 $\Delta_{1,2}$ )の確率分布を調べます。
生成された $\Delta_{1,2}$ 100個の内、4.5以上が12個であれば $P=0.12$ 、と計算します。αを0.05とすると、P<α (0.12<0.05)のため、帰無仮説は棄却できないと評価されます。

メリット
- どのような統計モデルでも適用できる
- 元のデータ数が小さくても、サンプルサイズを大きくすれば高い精度が得られる
デメリット
- 生成するデータ数を適切に設定する必要がある

Pythonでは以下のような実装となります。

def pb(data, n):
    # パラメトリックブートストラップ法(PB法)によるP値の計算
    const_model_data = data.copy() # 種子数データ100行(各行は種子数yと体サイズxの2列)をコピー
    pb_d = []

    # n個のデータセットからn個の逸脱度の差を計算
    for i in range(n):
        # 一定モデルから目的変数100個生成
        const_model_data['y'] = np.random.poisson(const_model_results.params['Intercept'], len(data))

        # 逸脱度の差を計算
        const_model_llf = smf.poisson('y ~ 1', data=const_model_data).fit(disp=0).llf
        x_model_llf = smf.poisson('y ~ x', data=const_model_data).fit(disp=0).llf
        pb_d.append(-2 * (const_model_llf - x_model_llf))
        
    return pb_d
   
# 1000個のデータセットで計算 
pb_d = pb(data, 1000)

# P値の計算
p = (pb_d > d).sum() / len(pb_d)

逸脱度の差は以下の分布となります。この時、4.5を超えるのは1000個中41個で、P値は0.041と計算されます。
α=0.05とすると、P<αのため、帰無仮説(一定モデル)は棄却され、対立仮説(xモデル)が支持されたと言えます。

f:id:ohke:20180202085436p:plain

ばらつきを考慮したデータセット数を適切に設定する必要があります。
100個のデータセットでP値を計測したところ、0.05以上(つまり帰無仮説が棄却できない)が29個でした。

f:id:ohke:20180203102532p:plain

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

検定統計量 $D_{1,2}$ のばらつきを、自由度1(パラメータ数1のため)の $\chi^2$ 分布で近似する方法です。

$\chi^2$ 分布近似は、元のデータ数が小さい場合に精度が高くありません。
その場合、PB法にするか、データのばらつきが(ポアソン分布ではなく)等分散正規分布に従うならば別の確率分布を使う方が、正確なP値が得られます。適切な確率分布、検定統計量の特性によって異なります。

平均の差を検定統計量とする場合は、t分布
分散比を検定統計量とする場合は、F分布

この近似計算ではANOVAが使われます。ANOVAについては、こちらの解説(PDF)が詳しいです。

https://www.jstage.jst.go.jp/article/kagakutoseibutsu/51/7/51_483/_pdf

Rでは、構築済みの一定モデルとxモデルを使って計算できます。

> fit1 <- glm(y ~ 1, data = d family = poisson)
> fit2 <- glm(y ~ x, data = d family = poisson)
> anova(fit1, fit2, test="Chisq")

statsmodelsでも類似のanova_lmが提供されていますが、GLMには適用できません(以下のコードはエラーとなります。おそらく線形回帰で求める残差二乗和ssrがプロパティに無いよ、と怒られてます[参考])。

import statsmodels.api as sm

# エラーになります
sm.stats.anova_lm(const_model_results, x_model_results, test='Chisq')

要するに、自由度1の $\chi^2$ 分布において $D_{1,2}$ が4.5(これはstatsmodelsで計算済み)を超える確率を求めれば良いので、以下で十分かと思います(自信がなく間違っているかもです)。
以下では、自由度1の時の逸脱度の差(ここでは4.5)の生存関数の値を求めています。つまり4.5以上の値が得られる確率です。Pは0.034で、α以下のため、棄却できると結論付けられます。

from scipy.stats import chi2

p = chi2.sf(x=d, df=1) # 自由度1のカイ二乗分布における逸脱度の差(4.5)以上の確率を計算

5.5 「帰無仮説を棄却できない」は「差がない」ではない

$P \ge \alpha$ の場合は、帰無仮説を棄却できないだけで、「帰無仮説が正しい」ということにはならない。

Neyman-Pearsonの検定では、第二種の過誤の確率 $\beta$ を評価しません。
第二種の過誤を検討する時は、帰無仮説を棄却できる確率 $1-\beta$ がよく使われます( $1-\beta$ は検定力と呼ばれます)。サンプルサイズを大きくすることなどで、この検定力を高めることができます。

5.6 検定とモデル選択，そして検定された統計モデルの解釈

尤度比検定では、帰無仮説の棄却のみに着目しているため、帰無仮説が棄却された後に残る対立仮説が良いモデルなのかどうかは解釈が必要となる。