GLMの基礎：二項分布を使ったGLM

「ロジスティック回帰分析」からの続き。

以下は、前回示した生存種子数（生存率）と「体サイズ、肥料の有無」の分析結果（ 本書図 6.7 に相当）。

二項分布とポアソン分布は、このブログでも何度か取り上げているが、本書の解説にそって定義すると

二項分布は上限のある「あり／なし」カウントデータを表現、例えば「N 回のコイン投げで表の回数が y, 裏は N - y 」。一方、ポアソン分布は上限のないカウントデータを表現。N = 1 かつ y ∈ {0, 1} の場合は、ベルヌーイ分布（Bernoulli distribution）。

よって、今回の観測データ「種子数 8 個のうちの生存数」を表現する統計モデルでは二項分布を使う。

右図は二項分布の一般式。

分析結果の図 6.7 のS字曲線は「ロジスティック曲線」で表現できる（参照：ロジスティック回帰分析）。二項分布を使った GLM の一種が、ロジスティック回帰分析で

確率分布：二項分布、リンク関数：ロジットリンク関数

を指定する。

ロジスティック関数は

> logistic <- function(z) {1/(1+exp(-z))}
> z<-seq(-6,6,0.1)
> ggplot(data.frame(x=z,y=logistic(z)),aes(x=x,y=y))+geom_line()+geom_vline(xintercept = 0,linetype="longdash")+ylab("確率 q")+xlab("線形予測子 z")


exp(0) = 1 なので、z = 0 で q = 0.5 の曲線となる。

線形予測子は z_i = β₁ + β₂x_i +...

ロジスティック関数を以下のように定義して

> logistic <- function(z) {1/(1+exp(-z))}

y 軸の確率 q  が、体サイズ x_i だけに依存すると仮定した曲線例は、左図：傾き一定（β₂=2）、右図：切片一定（β₁=0）で描画（図 6.5 に相当、末尾「Rコード」参照）

左図：logistic(0 + 2x), logistic(2 + 2x), logistic(-3 + 2x)
右図：logistic(0 + 4x), logistic(0 + 2x), logistic(0 - 1x)

線形予測子は z_i = β₁ + β₂x_i で、q と x_i の関係が、パラメーターβ₁, β₂ に依存している。

ロジスティック関数と線形予測子、ロジスティック関数を変形したロジット関数（logit function）は以下とおり。

ロジット関数はロジスティック関数の逆関数で、オッヅ q_i / (1-q_i) を対数変換したもの。以上が、確率分布、線形予測子、リンク関数。

次に、この統計モデルにデータをあてはめてパラメーター推定する。

尤度関数と対数尤度関数は

q_i は {β_j} の関数なので、logL を最大にする {β_j} を探すのが最尤推定。

具体的な分析例は前回の「ロジスティック回帰分析」を参照。


「オフセット項（割算値の回避）」に続く。


Rコード

図 6.5 の描画スクリプト。

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

> logistic <- function(z) {1/(1+exp(-z))}

> x<-seq(-3,3,0.1)

> p1<-ggplot(data.frame(x=x,y=logistic(0+2*x)),aes(x=x,y=y))+geom_line()+geom_vline(xintercept = 0,linetype="longdash")+ylab("確率 q")+xlab("説明変数 x")
> p2<-p1+geom_line(data = data.frame(x=x,y=logistic(2+2*x)),aes(x=x,y=y),color="red")
> p3<-p2+geom_line(data = data.frame(x=x,y=logistic(-3+2*x)),aes(x=x,y=y),color="blue")
> pA<-p3+annotate("text",label="β1=0",x=0.3,y=0.5)+annotate("text",label="β1=2",x=-0.3,y=0.9,color="red")+annotate("text",label="β1=-3",x=1.2,y=0.25,color="blue")
> p1<-ggplot(data.frame(x=x,y=logistic(0+4*x)),aes(x=x,y=y))+geom_line()+geom_vline(xintercept = 0,linetype="longdash")+ylab("確率 q")+xlab("説明変数 x")
> p2<-p1+geom_line(data = data.frame(x=x,y=logistic(0+2*x)),aes(x=x,y=y),color="red")
> p3<-p2+geom_line(data = data.frame(x=x,y=logistic(0-1*x)),aes(x=x,y=y),color="blue")
> pB<-p3+annotate("text",label="β2=4",x=0.5,y=1.0)+annotate("text",label="β2=2",x=-1,y=0.25,color="red")+annotate("text",label="β2=-1",x=1.7,y=0.25,color="blue")

> grid.newpage()

> pushViewport(viewport(layout = grid.layout(1,2)))

> print(pA,vp=viewport(layout.pos.row=1,layout.pos.col = 1))

> print(pB,vp=viewport(layout.pos.row=1,layout.pos.col = 2))

/**Tex **/