データ解析のための統計モデリング入門 一般化線形モデル（GLM）


前回（[2章]()）に引き続いて，３章．
**説明変数**を組み込んだ統計モデルについての勉強．

# 3.1 例題：個体ごとに平均種子数が異なる場合

前回との差分として，個体の属性のひとつである体サイズ$x_i$とコントロールの有無（何もしていない：**C**，肥料を加える：**T**）を説明変数として用いる．

### データの読み込み

[ここ]()の３章のところの`data3a.csv`をダウンロードして使う．
今回使うデータはDataFrame型と呼ばれるやつ．

```pycon
# d <- read.csv("data3a.csv")
>>> import pandas
>>> d = pandas.read_csv("data3a.csv")
>>> d
     y      x  f
0    6   8.31  C
1    6   9.44  C
2    6   9.50  C
3   12   9.07  C
...（中略）...
98   7  10.86  T
99   9   9.97  T

[100 rows x 3 columns]

# d$x
>>> d.x
0      8.31
1      9.44
2      9.50
...（中略）...
98    10.86
99     9.97
Name: x, Length: 100, dtype: float64


# d$y
>>> d.y
0      6
1      6
2      6
...（中略）...
98     7
99     9
Name: y, Length: 100, dtype: int64

# d$f
>>> d.f = d.f.astype('category')
>>> d.f
0     C
1     C
2     C
3     C
...（中略）...
98    T
99    T
Name: f, Length: 100, dtype: category
Categories (2, object): [C < T]
```

Rではf列のことを**因子**（factor）と呼ぶらしいけど，
python(pandas)ではcategoryになってる．

データの型を確認するには

```pycon
# class(d)はちょっと確認できず．
# class(d$y)やらclass(d$x)やらclass(d$f)
>>> d.y.dtype
dtype('int64')
>>> d.x.dtype
dtype('float64')
>>> d.f.dtype
dtype('O')
```

データフレームの概要は，Series（前回のと一緒）

```pycon
# summary(d)
>>> d.describe()
                y           x
count  100.000000  100.000000
mean     7.830000   10.089100
std      2.624881    1.008049
min      2.000000    7.190000
25%      6.000000    9.427500
50%      8.000000   10.155000
75%     10.000000   10.685000
max     15.000000   12.400000

>>> d.f.describe()
count     100
unique      2
top         T
freq       50
Name: f, dtype: object
```

### 可視化
Rみたいに一行で散布図を色分ける方法が見つからず．．．

```pycon
# plot(d$x, d$y, pch=c(21, 19)[d$f])
>>> plt.plot(d.x[d.f == 'C'], d.y[d.f == 'C'], 'bo')
[<matplotlib.lines.Line2D at 0x1184d0490>]
>>> plt.plot(d.x[d.f == 'T'], d.y[d.f == 'T'], 'ro')
[<matplotlib.lines.Line2D at 0x111b58dd0>]
>>> plt.show()
>>> d.boxplot(column='y', by='f')
<matplotlib.axes._subplots.AxesSubplot at 0x11891a9d0>
```

