GLMの基礎：GLMの詳細

「GLMの目的と特徴」からの続き。

ここでは一般化線形モデル（GLM）の詳細を取り上げるが、これ以上は本書を参照。

「正規分布、直線回帰の限界」に記したように

GLM は 確率分布、リンク関数、線形予測子を指定する統計モデル

この 3 つの関係を GLM では理解すること。


確率分布

今回の 確率分布 はポアソン分布で

これは、左の第 2 章のモデルと、添え字 i を除いて同じ。

個体ごとの平均種子数 λi と、観測データである「応答変数 y_i」の関係 p(y_i | λi) は、λ_i の場合の y_i の条件付き確率で、右辺のポアソン分布に従うとしたモデル。


リンク関数、線形予測子

λ_i を、今回は植物のサイズである「説明変数 x_i 」の関数で定義すると

λ_i = exp(β₁ + β₂x_i)

この式を、β₁= {-1,-2,-1,-1}, β₂ = {0.4, 0.4, 0.5, -0.4} でプロットしたのが以下（本書の図3.4に相当）で、β₁ が切片、β2 が傾きに相当するのが分かる。

> df1 <- data.frame(type="1", x= x,y=exp(-1 + 0.4*x))
> df2 <- data.frame(type="2", x= x,y=exp(-2 + 0.4*x))
> df3 <- data.frame(type="3", x= x,y=exp(-1 + 0.5*x))
> df4 <- data.frame(type="4", x= x,y=exp(-1 + -0.4*x))
> df <- rbind(df1,df2,df3,df4)
> ggplot(df,aes(x,y,color=type))+geom_line()+ylim(0,4.5)+xlim(-4,5)

この式の両辺を対数変換

log λ_i = β₁ + β₂x_i

右辺 β₁ + β₂x_i が 線形予測子：二次関数のような曲線でも線形予測子 

左辺 log λ_i が リンク関数：今回は対数関数で「対数リンク関数（log link function）」と呼ぶ

図3.4 からも分かるように λ_i  = exp(線形予測子) ≥ 0 で、ポアソン分布のパラメータ λ ≥ 0 の条件と一致する、これがポアソン回帰の GLM で対数リンク関数を使う理由。

なお、第 6 章「GLM の応用範囲をひろげる」では（P.120）、確率分布は二項分布、リンク関数はロジットリンク関数（logit link function）。二項分布のパラメータは確率 q なので 0 ≤ q ≤ 1で、ロジットリンク関数はこの q の制約と、線形予測子をうまく関連づけるリンク関数となっている。


パラメーター最尤推定

ポアソン回帰とは、観測データに対するポアソン分布を使った統計モデルのあてはめ（fitting）であり、この統計モデルの対数尤度 logL が最大になるパラメータ β₁ と β₂ の推定値を求めること P.49

このモデルの対数尤度は

線形予測子 log λ_i = β₁ + β₂x_i から、λ_i が β₁, β₂ の関数であることに注意。

「パラメータの最尤推定法」は本書第 2 章の通りだが、ここで異なるのは {β₁, β₂} と複数パラメータの最尤推定値を求めること。ここでは、手動ではなく glm 関数で求める。

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

注意：family = possion(link = "log") が初期値のため、対数リンク関数の指定は省略可能。

summary(fit) より、パラメーター推定値は

Coefficients:

            Estimate Std. Error z value Pr(>|z|)    

(Intercept)  1.29172    0.36369   3.552 0.000383 ***

x            0.07566    0.03560   2.125 0.033580 * 

注意：「パラメーターの標準誤差の推定値」の「Std. Error」, z value, pr(>|z|) については割愛。今後「モデル選択」等を扱う際に、必要に応じて追記予定。

結果分析

パラメーター推定値を式に代入して、図示する。

λ = exp(1.29 + 0.0757x)

> x<-seq(min(d$x),max(d$x),length=100)
> df<-data.frame(x=x,y=exp(1.29+0.075*x))
> ggplot()+geom_point(data=d,aes(x=x,y=y,shape=f,colour=f))+geom_line(data=df,aes(x,y))

大きな傾向ではないが、「体サイズが大きいほど、平均種子の数が多い」が見て取れる。

同様に説明変数を「肥料の有無」にした場合

 λ_i = exp(β₁ + β₃d_i)

> fit.f<-glm(y~f,data=d,family=poisson())
> summary(fit.f)

Coefficients:

(Intercept)            x  

    1.29172      0.07566  

λ_i を C（肥料なし）、T（肥料あり）で求めると

> exp(as.numeric(fit.f$coefficients[1]))
[1] 7.78
> exp(as.numeric(fit.f$coefficients[1]) + as.numeric(fit.f$coefficients[2]))
[1] 7.88

わずかだが「肥料ありが平均種子数が多い」。

体サイズとの最大対数尤度の差は
> logLik(fit) - logLik(fit.f)

'log Lik.' 2.241006 (df=2)

「体サイズ」だけのモデルよりも「肥料有無」の方が小さく、当てはまりは悪くなっている。しかし、このような数値だけの「検定」で、確率モデルを評価できないのは「GLMの目的と特徴」に記した通り。


「逸脱度、残差逸脱度、逸脱残差」に続く。

/*tex */