GLMの基礎：GLMの目的と特徴

「正規分布、直接回帰の限界」からの続き。

本書の流れとは逆に、一般化線形モデル（GLM）の詳細の前に、本分析の目的や結果を先に示す。


背景

第 2 章「確率分布と統計モデルの最尤推定」のデータも、架空植物の種子数をポアソン分布で分析した。そこでは「λ は全個体で共通」だが、この第 3 章で扱うのは

個体ごとに λ が異なる、つまり個体ごとに異なる属性（説明変数）によって平均種子数が変化する統計モデル

そして

統計モデルを観測データにあてはめることをポアソン回帰（Poisson regression） 

同様の構造の統計モデルを総称して、一般化線形モデル（generalized linear model, GLM）

本分析で分かること

以下は、本章の観測データ

> d <- read.csv("http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/poisson/data3a.csv")

> summary(d)

       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         

x: 植物のサイズ、y: 採取した種子の数、f: 肥料を与えた／与えない「T/C」

> p1<-ggplot(d, aes(x=x,y=y,shape=f,colour=f))+geom_point()

> p2<-ggplot(d,aes(x=f,y=y))+geom_boxplot()
> multiplot(p1,p2,cols = 2) # For definiton, see Multiple graphs on one page (ggplot2)

この二つのプロットから、次のことが読み取れる

• サイズと種子数の関係性は不明（無相関に見える）
• 肥料を与えた効果もはっきりしない

本分析の主な結論である、体サイズ、肥料の有無と平均種子数の関係（予測）は、以下のように図示できる（本書図 3.8 の左図に相当）。

> fit.A <- glm(y ~ x+f,data=d,family = poisson(link = "log"))
> fit.A$coefficients
(Intercept)           x          fT 
 1.26310504  0.08007260 -0.03199939 
> x<-seq(5,20,by = 0.1)
> df1<-data.frame(x=x,y=exp(1.26+0.08*x),f="C")
> df2<-data.frame(x=x,y=exp(1.26+0.08*x-0.032),f="T")
> df<- rbind(df1,df2)

> theme_set(theme_gray(base_size=12, base_family="HiraKakuProN-W3"))

> ggplot(df,aes(x=x,y=y,colour=f))+geom_line()+xlab("体サイズ xi")+ylab("平均種子数 λi")

GLM の特徴であるリンク関数は λ_i = exp(β₁ + β₂x_i + β3d_i)  で、glm 関数で「切片と係数」 β₁,  β₂, β₃ を求める。

予測（or 分析結果）は先の図 3.8 と、以下の通り。

　肥料なし：λ_i = exp(1.26₁ + 0.08x_i)  = exp(1.26) × exp(0.08x_i)

　肥料あり：λ_i = exp(1.26₁ + 0.08x_i - 0.032) = exp(1.26) × exp(0.08x_i) × exp(-0.032)

• 体サイズが 1 増加すると、種子数は exp(0.08 × 1) = 1.08 倍と予測
• 肥料を与えると、種子数が exp(-0.032) = 0.969 倍と予測

前回同様に、直線回帰（恒等リンク関数 identity link function）は

> fit.B <- glm(y ~ x+f,data=d,family = poisson(link = "identity"))

> fit.B$coefficients

(Intercept)           x          fT 

  1.2669825   0.6606528  -0.2047486 

よって回帰式 λ_i = 1.27 + 0.661d_i  - 0.205d_i より

平均種子数が 0.1 個や 1,000個の場合、肥料を与えると 0.205個減り「-0.1個や999.8個になる」と予測

先の GLM の予測とは異なる理由は、対数リンク関数では、複数の効果が かけ算で影響 するため。ここでは GLM 予測の方が「妥当」。

なお、GML の場合と、この直線回帰（LM）との当てはまり具合に差（最大対数尤度の差）は、ほぼない
> logLik(fit.A); logLik(fit.B)
'log Lik.' -235.2937 (df=3)
'log Lik.' -235.159 (df=3)

つまり

どちらのリンク関数（今回は対数リンク関数、恒等リンク関数）が「妥当なモデル」かは、当てはまりの良しあしだけでは決まらない

「GLMの詳細」に続く。