![Screen Shot 2015-06-01 at 18.53.49.png](https://qiita-image-store.s3.amazonaws.com/0/48856/5a921784-c781-7c16-af78-46d6d05447b7.png)
![Screen Shot 2015-06-01 at 18.51.47.png](https://qiita-image-store.s3.amazonaws.com/0/48856/5b74fdba-34d7-f2e4-56dd-3a123629f970.png)

これらの図からわかることは

- 体サイズ$x_i$が増加するにつれ，種子数$y$が増加しているように見える
- 肥料の効果$f$はないように思われる

# ポアソン回帰の統計モデル

カウントデータである種子数データをうまく表現できそうな統計モデルをつくる．

## 個体$i$の体サイズ$x_i$に依存する統計モデル

**説明変数**：$x_i$
**応対変数**：$y_i$

ある個体$i$において種子数が$y_i$である確率$p(y_i|\lambda_i)$はポアソン分布にしたがっていて，

```math
p(y_i|\lambda_i) = \frac{\lambda _i ^{y_i} \exp (-\lambda_i)}{y_i !}

```

と仮定する．
ここで，ある個体$i$の平均種子数$\lambda_i$が

```math
\begin{eqnarray}
\lambda_i &=& \exp(\beta_1 + \beta_2 x_i )
\log \lambda_i &=& \beta_1 + \beta_2 x_i 
\end{eqnarray}
```

であるとする．
$\beta_1$を**切片**，$\beta_2$は**傾き**とよび，
右辺$\beta_1 + \beta_2 x_i$を**線形予測子**とよぶ．

また（$\lambda_i$の関数）＝（線形予測子）となっている場合，
左辺の関数は**リンク関数**と呼ばれる．

## ポアソン回帰

ポアソン回帰は以下の式が最大となる$\beta_1, \beta_2$を見つけること．

```math

\log L(\beta_1, \beta_2) = \sum_i \log \frac{\lambda_i^{y_i}\exp(-\lambda_i)}{y_i !}
```

Rだと`glm`とかいうのを使えば一発らしい．
pythonだと`statsmodels`ライブラリを使って延々と頑張らないといけないとか，，，（[参考](http://breakbee.hatenablog.jp/entry/2014/05/13/114341)）

```pycon
# fit <- glm(y ~ x, data=d, familiy=poisson)
>>> import statsmodels.api as sm
>>> import statsmodel.formula.api as smf
>>> model = smf.glm('y ~ x', data=d family=sm.families.Poisson())
>>> results = model.fit()
>>> print(results.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:                 月, 01  6 2015   Deviance:                       84.993
Time:                        23:06:45   Pearson chi2:                     83.8
No. Iterations:                     7
==============================================================================
                 coef    std err          z      P>|z|      [95.0% Conf. Int.]
------------------------------------------------------------------------------
const          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
==============================================================================
>>> results.llf
-235.38625076986077
```

RではInterceptってなってるらしいけど，
Pythonでは一番下のconstが$\beta_1$（切片），xが$\beta_2$（傾き）を表しているらしい．

つまり，最尤推定量は$\beta_1=1.2917$と$\beta_2=0.0757$ということ．

std errは**標準誤差**を示していて，$\beta_1, \beta_2$のばらつきを示している．

> 同じ調査方法で同数の別データをサンプリングしなおすと，最尤推定量も変わり，そのばらつき具合

zは**z値**という統計量で，最尤推定量を標準誤差で割ったもの．

このzによって**Wald信頼区間**というものが推定できるらしい
この信頼区間はどうやら，最尤推定量が取りうる値に０が入らないということっぽい[参照](http://hosho.ees.hokudai.ac.jp/~kubo/log/2009/0101.html#08)




## ポアソン回帰モデルによる予測

ポアソン回帰の推定結果を使って，体サイズ$x$における平均種子数$\lambda$を予測する．

```math
\lambda = \exp(1.29 + 0.0757x)
```

を用いて図示する．

```pycon
# plot(d$x, d$y, pch=c(21, 19)[d$f])
# xx <- seq(min(d$x), max(d$x), length=100)
# lines(xx, exp(1.29 + 0.0757 * xx), lwd = 2)
>>> plt.plot(d.x[d.f == 'C'], d.y[d.f == 'C'], 'bo')
>>> plt.plot(d.x[d.f == 'T'], d.y[d.f == 'T'], 'ro')
>>> xx = [d.x.min() + i * ((d.x.max() - d.x.min()) / 100.) for i in range(100)]
>>> yy = [numpy.exp(1.29 + 0.0757 * i )for i in xx]
>>> plt.plot(xx, yy)
>>> plt.show()
```

![Screen Shot 2015-06-04 at 16.52.44.png](https://qiita-image-store.s3.amazonaws.com/0/48856/6bf9b567-132a-c223-d3a5-0ab220c40146.png)

# 3.5 説明変数が因子型の統計モデル

施肥処理$f_i$を説明変数として組み込んだ統計モデルを考える．

因子型の説明変数はGLMの中では，**ダミー変数**を用いる．
つまり，モデルの平均値を，

```math
\begin{eqnarray}
\lambda_i &=& \exp (\beta_1 + \beta_3 d_i) \\

 \therefore d_i &=& \left\{ \begin{array}{ll}
    0 & (f_i = Cの場合) \\
    1 & (f_i = Tの場合)
  \end{array} \right.

\end{eqnarray}
```

と書くことになる．

```pycon
# fit.f <- glm(y ~ f, data=d, family=poisson)
>>> model = smf.glm('y ~ f', data=d, familiy=sm.families.Poisson())
>>> model.fit().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:                -237.63
Date:                 木, 04  6 2015   Deviance:                       89.475
Time:                        17:20:11   Pearson chi2:                     87.1
No. Iterations:                     7
==============================================================================
                 coef    std err          z      P>|z|      [95.0% Conf. Int.]
------------------------------------------------------------------------------
Intercept      2.0516      0.051     40.463      0.000         1.952     2.151
f[T.T]         0.0128      0.071      0.179      0.858        -0.127     0.153
==============================================================================
>>> model.llf
-237.62725696068682
```


この結果から，もし個体$i$の$f_i$がCならば，

```math
\lambda_i = \exp(2.05 + 0.0128 * 0) = \exp(2.05) = 7.77
```

もしTならば，

```math
\lambda_i = \exp(2.05 + 0.0128 * 1) = \exp(2.0628) = 7.87
```

となり，**肥料を与えると平均種子数がほんの少しだけ増加する**と予測している．

対数尤度（llf）については，小さくなってしまっているため，あてはまりが悪いといえる．


# 3.6 説明変数が数量型 ＋ 因子型の統計モデル

個体の体サイズ$x_i$と施肥処理$f_i$の両方を組み込んだ統計モデルを考える．

```math
\log \lambda_i = \beta_1 + \beta_2 x_i + \beta_3 d_i

```

```pycon
>>> model = smf.glm('y ~ x + f', data=d, family=sm.families.Poisson())
>>> result = model.fit()
>>> result.summary()
                 Generalized Linear Model Regression Results
==============================================================================
Dep. Variable:                      y   No. Observations:                  100
Model:                            GLM   Df Residuals:                       97
Model Family:                 Poisson   Df Model:                            2
Link Function:                    log   Scale:                             1.0
Method:                          IRLS   Log-Likelihood:                -235.29
Date:                 木, 04  6 2015   Deviance:                       84.808
Time:                        17:31:34   Pearson chi2:                     83.8
No. Iterations:                     7
==============================================================================
                 coef    std err          z      P>|z|      [95.0% Conf. Int.]
------------------------------------------------------------------------------
Intercept      1.2631      0.370      3.417      0.001         0.539     1.988
f[T.T]        -0.0320      0.074     -0.430      0.667        -0.178     0.114
x              0.0801      0.037      2.162      0.031         0.007     0.153
==============================================================================
>>> result.llf
-235.29371924249369
```

最大対数尤度は最も大きいことから，上記２つの統計モデルよりあてはまりが良いといえる．

## 対数リンク関数のわかりやすさ：掛け算される効果

上記の数量型 ＋ 因子型の統計モデルの数式を変形すると，

```math
\begin{eqnarray}
\log \lambda_i &=& \beta_1 + \beta_2 x_i + \beta_3 d_i \\
\lambda_i &=& \exp(\beta_1 + \beta_2 x_i + \beta_3 d_i) \\
\lambda_i &=& \exp(\beta_1) * \exp(\beta_2 x_i) * \exp(\beta_3 d_i) \\
\lambda_i &=& \exp(定数) * \exp(体サイズの効果) * \exp(施肥処理の効果) 
\end{eqnarray}
```

と表せ，効果が掛け合わされているのがわかる．

> リンク関数がなにもない場合，**恒等リンク関数**と呼ぶ．

GLMにおいて

正規分布かつ，口頭リンク関数の場合を **線形モデル**(**linear model, LM**)もしくは**一般線形モデル**(**general linear model**)と呼ぶ．