R で個体差のあるロジスティック回帰1 - glmmML

R GLM Bayes

前回 のロジスティック回帰に続き、書籍 「 データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学) 」のサンプルを使って個体差を考慮したロジスティック回帰を GLMM と階層ベイズモデルで試してみます。

  • (1) GLMM によるロジスティック回帰 （glmmML 関数）
  • (2) MCMCglmm を使った階層ベイズモデルによるロジスティック回帰（MCMCglmm 関数）

ただし、今回は (1) だけで (2) に関しては次回に書く予定です。

サンプルソースは http://github.com/fits/try_samples/tree/master/blog/20140209/

はじめに

パッケージのインストール

今回使用するパッケージを R へインストールしておきます。

install.packages("glmmML")

データ

データは書籍のサポート web サイトから取得します。 今回は 7章「一般化線形混合モデル（GLMM）- 個体差のモデリング -」のサンプルデータを使います。

データ data7.csv
"N","y","x","id"
8,0,2,1
8,1,2,2
8,2,2,3
・・・

データ内容は下記の通りです。

項目 内容
N 調査種子数
y 生存種子数
x 植物の葉数
id 個体のID

葉数 { x_i } の値が 2 ～ 6 までとなっており、 葉数毎に 20 個体ずつの生存種子数 { y_i } をサンプリングしたデータとなっています。 （合計 100 個体）

二項分布の性質と過分散

ところで、二項分布の性質として下記があります。
今回のケースでは n が調査種子数（= 8）、p が生存確率です。

  • (1) 平均（期待値とするのが適切かも）が { np }
  • (2) 分散が { np(1 - p) }

ここで、今回のデータは下記のようになっています。

葉数 4（x = 4）場合の平均と分散
> mean(d[d$x == 4,]$y)
[1] 4.05

> var(d[d$x == 4,]$y)
[1] 8.365789

np = 4.05 とすると p = 4.05 / 8 = 約 0.5 となり、二項分布として期待される分散は np(1 - p) = 8 × 0.5 × (1 - 0.5) = 2 で、実際の分散 8.36 の方が明らかに大きい値となっています。 （過分散）

二項分布として考えると、原因不明な個体差等による効果を考慮した統計モデルを使ったロジスティック回帰が必要となり、GLMM や階層ベイズを使用するという流れとなるようです。

(1) GLMM によるロジスティック回帰 （glmmML 関数）

それでは glmmML() 関数を使った推定値の算出とグラフ化を実施してみます。

推定方法

個体差を考慮するため、GLMM では個体差のパラメータ { r_i } を追加した線形予測子 { logit(q_i) = \beta_1 + \beta_2 x_i + r_i } を使う事になるようです。
ここで、 { q_i } は生存確率です。

そして、{ r_i } が正規分布の確率分布に従っていると仮定し、
確率密度関数 { p(r_i | s) = \frac{1}{\sqrt{ 2 \pi s_^2} } \exp(- \frac{r_i ^2}{ 2s ^2 }) } を加味して { r_i } を積分した
尤度 { L_i = \int p(y_i | \beta_1, \beta_2, r_i) p(r_i | s) dr_i } に対して
全データの積 { L(\beta_1, \beta_2, s) = \prod_i L_i } の対数尤度 { \log L(\beta_1, \beta_2, s) } が最大になるようなパラメータ { \beta_1, \beta_2, s } の最尤推定値を求めます。

なお、{ p(y_i | \beta_1, \beta_2, r_i) } は二項分布の確率密度関数で、 { s } は個体差 { r_i } のばらつき（標準偏差）です。

glmmML() による推定

glmmML() 関数へ指定するオプションは、前回 の glm() とほとんど同じですが 、下記の点が異なります。

  • 個体差のパラメータ { \r_i } 部分を cluster オプションで指定

今回は個体毎に異なる値が設定されている id 列を cluster オプションへ指定します。

logiGlmmML.R
d <- read.csv('data7.csv')

library(glmmML)

d.res <- glmmML(cbind(y, N - y) ~ x, data = d, family = binomial, cluster = id)

summary(d.res)
・・・

