Blog

Pythonで「データ解析のための統計モデリング入門」:第5章

Pythonで緑本(「データ解析のための統計モデリング入門」)の演習をしていきます。前記事はこちら

第5章は「GLMの尤度比検定と検定の非対称性」ということで、統計モデルの逸脱度の差に注目する尤度比検定についてです。

統計モデルの検定における「帰無仮説・対立仮説」が、AICによるモデル選択では「単純モデル・複雑モデル」となります。一定モデル(単純モデル)とxモデル(複雑モデル)の逸脱度の差が例えば4.5だったとして、尤度比検定では、検定統計量であるこの逸脱度の差が「4.5ぐらいでは改善されていない」と言ってよいのかどうかを調べます。逸脱度の差は以下の式で表せます

\[ \Delta D_{1,2} = -2 \times (\log L_{1} – \log L_{2}) \]

検定の全体の流れは以下のようになります。

(1)まずは帰無仮説である一定モデルが正しいものだと仮定する

(2)観測データに一定モデルをあてはめると、 (\hat{\beta_{1}} = 2.06 ) となったので、これは真のモデルとほぼ同じと考えよう

(3)この真のモデルからデータを何度も生成し、そのたびに ( \beta_{2} = 0 (k = 1) ) と( \beta_{2}\neq 0 (k = 2) ) のモデルをあてはめれば、たくさんの ( \Delta D_{1,2} ) が得られるから、 ( \Delta D_{1,2} ) の分布がわかるだろう

(4)そうすれば、一定モデルとxモデルの逸脱度の差が( \Delta D_{1,2}\ge 4.5 )となる確率 ( P ) が評価できるだろう

 

5.4.1 方法(1) 汎用性のあるパラメトリックブートストラップ法

パラメトリックブートストラップ法は、前述した検定の流れの(3)における「データを何度も生成し」の仮定を、乱数発生のシミュレーションによって実施する方法です。

この章の例題データのGLMによる推定結果は、一定モデルとxモデルそれぞれfit1とfit2に格納されているとします。

一定モデルとxモデルの逸脱度の差 ( \Delta D_{1,2}\ge 4.5 ) を計算します。

from __future__ import division, print_function, absolute_import, unicode_literals
import pandas as pd

d = pd.read_csv('data/data3a.csv')

import statsmodels.api as sm
import statsmodels.formula.api as smf 

fit1 = smf.glm(formula='y ~ 1', data=d, family=sm.families.Poisson()).fit()
fit2 = smf.glm(formula='y ~ x', data=d, family=sm.families.Poisson()).fit()

print(fit1.deviance - fit2.deviance)

4.51394107885

真のモデルから100個体ぶんのデータを新しく生成し、一定モデルとxモデルをこの新データにあてはめてみます。

import numpy as np

d['y_rnd'] = np.random.poisson(d.y.mean(), 100)

fit1 = smf.glm(formula='y_rnd ~ 1', data=d, family=sm.families.Poisson()).fit()
fit2 = smf.glm(formula='y_rnd ~ x', data=d, family=sm.families.Poisson()).fit()

print(fit1.deviance - fit2.deviance)

0.938292290633

今回は、逸脱度の差が0.94となりました。この1ステップによって「一定モデルが真のモデルである世界」での逸脱度の差がひとつ得られます。このステップを1000回ほど繰り返すと、「統計検定量の分布」、この例題でいうと「逸脱度の差 ( \Delta D_{1,2} ) の分布」を予測できます。

このPB法を実行するために、自作関数pb()を定義します。

def get_dd(d):
    d['y_rnd'] = np.random.poisson(d.y.mean(), len(d.y))
    fit1 = smf.glm(formula='y_rnd ~ 1', data=d, family=sm.families.Poisson()).fit()
    fit2 = smf.glm(formula='y_rnd ~ x', data=d, family=sm.families.Poisson()).fit()
    return fit1.deviance - fit2.deviance

def pb(d, n_bootstrap):
    return pd.Series([get_dd(d) for i in range(n_bootstrap)])

dd12 = pb(d, 1000)

以上の操作によって、逸脱度の差 ( \Delta D_{1,2} ) のサンプルが1000個つくられてdd12に格納されました。実は ( \Delta D_{1,2} ) のサンプル個数は ( 10^3 ) ぐらいでは十分ではなく、( 10^4 ) かそれ以上にしたほうが良いとのことです。

概要を調べます。

print(dd12.describe())

count 1000.000000
mean 0.983494
std 1.508326
min 0.000002
25% 0.092785
50% 0.456521
75% 1.259388
max 16.785845
dtype: float64

ヒストグラムで図示、そして逸脱度の差が4.5がどのあたりにくるのかを示します。

import matplotlib.pyplot as plt

fig = plt.figure()
ax = fig.add_subplot(111)
ax.hist(dd12, 100)
ax.set_xlabel(r'$\Delta D_{1,2}$')
plt.axvline(x=4.5, linestyle='--')
fig.savefig("itsudatsu.png")
plt.close(fig)

itsudatsu

合計1000個ある ( \Delta D_{1,2} ) のうちいくつぐらいが、この4.5より右にあるのでしょうか?

print((dd12 >= 4.5).sum())

31

1000個中の31個が4.5より大きい、つまり確率は31/1000、すなわち ( P = 0.0031 ) ということになります。ついでに、( P = 0.05 ) となる逸脱度の差 ( \Delta D_{1,2} ) を調べてみると、

print(dd12.quantile(0.95))

3.77686899182

となり、有意水準5%の統計学的検定のわくぐみのもとでは、 ( \Delta D_{1,2}\leq 3.78) ぐらいまでは「よくある差」とみなされます。

この尤度比検定の結論としては、「逸脱度の差4.5の ( P ) 値は0.031だったので、これは有意水準0.05よりも小さい」ので有意差が有り、「帰無仮説(一定モデル)は棄却され、xモデルが残るのでこれを採択」と判断します。

 

5.4.2 方法(2)  \( \chi^2 \) 分布を使った近似計算法

近似計算法を使うと、もっとお手軽に尤度比検定ができる場合があります。

逸脱度の差 ( \Delta D_{1,2} ) の確率分布は、自由の1の ( \chi^2 ) 分布に近似できる場合があります。Rではanova()関数を使い、anova(fit1, fit2, test=”Chisq”) といったふうにすれば ( P ) 値が得られますが、statsmodelsのanova関数(statsmodels.stats.anova.anova_lm())はGLMの結果に対応していません。

ですが、つまり5.4.1で得た逸脱度の差の確率分布を自由の1の ( \chi^2 ) 分布に近似して作って、その後は同様に計算すればいけるってことですよね・・・?

fit1 = smf.glm(formula='y ~ 1', data=d, family=sm.families.Poisson()).fit()
fit2 = smf.glm(formula='y ~ x', data=d, family=sm.families.Poisson()).fit()

chi = pd.Series(np.random.chisquare(1, 1000))

print((chi >= 4.5).sum() / 1000)

print(chi.quantile(0.95))

0.036
3.99139774575

得られた( P ) 値は0.0036です。ただ、( \chi^2 ) 分布近似はサンプルサイズが大きい場合に有効な近似計算なので、( P ) 値はあまり正確ではない可能性があります。

調査した個体数が多くない小標本のもとでは、PB法を使って逸脱度差の分布をシミュレーションで生成するのがいいでしょう、とのことです。

This Post Has 0 Comments

Leave A Reply