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


# 良いモデルとは

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

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

```pycon
>>> 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モデル）

```pycon
>>> 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$について

```math
D = -2 \log L ^\ast
```

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

### 最小逸脱度

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

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

xモデルで考えると，

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

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

```pycon
>>> 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$であるとき，

```math
\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が小さくなる．
この日本語が理解しきれていない，，，