データ解析のための統計モデリング入門 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をかけた逸脱度の差を用いる．本例題の場合は，

\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モデルの結果が格納されているとする．

>>> 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個のポアソン乱数」となる．

# 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}$の分布」を予測できる．

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：「良い予測をするモデル」の選択が目的
• 尤度比検定：帰無仮説の安全な棄却が目的

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