ポアソン回帰とモデル選択 - データサイエンス時代で活躍できる人材になるために

新年明けましておめでとうございます．
前回G検定対策の記事を書いて以降，論文を書いていたりなので全然更新ができていませんでした．G検定は無事合格したので，今後は更新頻度を上げて，E資格のための記事も書いていこうと思います．

今回は下記の本の3章について取り上げたいと思います．

データ解析のための統計モデリング入門一般化線形モデル・階層ベイズモデル・MCMC （シリーズ確率と情報の科学） [ 久保拓弥 ]

ジャンル: 本・雑誌・コミック > 科学・医学・技術 > 数学
ショップ: 楽天ブックス
価格: 4,180円

統計モデリング入門

植物100個体の種子数のデータを用います. データの詳細として下記に６個体のデータを示します.
また，
y：種子の数
x：個体の大きさ,
f：肥料を与えた場合T，肥料を与えてない場合にC
です．

> d<-read.csv("data.csv")

   y     x f
1  6  8.31 C
2  6  9.44 C
3  6  9.50 C
4 12  9.07 C
5 10 10.16 C
6  4  8.32 C

個体の大きさや肥料の影響は種子数の数に影響あるのかな??

またどのような統計モデルを当てはめれば，データの特徴をより表すことができているかについて考えます．
f:id:koki12070930:20191118223755p:plain

個体の大きさが影響する場合

データの性質上，ポアソン分布を仮定します．
ある個体iの平均種子数を $\lambda_i$ とします．
ポアソン分布のパラメータ $\lambda_i$ は
$\lambda_i=\mathrm{exp}(\beta_1+\beta_2x_i)$
となります．
ただし， $x_i$ は個体 $i$ の大きさです．
$x_i$ の項があることから, 個体の大きさを考慮してることになっています. パラメータ $\beta_2$ の値が大きければ, 平均種子数は大きくなるように $\beta_1, \beta_2$ の値はとても重要であることが分かります. これらパラメータ $\beta_1, \beta_2$ を求める代表的な方法は最尤法です. Rではこの計算を自動でやってくれるパッケージがあるので今回はそれを用います!

> fit1<-glm(y~x,data=d,family=poisson(link="log"))
> fit1

Call:  glm(formula = y ~ x, family = poisson(link = "log"), data = d)

Coefficients:
(Intercept)            x  
    1.29172      0.07566  

Degrees of Freedom: 99 Total (i.e. Null);  98 Residual
Null Deviance:	    89.51 
Residual Deviance: 84.99 	AIC: 474.8

したがって,
$\lambda_i=\mathrm{exp}(1.29172+0.07566 x_i)$ と求めることができます.
このパラメータを用いてポアソン分布をを可視化します．
f:id:koki12070930:20200102001122p:plain

施肥処理が影響する場合

因子型説明変数はダミー変数に変換します．
\begin{eqnarray}
d_i&=&
\begin{cases}
0~(f_i=C)\\
1~(f_i=T)
\end{cases}
\end{eqnarray}
$\lambda_i=\mathrm{exp}(\beta_1+\beta_3d_i)$ であるので， $\beta_1$ ， $\beta_3$ のパラメータ推定をします．

> fit2<-glm(y~f,data=d,family=poisson)
> fit2

Call:  glm(formula = y ~ f, family = poisson, data = d)

Coefficients:
(Intercept)           fT  
    2.05156      0.01277  

Degrees of Freedom: 99 Total (i.e. Null);  98 Residual
Null Deviance:	    89.51 
Residual Deviance: 89.48 	AIC: 479.3

したがって,
肥料なし $\to \lambda_i=\mathrm{exp}(2.05156+0)$
肥料あり $\to \lambda_i=\mathrm{exp}(2.05156+ 0.01277 \times 1)$ 　

個体の大きさと施肥処理が影響する場合

ここでは $\lambda_i=\mathrm{exp}(\beta_1+\beta_2x_i+\beta_3d_i)$ です．

>fit3<-glm(y~x+f,data=d,family=poisson)
> fit3

Call:  glm(formula = y ~ x + f, family = poisson, data = d)

Coefficients:
(Intercept)            x           fT  
    1.26311      0.08007     -0.03200  

Degrees of Freedom: 99 Total (i.e. Null);  97 Residual
Null Deviance:	    89.51 
Residual Deviance: 84.81 	AIC: 476.6

したがって,
肥料なし $\to \lambda_i=\mathrm{exp}(1.26311+0.08007x_i)$
肥料あり $\to \lambda_i=\mathrm{exp}(1.26311+0.08007x_i -0.032)$ 　
となります．

先ほど同様，これらを用いてポアソン分布を可視化します．
f:id:koki12070930:20200102003328p:plain

当てはまりが良いモデルはどれ??

ここまで，様々な状況下でパラメータ推定を行なって来ました．
結局どのモデルが一番当てはまりが良いのでしょうか??

最尤法は，対数尤度の最大化問題を考えてパラメータ推定を行うのでした．
では問題です．

これまでのモデルの中で対数尤度が最大のものが一番良いモデルであると結論づけて良いでしょうか??

答えは"良くない"です．
モデルの選択の基準として，最大対数尤度を用いることは適切ではなく，AICを用いたモデル選択が重要になります．
なぜ最大対数では適さないのかは下記の記事を参考にしてください．
koki12070930.hatenablog.com

今回はモデルの選択には情報量基準AICを推奨しましたが，目的によっては尤度比検定を行うことがあります．
ポアソン分布以外のGLMの話題や混合モデルへの拡張について，次回以降に書きたいと思います!

ではまた!