データ解析のための統計モデリング入門 GLMのモデル選択

良いモデルとは

線形予測子が切片だけのモデル（$\log \lambda = \beta_1$）は簡単すぎる．
パラメータ数を増やして，$\log \lambda = \beta_1 + \beta_2 x + \dots + \beta_6x^6$のようにすることで
当てはまりは改善されていく．

以下に多項式にした場合のソースコードを書いておく．

>>> import numpy
>>> import pandas
>>> import matplotlib.pyplot as plt
>>> import statsmodels.api as sm
>>> import statsmodels.formula.api as smf

>>> d = pandas.read_csv("data3a.csv")

>>> model = smf.glm('y ~ x + I(x**2) + I(x**3) + I(x**4) + I(x**5) + I(x**6)', data=d, family=sm.families.Poisson())
>>> result = model.fit()
>>> result.summary()

>>> params = result.params

>>> fit = lambda x: numpy.exp(params[0] + params[1] * x + params[2] * x ** 2 + params[3] * x ** 3 + params[4] * x ** 4 + params[5] * x ** 5 + params[6] * x ** 6)

>>> xx = numpy.linspace(min(d.x), max(d.x), 100)
>>> plt.plot(xx, fit(xx))

>>> plt.plot(d.x[d.f == 'T'], d.y[d.f == 'T'], 'ro')
>>> plt.plot(d.x[d.f == 'C'], d.y[d.f == 'C'], 'bo')
>>> plt.show()

AIC：「良い予測をするモデルが良いモデルである」

今まで種子のデータに対して4種類のモデルをあてはめた．

• 何も影響しないモデル（一定モデル）
• 体サイズの効果が影響するモデル（xモデル）
• 施肥効果が影響するモデル（fモデル）
• 体サイズの効果と施肥効果が影響するモデル（x+fモデル）

>>> plt.subplot(221)
>>> model = smf.glm('y ~ 1', data=d, family=sm.families.Poisson())
>>> result = model.fit()
>>> plt.plot(xx, [numpy.exp(result.params[0]) for _ in range(100)])

>>> plt.subplot(222)
>>> model = smf.glm('y ~ f', data=d, family=sm.families.Poisson())
>>> result = model.fit()
>>> plt.plot(xx, [numpy.exp(result.params[0] + result.params[1]) for _ in range(100)], 'r')
>>> plt.plot(xx, [numpy.exp(result.params[0]) for _ in range(100)], 'b')

>>> plt.subplot(223)
>>> model = smf.glm('y ~ x', data=d, family=sm.families.Poisson())
>>> result = model.fit()
>>> fit = lambda x: numpy.exp(result.params[0] + result.params[1] * x)
>>> plt.plot(xx, fit(xx))

>>> plt.subplot(224)
>>> model = smf.glm('y ~ x + f', data=d, family=sm.families.Poisson())
>>> result = model.fit()
>>> fit = lambda x: numpy.exp(result.params[0]+ result.params[2] * x)
>>> plt.plot(xx, fit(xx), 'b')
>>> fit = lambda x: numpy.exp(result.params[0]+ result.params[1] + result.params[2] * x)
>>> plt.plot(xx, fit(xx), 'r')

統計モデルのあてはまりの悪さ：逸脱度

当てはまりの良さである最大対数尤度$\log L ^ \ast$を変形した統計量逸脱度（deviance）$D$について

D = -2 \log L ^\ast

と定義する．
今までのresult.summary()にあったDevianceは，残差逸脱度（=そのモデルの逸脱度ー最小逸脱度，residual deviance）のこと．

最小逸脱度

フルモデルの逸脱度のこと．
フルモデルとは，$\lambda_i=y_i$とするモデルを指す．
つまり最大対数尤度は．

\log L = \sum_i \log \frac{y_i ^{y_i} \exp(-y_i}{y_i !}

xモデルで考えると，

# xモデルの最小逸脱度
>>> dpois = lambda x:sct.poisson.logpmf(x, x)
>>> sum(dpois(d.y))
-192.8897525244958
>>> sum(dpois(d.y)) * (-2)
385.77950504899161

以上のことから，このモデルの残差逸脱度は，

>>> result.llf * (-2) -  sum(dpois(d.y))*(-2)
84.992996490729922
>>> result.summary()
                 Generalized Linear Model Regression Results
==============================================================================
Dep. Variable:                      y   No. Observations:                  100
Model:                            GLM   Df Residuals:                       98
Model Family:                 Poisson   Df Model:                            1
Link Function:                    log   Scale:                             1.0
Method:                          IRLS   Log-Likelihood:                -235.39
Date:                 金, 05  6 2015   Deviance:                       84.993
Time:                        13:12:50   Pearson chi2:                     83.8
No. Iterations:                     7
==============================================================================
                 coef    std err          z      P>|z|      [95.0% Conf. Int.]
------------------------------------------------------------------------------
Intercept      1.2917      0.364      3.552      0.000         0.579     2.005
x              0.0757      0.036      2.125      0.034         0.006     0.145
==============================================================================

最大逸脱度

最小逸脱度の反対である最大逸脱度は，Nullモデルの逸脱度のこと．
Nullモデルとは，線形予測子が切片だけのモデルのこと．

まとめると，残渣逸脱度が大きいほど，当てはまりが悪く，小さいほど，あてはまりがよくなる（過学習？）

モデル選択基準：AIC

AIC（Akaike's information criterion）は予測の良さを重視する規準．
最尤推定したパラメータの数が$k$であるとき，

\begin{eqnarray}
AIC &=& -2 \{(最大対数尤度) - (最尤推定したパラメータ数) \} \\
&=& -2(\log L ^ \ast - k) \\
&=& D + 2k
\end{eqnarray}

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

表からxモデルがAIC最小の統計モデルとなる．

個人的なまとめ

一般に説明変数（パラメータ数）を増やすと，学習データに対して精度良くなるが，予測が悪くなる（ネスト関係にあるモデル間のバイアス差は一定）
→うまく説明できるものを選ぶと，バイアスは変わらないが，全体が良くなる→AICが下がる

最大対数尤度$\log L^ \ast$と平均対数尤度の差はパラメータ数増加（k=1 -> 2）による前者の改善が１より大きければ，一定モデル（k=1）よりAICが小さくなる．
この日本語が理解しきれていない，，，