ぺ ん ぎ ん の 閃 き

閃き- blog

きらびやかに、美しく、痛烈に.

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, ...)の値を求めれば良い。



以下で異なる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

・推定結果
β1の推定値 : 1.29172
β2の推定値 : 0.07566

(モデルC)「yはxとfによって説明される。」