RでGLMあれこれ（ポアソン回帰）

glm()関数を使ってモデリング

データセット：http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/poisson/data3a.csv

0. 統計モデリングの流れ

1. データの確認とモデルの選択
　要約統計量とグラフから、どの統計モデルを使うべきか決める。
2. 統計モデルの作成：あてはめ（fitting）
　GLMであれば、「説明変数」は、いくつかの「母数」によって説明される。
　「母数」は、リンク関数（「母数」の関数）と線形予測子（「応答変数」とその「パラメータ」によ
　　る線形結合）の関係によって推定される。
　　つまり、リンク関数・線形予測子を決定することで、モデリングが可能に。
3. 統計モデルの検証
　3-1. モデルによって「説明」が可能か？（あてはまりの良さ）
　3-2. モデルによって「予測」が可能か？（偶然性の小ささ）

1. データの確認とモデルの選択

> d   <- read.csv("data3a.csv")#ファイル読み込み
> head(d)#データフレームの確認
   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
> summary(d)#f列はCとTの2水準で構成
       y               x          f     
 Min.   : 2.00   Min.   : 7.190   C:50  
 1st Qu.: 6.00   1st Qu.: 9.428   T:50  
 Median : 8.00   Median :10.155         
 Mean   : 7.83   Mean   :10.089         
 3rd Qu.:10.00   3rd Qu.:10.685         
 Max.   :15.00   Max.   :12.400  
> class(d$y)#yはinteger型
[1] "integer"
> class(d$x)#xはinteger型
[1] "numeric"
> class(d$f)#fはfactor型
[1] "factor"
> plot(d$x, d$y, pch = c(21, 19)[d$f])#CとTで色分けして散布図
> legend("topleft", legend = c("C", "T"), pch = c(21, 19))#CとTの凡例を追加

f:id:yumaloop:20180320232343p:plain — 横軸にx列、縦軸にy列をとった散布図

　ポアソン分布に似ている。→ポアソン回帰でGLMが使えそう

この散布図に置いて、種子数yを説明変数、体サイズxと施肥処理fを応答変数として、統計モデル（ここではポアソン回帰）を作りたい。とりあえず、いくつかモデルを作りそれぞれの尤度をさぐる。

2. 統計モデルのあてはめ（fitting）

それぞれの仮定した統計モデルに対して、glm()関数で得られたパラメータ（ $\beta_1, \beta_2, \beta_3, ...$ ）の推定値により、母数（ $\lambda_i,$ ）を推定し、さらに得られた母数（ $\lambda_i$ ）から尤度（（ $L(\lambda_i)$ ））を求める。

※パラメータ（ $\beta_1, \beta_2, \beta_3, ...$ ）の推定値を求める過程
・ポアソン分布 $Po(\lambda)$ に従う確率変数 $x$ の確率密度関数 $f(x)$

　　 $f(x)={\large\frac{\lambda^{x}exp(-\lambda)}{x!}} {\scriptsize (x = 0, 1, 2, ...)}$

ある個体 $i$ において、種子数が $y_i$ である確率 $p(y_i | \lambda_i)$ はポアソン分布 $Po(\lambda_i)$ に従っていると仮定すると、以下の式が成り立つ。（ポアソン回帰の一般式）

　　 $p(y_i | \lambda_i)={\large\frac{\lambda_{i}^{y_i}exp(-\lambda_{i})}{y_i!}} {\scriptsize (y_i = 1, 2, 3, ...)}$

この式から、パラメータ（ $\beta_1, \beta_2, \beta_3, ...$ ）推定値は、ポアソン回帰においては以下の式によって導かれる。

　　 $\log{L(\beta_1, \beta_2, \beta_3, ...)} = \sum_{i=1}^{n} \log{\large{\frac{\lambda_{i}^{y_i}exp(-\lambda_{i})}{y_i!}}}$

すなわち、対数尤度関数 $\log{L(\beta_1, \beta_2, \beta_3, ...)}$ を最大にするようなパラメータ（ $\beta_1, \beta_2, \beta_3, ...$ ）の値を求めれば良い。

以下で異なる３種類のポアソン回帰モデルを立てて、推定結果を比較する。

（モデルA）「yはxによって説明される。」
　　　 ※各個体 $i$ について種子数 $x_i$ は、integer型で1以上の整数。
　　 $\lambda_i=exp(\beta_1+\beta_2x_i)$
　　i.e. $\log{\lambda_i}=\beta_1+\beta_2x_i$

> fit <- glm(y ~ x, data = d, family = poisson)#glm()関数でモデルの最尤推定値を得る
> fit#fitには推定結果が格納されている

Call:  glm(formula = y ~ x, family = poisson, 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

（モデルB）「yはfによって説明される。」
　　　 ※各個体 $i$ について施肥処理 $f_i$ は、factor型でC or Fをとる。
　　　 ※ $f_i$ のダミー変数 $d_i$ は、factor型で0 or 1をとる。

　　 $\lambda_i=exp(\beta_1+\beta_3d_i)$
　　i.e. $\log{\lambda_i}=\beta_1+\beta_3d_i$

> fit <- glm(y ~ x, data = d, family = poisson)
> fit

Call:  glm(formula = y ~ x, family = poisson, 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