glmmML() による推定結果は下記の通りです。

実行結果
Call:  glmmML(formula = cbind(y, N - y) ~ x, family = binomial, data = d, cluster = id) 


              coef se(coef)      z Pr(>|z|)
(Intercept) -4.190   0.8777 -4.774 1.81e-06
x            1.005   0.2075  4.843 1.28e-06

Scale parameter in mixing distribution:  2.408 gaussian 
Std. Error:                              0.2202 

        LR p-value for H_0: sigma = 0:  2.136e-55 

Residual deviance: 269.4 on 97 degrees of freedom   AIC: 275.4

{ \hat \beta_1 = -4.190, \hat \beta_2 = 1.005, \hat s = 2.408 } と推定されました。

なお、glmmML() の結果から、
{ \hat \beta_1, \hat \beta_2 } の値は $coefficients で、{ \hat s } の値は $sigma で取得できます。

葉数と生存種子数

推定したパラメータを使って葉数 x と生存種子数 y の予測線をグラフ化する処理は下記のようになります。

生存種子数の予測値は 調査種子数（= max(d$N) = 8）× 生存確率 で算出し、
生存確率 { q_i } は { logit(q_i) = \beta_1 + \beta_2 x_i + r_i } の式を元に算出します。

logiGlmmML.R
・・・
# 生存確率の算出
calcProb <- function(x, b, r)
    1.0 / (1.0 + exp(-1 * (b[1] + b[2] * x + r)))

png("logiGlmmML_1.png")

# x と y をプロット
plot(d$x, d$y)

# x の範囲（2 ～ 6）を 100分割
xx <- seq(min(d$x), max(d$x), length = 100)

beta <- d.res$coefficients

# 個体差 r = 0 の葉数と生存種子数との予測線
lines(xx, max(d$N) * calcProb(xx, beta, 0), col="green")
# 個体差 r = -2.408 の葉数と生存種子数との予測線
lines(xx, max(d$N) * calcProb(xx, beta, -1 * d.res$sigma), col="blue")
# 個体差 r = 2.408 の葉数と生存種子数との予測線
lines(xx, max(d$N) * calcProb(xx, beta, d.res$sigma), col="blue")

dev.off()
・・・

個体差 r = 0 （個体差なし）とした場合の予測線を緑色で、 r = ±2.408 の予測線を青色で描画しています。

f:id:fits:20140209172056p:plain

葉数 { x_i = 4 } の生存種子数分布

次に、葉数 x が 4 の場合の生存種子数 y と個体数の分布と予測線をグラフ化する処理は下記のようになります。

glmmML() の結果から、x = 4 の場合の y の確率分布を算出し、x = 4 のサンプル数（= 20）を乗算して予測線を描画しています。

y の確率分布は、0 ～ 8 のそれぞれの値に対して { \int p(y_i | \beta_1, \beta_2, r_i) p(r_i | s) dr_i } で算出します。 （二項分布 × 正規分布 を個体差 { r_i } で積分）

下記では sapply() と integrate() 関数を使って上記の算出処理を実装しています。

なお、lower と upper オプションを使って積分範囲を { -10 \times \hat s } から { 10 \times \hat s } としています。

logiGlmmML.R
・・・
png("logiGlmmML_2.png")

# 0 ～ 8
yy <- 0:max(d$N)

# 生存種子数 y と個体数の分布
plot(yy, table(d[d$x == 4,]$y), xlab="y", ylab="num")

# 葉数 x を固定した場合の生存種子数 y の確率分布を算出
calcL <- function(ylist, xfix, n, b, s)
  sapply(ylist, function(y) integrate(
      f = function(r) dbinom(y, n, calcProb(xfix, b, r)) * dnorm(r, 0, s),
      lower = s * -10,
      upper = s * 10
    )$value
  )

# 予測線
lines(yy, calcL(yy, 4, max(d$N), beta, d.res$sigma) * length(d[d$x == 4,]$y), col="red", type="b")

dev.off()

f:id:fits:20140209172110p:plain

今回はここまで。 MCMCglmm を使った階層ベイズは次回。