GLMMを使って過分散に対処する

　以前紹介したが、仮にデータが厳密にポアソン分布や二項分布に従っていたとしても、データの内部に階層構造が存在するとき、過分散は生じえる。ゆえに、その階層構造をうまいこと考慮して解析を行うのが、一般化線形混合モデルGLMMだったわけである。本項では、GLMMを使って過分散なデータの解析を行ってみる。

具体的な解析：ポアソン分布

　では、以下のようなデータが得られているとき、どうやって解決していくかを確認してゆこう。group変数には、測定した各地点A,B,C……が格納されていると考えてほしい。

------------------------------------------------------

library(plotn)

library(glmmTMB)

x <- c(5.1, 23.8, 26.9, 16.1, 22.1, 27.6, 11.6, 9.2,

       20.2, 9.8, 26, 27.3, 13, 16.3, 25.3, 27.2, 

       23.6, 12.2, 10.4, 23.8, 9.7, 24, 29.4, 16.9,

       12.1, 26, 20.9, 20.1, 16.8, 12.6, 18.2, 18.1,

       16.4, 24.7, 6.6, 3.3, 3.6, 27.2, 18.2, 6.6)

y <- c(4, 6, 10, 5, 9, 6, 3, 6, 8, 6, 10, 9, 10, 10, 

       10, 12, 13, 10, 19, 9, 16, 7, 11, 16, 15, 13, 

       16, 11, 6, 10, 2, 7, 5, 1, 13, 7, 6, 5, 4, 8)

g <- 4

n <- 10

group <- rep(LETTERS[1:g], each = n)

d <- data.frame(x = x, group = group, y = y)

plotn(y ~ x, data = d)#図1の描画

------------------------------------------------------

0以上の離散値だから、ポアソン回帰を行いたくなるが、下記に示すように、このデータはポアソン分布にしては過分散である。

------------------------------------------------------

mean(y)

## [1] 8.85

var(y)

## [1] 17.15641

------------------------------------------------------

そこでgroupごとに色分けして、図示してみる。すると、group内、つまり同じ色内ではそれほどばらつきが大きいわけではないことに気が付く。

------------------------------------------------------

plotn(y ~ x, data = d, mode = "m", group = "group")#図2の描画

------------------------------------------------------

このような時が、まさに、階層構造によって引き起こされる過分散である。こういう時には、GLMMによっての解析が活きるときである。glmmTMBパッケージ内の、glmmTMB関数を使って解析してみよう。下記のようにコマンドの表記法はglm関数とほとんど変わらない。今回は、分布にポアソン分布、リンク関数は対数リンクを想定する。

------------------------------------------------------

fit1 <- glmmTMB(y ~ x + (1|group), data = d, family = poisson(link = "log"))

#一般化線形混合モデル

------------------------------------------------------

唯一、目新しい項目が、(1|group)という表記で、これがランダム効果に相当する。groupをランダム効果が生ずる要因だと考えており、この表記だと切片にのみランダム効果が存在することを仮定するモデルになっている。もし、(x1 + x2 + ……|group)とすれば、切片及び、含まれる固定効果の係数にランダム効果が付与されることを考慮したモデルになる（GLMMの項参照）。将来的には、いわゆる傾き（つまり固定効果）にランダム効果が付与されたモデルを検討するが、今回は切片だけにランダム効果を仮定しよう。

　では、解析結果を眺めてみる。

------------------------------------------------------

fit2 <- summary(fit1)

fit2

##  Family: poisson  ( log )

## Formula:          y ~ x + (1 | group)

## Data: d

## 

##      AIC      BIC   logLik deviance df.resid 

##    215.6    220.6   -104.8    209.6       37 

## 

## Random effects:

## 

## Conditional model:

##  Groups Name        Variance Std.Dev.

##  group  (Intercept) 0.1125   0.3354  

## Number of obs: 40, groups:  group, 4

## 

## Conditional model:

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

## (Intercept)  2.321892   0.220752  10.518   <2e-16 ***

## x           -0.010882   0.007732  -1.407    0.159    

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

------------------------------------------------------

今回、結果の出力で目新しいものは、Random effectsの項だろう。glmmTMBではランダム効果は線形予測子に正規分布の形で付与されることを想定し、Random effects項にはランダム効果をもたらす正規分布の分散（標準偏差）が推定されている。今回だと、groupランダム変数は切片に対して標準偏差0.34程度の正規分布から生成されている、と推定したようだ。

　今まで明示的に扱わなかったが、Conditional modelが固定効果の推定値である。切片は2.32、xの係数は-0.011程度のようだ。切片は有意だが、xの係数は有意ではないので、このデータはxとはあまり関係のないデータのようだ。このあたりの解釈は今までと変わらない。

　では、回帰曲線を描いてみよう。固定効果だけを使った回帰曲線は下記の通り。

------------------------------------------------------

b0e <- coef(fit2)$cond[1,1]

b1e <- coef(fit2)$cond[2,1]

plotn(y ~ x, data = d)#図3の描画

xx <- seq(min(d$x), max(d$x), length = 200)

yy <- exp(b0e + b1e * xx)

overdraw(points(xx, yy, type = "l"))

------------------------------------------------------

下記のようにglmmTMBパッケージなどでは、ranef関数が実装されており、ランダム効果の推定値を取り出すことができる（固定効果を取り出したいときはfixef関数もある）。

------------------------------------------------------

ranef(fit1)

## $group

##   (Intercept)

## A  -0.2609927

## B   0.2913757

## C   0.3483375

## D  -0.3589841

------------------------------------------------------

これを使うことで、各グループごとの回帰曲線も描くことができる。

------------------------------------------------------

plotn(y ~ x, data = d, mode = "m", group = "group")

for(i in 1:g){

  tmp <- d[d$group == unique(d$group)[i],]

  xx <- seq(min(tmp$x), max(tmp$x), length = 200)

  yy <- exp(b0e + b1e * xx + ranef(fit1)$cond$group[i,1])

  overdraw(points(xx, yy, type = "l", col = .default_col[i]))

}

------------------------------------------------------

ランダム効果によって切片がばらつくことで、過分散なデータになっていたことがわかるだろう。

　実は今回のデータは下記のようにして生成した。

------------------------------------------------------

b0 <- 2

b1 <- 0

s <- 1 #ランダム効果：正規分布の標準偏差

g <- 4

n <- 10

x <- round(runif(n*g, 0, 30), digits = 1)

group <- rep(LETTERS[1:g], each = n)

sds <- rep(rnorm(g, mean = 0, sd = s), each = n)

y <- rpois(n = n*g, lambda = exp(b0 + b1*x + sds))

------------------------------------------------------

つまり、ランダム効果の正規分布の標準偏差は1が正解だったわけであるが、今回は0.34と過少に推定されている。これはモデルが悪いとかではなくて、そもそもランダム効果のグループの数が少ないため、ばらついているに過ぎない（サンプルサイズ4で分散を推定するようなものである、後述するように、おそらくここで計算される分散は標本分散で、不偏分散ではないと思われる）。

　では、上記の設定のもとで、シミュレーションを行ってみよう。係数だけでなく、推定された分散も収集しておく。

------------------------------------------------------

b0 <- 2

b1 <- 0

s <- 1

g <- 4

n <- 10

p1 <- NULL

p2 <- NULL

b0e <- NULL

b1e <- NULL

s_e <- NULL

group <- rep(LETTERS[1:g], each = n)

for (i in 1:1000){

  x <- round(runif(n*g, 0, 30), digits = 1)

  sds <- rep(rnorm(g, mean = 0, sd = s), each = n)

  y <- rpois(n = n*g, lambda = exp(b0 + b1*x + sds))

  d <- data.frame(x = x, group = group, y = y)

  fit1 <- glmmTMB(y ~ x + (1|group), data = d, family = poisson(link = "log"))

  #一般化線形混合モデル

  fit2 <- summary(fit1)

  p1 <- c(p1, coef(fit2)$cond[1,4])

  p2 <- c(p2, coef(fit2)$cond[2,4])

  b0e <- c(b0e, coef(fit2)$cond[1,1])

  b1e <- c(b1e, coef(fit2)$cond[2,1])

  s_e <- c(s_e, as.numeric(fit2$varcor$cond$group)) #分散

}

histn(p1, xlab = "P value", ylab = "頻度")#図4の描画

histn(p2, xlab = "P value", ylab = "頻度")#図5の描画

sum(p1 < 0.05)

## [1] 968

sum(p2 < 0.05)

## [1] 51

------------------------------------------------------

切片の検出力は100%、xの係数の危険率は5%程度となった。過分散なデータでは危険率が増加する傾向があるので、適切な水準に収まったと言えるだろう。

　係数の推定値を確認してみる。

------------------------------------------------------

histn(b0e, xlab = "b0", ylab = "頻度")#図6の描画

histn(b1e, xlab = "b1", ylab = "頻度")#図7の描画

mean(b0e)

## [1] 1.977354

mean(b1e)

## [1] 0.0002803408

------------------------------------------------------

切片の推定値の平均は約2、xの係数の推定値の平均は約0なので、適切に推定できている。

　さらに推定された分散を確認してみる。

------------------------------------------------------

histn(s_e, xlab = "b1", ylab = "頻度")#図8の描画

mean(s_e)

## [1] 0.7536311

------------------------------------------------------

分散の平均は0.75で、母分散1の3/4倍程度になっている。ゆえに、推定されている分散は標本分散である。このように、GLMMを使うことで、固定効果をランダム効果を識別し、過分散なデータを解析することができた。

　次は、下記のデータを解析してみる。グループごとに識別してプロットしてみよう。

------------------------------------------------------

x <- c(1, 19.3, 23.4, 13.9, 25.4, 19, 13.3, 1.8, 6.1, 

       18.6, 23.9, 8.4, 20.8, 14.9, 29.4, 5.3, 26.8, 

       1.4, 11.7, 22.9, 21.1, 2.3, 9.1, 28.3, 28.3, 

       25.8, 17.9, 28.3, 7.1, 8.2, 4, 19.6, 18.7, 9.9,

       7.4, 6.7, 24.4, 21.1, 20.4, 24.7)

y <- c(1, 6, 12, 8, 16, 14, 5, 2, 2, 9, 26, 5, 26, 15,

       44, 2, 51, 4, 11, 31, 11, 1, 4, 26, 20, 20, 8, 

       15, 1, 4, 1, 3, 1, 2, 0, 1, 4, 2, 1, 1)

g <- 4

n <- 10

group <- rep(LETTERS[1:g], each = n)

d <- data.frame(x = x, group = group, y = y)

plotn(y ~ x, data = d, mode = "m", group = "group")#図9の描画

------------------------------------------------------

全体的には大きくばらつくが、グループ内ではそれほどばらつかない。実際、GLMによる解析で、ピアソンのχ^2と自由度の比をとると、普通のGLMでは過分散である。

------------------------------------------------------

sum(residuals(glm(y ~ x, family = poisson(link = "log")), type = "pearson")^2)/length(y-2)

## [1] 4.667615

------------------------------------------------------

ゆえに、GLMMをしてみる価値があるだろう。GLMMを行って、固定効果の回帰曲線を描く。

------------------------------------------------------

fit1 <- glmmTMB(y ~ x + (1|group), data = d, family = poisson(link = "log")) #一般化線形混合モデル

fit2 <- summary(fit1)

fit2

##  Family: poisson  ( log )

## Formula:          y ~ x + (1 | group)

## Data: d

## 

##      AIC      BIC   logLik deviance df.resid 

##    194.1    199.2    -94.1    188.1       37 

## 

## Random effects:

## 

## Conditional model:

##  Groups Name        Variance Std.Dev.

##  group  (Intercept) 0.7222   0.8498  

## Number of obs: 40, groups:  group, 4

## 

## Conditional model:

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

## (Intercept)   0.1676     0.4606   0.364    0.716    

## x             0.0960     0.0075  12.799   <2e-16 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

b0e <- coef(fit2)$cond[1,1]

b1e <- coef(fit2)$cond[2,1]

plotn(y ~ x)#図10の描画

xx <- seq(min(x), max(x), length = 200)

yy <- exp(b0e + b1e * xx)

overdraw(points(xx, yy, type = "l"))

------------------------------------------------------

固定効果の切片は0.17、xの係数は0.096である。各グループごとの回帰曲線は下記の通り。

------------------------------------------------------

ranef(fit1)

## $group

##   (Intercept)

## A  -0.2609927

## B   0.2913757

## C   0.3483375

## D  -0.3589841

plotn(y ~ x, data = d, mode = "m", group = "group")#図11の描画

for(i in 1:g){

  tmp <- d[d$group == unique(d$group)[i],]

  xx <- seq(min(tmp$x), max(tmp$x), length = 200)

  yy <- exp(b0e + b1e * xx + ranef(fit1)$cond$group[i,1])

  overdraw(points(xx, yy, type = "l", col = .default_col[i]))

}

------------------------------------------------------

指数関数は線形予測子の切片項が変わると、回帰曲線の接線（傾き）もそれに合わせて大きく変化する。なので、見た目的にはあたかもxの係数にランダム効果があるように見えるが、実はそれはないことに注意しよう（二項回帰だとより顕著なのでさらに注意）。データにはよく合った解析になっていることはわかると思う。

　実は今回のデータは下記のようにして生成した。

------------------------------------------------------

b0 <- 0

b1 <- 0.1

s <- 1.5 #ランダム効果：正規分布の標準偏差

g <- 4

n <- 10

x <- round(runif(n*g, 0, 30), digits = 1)

group <- rep(LETTERS[1:g], each = n)

sds <- rep(rnorm(g, mean = 0, sd = s), each = n)

y <- rpois(n = n*g, lambda = exp(b0 + b1*x + sds))

------------------------------------------------------

つまり、ランダム効果の正規分布の標準偏差は1.5が正解である。

　では、上記の設定のもとで、シミュレーションを行ってみよう。係数だけでなく、推定された分散も収集しておく。

------------------------------------------------------

b0 <- 0

b1 <- 0.1

s <- 1.5

g <- 4

n <- 10

p1 <- NULL

p2 <- NULL

b0e <- NULL

b1e <- NULL

s_e <- NULL

group <- rep(LETTERS[1:g], each = n)

for (i in 1:1000){

  x <- round(runif(n*g, 0, 30), digits = 1)

  sds <- rep(rnorm(g, mean = 0, sd = s), each = n)

  y <- rpois(n = n*g, lambda = exp(b0 + b1*x + sds))

  d <- data.frame(x = x, group = group, y = y)

  fit1 <- glmmTMB(y ~ x + (1|group), data = d, family = poisson(link = "log"))

  #一般化線形混合モデル

  fit2 <- summary(fit1)

  p1 <- c(p1, coef(fit2)$cond[1,4])

  p2 <- c(p2, coef(fit2)$cond[2,4])

  b0e <- c(b0e, coef(fit2)$cond[1,1])

  b1e <- c(b1e, coef(fit2)$cond[2,1])

  s_e <- c(s_e, as.numeric(fit2$varcor$cond$group)) #分散

}

histn(p1, xlab = "P value", ylab = "頻度")#図12の描画

histn(p2, xlab = "P value", ylab = "頻度")#図13の描画

sum(p1 < 0.05)

## [1] 157

sum(p2 < 0.05)

## [1] 1000

------------------------------------------------------

切片の危険率は15%とかなり高い。切片の推定は、あてにしにくいことがわかる。

　係数の推定値を確認してみる。

------------------------------------------------------

histn(b0e, xlab = "b0", ylab = "頻度")#図14の描画

histn(b1e, xlab = "b1", ylab = "頻度")#図15の描画

mean(b0e)

## [1] -0.008622309

mean(b1e)

## [1] 0.1001333

------------------------------------------------------

切片の推定値の平均は約0、xの係数の推定値の平均は約0.1なので、推定値自体は正しいと言える。おそらくは、ランダム効果が切片にかかるため、それのばらつきを受けて切片は有意になりやすい傾向があるのだと思われる。

　さらに推定された分散を確認してみる。

------------------------------------------------------

histn(s_e, xlab = "b1", ylab = "頻度")#図16の描画

mean(s_e)

## [1] 1.751688

------------------------------------------------------

分散の平均は1.75で、母分散1.5^2 = 2.25の3/4倍程度になっている。ゆえに、推定されている分散は標本分散である。このように、GLMMを使うことで、固定効果をランダム効果を識別し、過分散なデータを解析することができた。

　次は、下記のデータを解析してみる。極めて極端な場合である。グループごとに識別してプロットしてみよう。

------------------------------------------------------

x <- c(26.3, 24.2, 14.4, 22.6, 25.7, 20, 1.9, 10.3, 22.3,

       17.1, 17.5, 9, 4.7, 29.4, 24.8, 18.7, 5.1, 24.4, 

       8.7, 24.3, 6, 15.5, 7.5, 26.3, 14.8, 20.1, 7.2, 

       10.2, 3.1, 28.4)

y <- c(77, 65, 135, 1, 2, 1, 18, 11, 8, 232, 241, 355,

       31, 6, 21, 0, 0, 0, 714, 337, 793, 1, 2, 2, 4, 

       2, 8, 22, 29, 9)

g <- 10

n <- 3

group <- rep(LETTERS[1:g], each = n)

d <- data.frame(x = x, group = group, y = y)

plotn(y ~ x, data = d, mode = "m", group = "group")#図17の描画

------------------------------------------------------

全体的には大きくばらつくが、グループ内ではそれほどばらつかない。実際、GLMによる解析で、ピアソンのχ^2と自由度の比をとると、普通のGLMでは過分散である。

------------------------------------------------------

sum(residuals(glm(y ~ x, family = poisson(link = "log")), type = "pearson")^2)/length(y-2)

## [1] 286.5077

------------------------------------------------------

ゆえに、GLMMをしてみる価値があるだろう。GLMMを行って、固定効果の回帰曲線を描く。

------------------------------------------------------

fit1 <- glmmTMB(y ~ x + (1|group), data = d, family = poisson(link = "log")) #一般化線形混合モデル

fit2 <- summary(fit1)

fit2

##  Family: poisson  ( log )

## Formula:          y ~ x + (1 | group)

## Data: d

## 

##      AIC      BIC   logLik deviance df.resid 

##    216.4    220.6   -105.2    210.4       27 

## 

## Random effects:

## 

## Conditional model:

##  Groups Name        Variance Std.Dev.

##  group  (Intercept) 5.97     2.443   

## Number of obs: 30, groups:  group, 10

## 

## Conditional model:

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

## (Intercept)  3.243542   0.790976   4.101 4.12e-05 ***

## x           -0.047476   0.002899 -16.376  < 2e-16 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

b0e <- coef(fit2)$cond[1,1]

b1e <- coef(fit2)$cond[2,1]

plotn(y ~ x)#図18の描画

xx <- seq(min(x), max(x), length = 200)

yy <- exp(b0e + b1e * xx)

overdraw(points(xx, yy, type = "l"))

------------------------------------------------------

固定効果の切片は3.24、xの係数は-0.047である。一部の外れ値に引きずられることなく、推定がなされていることは注目に値するだろう。各グループごとの回帰曲線は下記の通り。

------------------------------------------------------

ranef(fit1)

## $group

##   (Intercept)

## A   2.2759582

## B  -1.8080876

## C  -0.2594625

## D   3.0482774

## E   0.5110935

## F  -4.0452131

## G   3.7281331

## H  -1.9535282

## I  -1.0561875

## J   0.2966789

plotn(y ~ x, data = d, mode = "m", group = "group")#図19の描画

for(i in 1:g){

  tmp <- d[d$group == unique(d$group)[i],]

  xx <- seq(min(tmp$x), max(tmp$x), length = 200)

  yy <- exp(b0e + b1e * xx + ranef(fit1)$cond$group[i,1])

  overdraw(points(xx, yy, type = "l", col = .default_col[i]))

}

------------------------------------------------------

また、各グループには3つしかデータ点がない。それでも、各グループの推定はうまくいっていることがわかる。

具体的な解析：二項分布

　　次は、二項分布に従うような下記のデータを解析してみる。グループごとに識別してプロットしてみよう。

------------------------------------------------------

x <- c(9.1, 11, 17.4, 11.3, 28.1, 28.2, 18.2, 4.6,

       14.9, 19.2, 13.2, 12.6, 7.6, 5.3, 19.4, 13,

       4.8, 14.3, 12.9, 18.6, 0.7, 9.8, 8, 3.7, 22.9, 

       23.8, 20.7, 29.1, 0.9, 1.1, 8.1, 12.1, 16.9, 

       20, 12.5, 25.1, 10, 0.8, 2.8, 23.9)

n <- c(39, 26, 30, 27, 24, 38, 39, 33, 35, 26, 21, 

       29, 27, 25, 38, 23, 24, 33, 32, 34, 32, 30,

       40, 29, 31, 20, 24, 37, 24, 31, 40, 39, 25,

       24, 32, 22, 34, 34, 26, 34)

y <- c(36, 25, 27, 25, 22, 38, 36, 32, 31, 26, 16,

       20, 20, 20, 31, 16, 17, 23, 25, 23, 29, 25,

       32, 25, 24, 19, 20, 29, 22, 24, 16, 16, 10,

       12, 19, 8, 15, 10, 13, 20)

g <- 4

n <- 10

group <- rep(LETTERS[1:g], each = n)

d <- data.frame(x = x, n = n, group = group, y = y, rate = y/n)

plotn(rate ~ x, data = d, mode = "m", group = "group")#図20の描画

------------------------------------------------------

全体的には大きくばらつくが、グループ内ではそれほどばらつかない。実際、GLMによる解析で、ピアソンのχ^2と自由度の比をとると、普通のGLMでは過分散である。

------------------------------------------------------

sum(residuals(glm(cbind(y,n-y) ~ 1, family = binomial(link = "logit")), type = "pearson")^2)/length(y-1)

## [1] 6.063804

------------------------------------------------------

ゆえに、GLMMをしてみる価値があるだろう。GLMMを行って、固定効果の回帰曲線を描く。

------------------------------------------------------

fit1 <- glmmTMB(cbind(y,n-y) ~ x + (1|group), data = d, family = binomial(link = "logit")) #一般化線形混合モデル

fit2 <- summary(fit1)

fit2

##  Family: binomial  ( logit )

## Formula:          cbind(y, n - y) ~ x + (1 | group)

## Data: d

## 

##      AIC      BIC   logLik deviance df.resid 

##    185.4    190.4    -89.7    179.4       37 

## 

## Random effects:

## 

## Conditional model:

##  Groups Name        Variance Std.Dev.

##  group  (Intercept) 1.107    1.052   

## Number of obs: 40, groups:  group, 4

## 

## Conditional model:

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

## (Intercept) 1.224964   0.547102   2.239   0.0252 *

## x           0.005371   0.009496   0.566   0.5717  

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

b0e <- coef(fit2)$cond[1,1]

b1e <- coef(fit2)$cond[2,1]

plotn(rate ~ x, data = d)#図21の描画

xx <- seq(min(x), max(x), length = 200)

yy <- 1/(1+exp(-b0e-b1e * xx))

overdraw(points(xx, yy, type = "l"))

------------------------------------------------------

固定効果の切片は1.2、xの係数は0.0054である。各グループごとの回帰曲線は下記の通り。

------------------------------------------------------

ranef(fit1)

## $group

##   (Intercept)

## A   1.3745488

## B  -0.2547913

## C   0.3303154

## D  -1.4824035

plotn(rate ~ x, data = d, mode = "m", group = "group")#図22の描画

for(i in 1:g){

  tmp <- d[d$group == unique(d$group)[i],]

  xx <- seq(min(tmp$x), max(tmp$x), length = 200)

  yy <- 1/(1+exp(-b0e-b1e * xx-ranef(fit1)$cond$group[i,1]))

  overdraw(points(xx, yy, type = "l", col = .default_col[i]))

}

------------------------------------------------------

ランダム効果によって切片がばらつくことで、過分散なデータになっていたことがわかるだろう。

　実は今回のデータは下記のようにして生成した。

------------------------------------------------------

b0 <- 2

b1 <- 0

s <- 2 #ランダム効果：正規分布の標準偏差

g <- 4

n <- 10

x <- round(runif(m*g, 0, 30), digits = 1)

n <- sample(x = 20:40, size = m*g, replace = T)

group <- rep(LETTERS[1:g], each = m)

sds <- rep(rnorm(g, mean = 0, sd = s), each = m)

y <- rbinom(n = m*g, size = n, prob = 1/(1+exp(-b0-b1*x-sds)))

------------------------------------------------------

つまり、ランダム効果の正規分布の標準偏差は2が正解だったわけであるが、今回は1.05と過少に推定されている。

　では、上記の設定のもとで、シミュレーションを行ってみよう。係数だけでなく、推定された分散も収集しておく。

------------------------------------------------------

b0 <- 2

b1 <- 0

s <- 2

g <- 4

n <- 10

p1 <- NULL

p2 <- NULL

b0e <- NULL

b1e <- NULL

s_e <- NULL

group <- rep(LETTERS[1:g], each = n)

for (i in 1:1000){

  x <- round(runif(m*g, 0, 30), digits = 1)

  n <- sample(x = 20:40, size = m*g, replace = T)

  sds <- rep(rnorm(g, mean = 0, sd = s), each = m)

  y <- rbinom(n = m*g, size = n, prob = 1/(1+exp(-b0-b1*x-sds)))

  d <- data.frame(x = x, n = n, group = group, y = y)

  fit1 <- glmmTMB(cbind(y,n-y) ~ x + (1|group), data = d, family = binomial(link = "logit"))

  #一般化線形混合モデル

  fit2 <- summary(fit1)

  p1 <- c(p1, coef(fit2)$cond[1,4])

  p2 <- c(p2, coef(fit2)$cond[2,4])

  b0e <- c(b0e, coef(fit2)$cond[1,1])

  b1e <- c(b1e, coef(fit2)$cond[2,1])

  s_e <- c(s_e, as.numeric(fit2$varcor$cond$group))

}

histn(p1, xlab = "P value", ylab = "頻度")#図23の描画

histn(p2, xlab = "P value", ylab = "頻度")#図24の描画

sum(p1 < 0.05)

## [1] 611

sum(p2 < 0.05)

## [1] 32

------------------------------------------------------

切片の検出力は61%、xの係数の危険率は3%程度となった。過分散なデータでは危険率が増加する傾向があるので、適切な水準に収まったと言えるだろう。

　係数の推定値を確認してみる。

------------------------------------------------------

histn(b0e, xlab = "b0", ylab = "頻度")#図25の描画

histn(b1e, xlab = "b1", ylab = "頻度")#図26の描画

mean(b0e)

## [1] 2.044522

mean(b1e)

## [1] 2.834862e-05

------------------------------------------------------

切片の推定値の平均は約2、xの係数の推定値の平均は約0なので、適切に推定できている。

　さらに推定された分散を確認してみる。

------------------------------------------------------

histn(s_e, xlab = "b1", ylab = "頻度")#図27の描画

mean(s_e)

## [1] 3.270719

------------------------------------------------------

分散の平均は3.3で、母分散2^2の3/4倍程度になっている。ゆえに、推定されている分散は標本分散である。このように、GLMMを使うことで、固定効果をランダム効果を識別し、過分散なデータを解析することができた。

　次は、下記のデータを解析してみる。グループごとに識別してプロットしてみよう。

------------------------------------------------------

x <- c(21.6, 3.4, 22.1, 16.1, 5.9, 20.3, 24.6, 8.2,

       3.8, 8.4, 10.3, 0.9, 8, 26.3, 21.8, 13.8, 1.7,

       20.6, 15.1, 8.2, 27.7, 20.2, 16.5, 0.1, 21, 

       24.6, 15.1, 20.4, 28.4, 8.9, 8.7, 15.8, 5.5, 

       10, 25.8, 15.4, 13.2, 5.2, 12.2, 5.4)

n <- c(28, 28, 25, 29, 31, 38, 28, 38, 40, 36, 

       32, 32, 29, 25, 22, 40, 20, 40, 27, 36, 

       28, 23, 31, 29, 33, 29, 34, 40, 36, 34, 

       25, 20, 21, 35, 28, 39, 25, 39, 33, 24)

y <- c(20, 10, 15, 8, 14, 23, 18, 7, 8, 15, 23, 

       9, 14, 24, 16, 29, 8, 30, 19, 16, 28, 22,

       26, 15, 29, 28, 29, 38, 36, 25, 16, 18, 

       18, 29, 24, 32, 19, 22, 24, 16)

g <- 4

n <- 10

group <- rep(LETTERS[1:g], each = n)

d <- data.frame(x = x, group = group, y = y)

plotn(rate ~ x, data = d, mode = "m", group = "group")#図28の描画

------------------------------------------------------

全体的には大きくばらつくが、グループ内ではそれほどばらつかない。実際、GLMによる解析で、ピアソンのχ^2と自由度の比をとると、普通のGLMでは過分散である。

------------------------------------------------------

sum(residuals(glm(cbind(y,n-y) ~ x, family = binomial(link = "logit")), type = "pearson")^2)/length(y-2)

## [1] 4.286451

------------------------------------------------------

ゆえに、GLMMをしてみる価値があるだろう。GLMMを行って、固定効果の回帰曲線を描く。

------------------------------------------------------

fit1 <- glmmTMB(cbind(y,n-y) ~ x + (1|group), data = d, family = binomial(link = "logit")) #一般化線形混合モデル

fit2 <- summary(fit1)

fit2

##  Family: binomial  ( logit )

## Formula:          cbind(y, n - y) ~ x + (1 | group)

## Data: d

## 

##      AIC      BIC   logLik deviance df.resid 

##    203.8    208.9    -98.9    197.8       37 

## 

## Random effects:

## 

## Conditional model:

##  Groups Name        Variance Std.Dev.

##  group  (Intercept) 0.5391   0.7342  

## Number of obs: 40, groups:  group, 4

## 

## Conditional model:

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

## (Intercept) -0.468731   0.391486  -1.197    0.231    

## x            0.098086   0.009555  10.265   <2e-16 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

b0e <- coef(fit2)$cond[1,1]

b1e <- coef(fit2)$cond[2,1]

plotn(rate ~ x, data = d)#図29の描画

xx <- seq(min(x), max(x), length = 200)

yy <- 1/(1+exp(-b0e-b1e * xx))

overdraw(points(xx, yy, type = "l"))

------------------------------------------------------

固定効果の切片は-0.47、xの係数は0.098である。各グループごとの回帰曲線は下記の通り。

------------------------------------------------------

ranef(fit1)

## $group

##   (Intercept)

## A  -1.0804197

## B  -0.2099821

## C   0.7703513

## D   0.5069592

plotn(rate ~ x, data = d, mode = "m", group = "group")#図30の描画

for(i in 1:g){

  tmp <- d[d$group == unique(d$group)[i],]

  xx <- seq(min(tmp$x), max(tmp$x), length = 200)

  yy <- 1/(1+exp(-b0e-b1e * xx-ranef(fit1)$cond$group[i,1]))

  overdraw(points(xx, yy, type = "l", col = .default_col[i]))

}

------------------------------------------------------

ロジスティック関数も同様に線形予測子の切片項が変わると、回帰曲線の接線（傾き）もそれに合わせて大きく変化する。なので、見た目的にはあたかもxの係数にランダム効果があるように見えるが、実はそれはないことに注意しよう。データにはよく合った解析になっていることはわかると思う。

　実は今回のデータは下記のようにして生成した。

------------------------------------------------------

b0 <- 0

b1 <- 0.1

s <- 1 #ランダム効果：正規分布の標準偏差

g <- 4

n <- 10

x <- round(runif(m*g, 0, 30), digits = 1)

n <- sample(x = 20:40, size = m*g, replace = T)

group <- rep(LETTERS[1:g], each = m)

sds <- rep(rnorm(g, mean = 0, sd = s), each = m)

y <- rbinom(n = m*g, size = n, prob = 1/(1+exp(-b0-b1*x-sds)))

------------------------------------------------------

つまり、ランダム効果の正規分布の標準偏差は1が正解である。

　では、上記の設定のもとで、シミュレーションを行ってみよう。係数だけでなく、推定された分散も収集しておく。

------------------------------------------------------

b0 <- 0

b1 <- 0.1

s <- 1

g <- 4

n <- 10

p1 <- NULL

p2 <- NULL

b0e <- NULL

b1e <- NULL

s_e <- NULL

group <- rep(LETTERS[1:g], each = n)

for (i in 1:1000){

  x <- round(runif(m*g, 0, 30), digits = 1)

  n <- sample(x = 20:40, size = m*g, replace = T)

  sds <- rep(rnorm(g, mean = 0, sd = s), each = m)

  y <- rbinom(n = m*g, size = n, prob = 1/(1+exp(-b0-b1*x-sds)))

  d <- data.frame(x = x, n = n, group = group, y = y)

  fit1 <- glmmTMB(cbind(y,n-y) ~ x + (1|group), data = d, family = binomial(link = "logit"))

  #一般化線形混合モデル

  fit2 <- summary(fit1)

  p1 <- c(p1, coef(fit2)$cond[1,4])

  p2 <- c(p2, coef(fit2)$cond[2,4])

  b0e <- c(b0e, coef(fit2)$cond[1,1])

  b1e <- c(b1e, coef(fit2)$cond[2,1])

  s_e <- c(s_e, as.numeric(fit2$varcor$cond$group))

}

histn(p1, xlab = "P value", ylab = "頻度")#図31の描画

histn(p2, xlab = "P value", ylab = "頻度")#図32の描画

sum(p1 < 0.05)

## [1] 156

sum(p2 < 0.05)

## [1] 1000

------------------------------------------------------

切片の危険率は15%とかなり高い。切片の推定は、あてにしにくいことがわかる。

　係数の推定値を確認してみる。

------------------------------------------------------

histn(b0e, xlab = "b0", ylab = "頻度")#図33の描画

histn(b1e, xlab = "b1", ylab = "頻度")#図34の描画

mean(b0e)

## [1] -0.02614901

mean(b1e)

## [1] 0.100223

------------------------------------------------------

切片の推定値の平均は約0、xの係数の推定値の平均は約0.1なので、推定値自体は正しいと言える。おそらくは、ランダム効果が切片にかかるため、それのばらつきを受けて切片は有意になりやすい傾向があるのだと思われる。

　さらに推定された分散を確認してみる。

------------------------------------------------------

histn(s_e, xlab = "b1", ylab = "頻度")#図35の描画

mean(s_e)

## [1] 0.7578497

------------------------------------------------------

分散の平均は0.75で、母分散1の3/4倍程度になっている。ゆえに、推定されている分散は標本分散である。このように、GLMMを使うことで、固定効果をランダム効果を識別し、過分散なデータを解析することができた。

　最後に、下記のデータを解析してみる。極めて極端な場合である。グループごとに識別してプロットしてみよう。

------------------------------------------------------

x <- c(17.6, 5.4, 15.5, 8.6, 24.6, 14.5, 27.9, 20.8, 

       13.4, 19.5, 6.8, 0.9, 24.4, 27.1, 17.6, 4.1, 

       1.6, 21.1, 23.6, 19.3, 10.4, 11.1, 19.2, 9.9,

       11.4, 4, 22.8, 21.3, 29.6, 21.2, 2.9, 17.7, 

       24.6, 21.3, 14.7, 14.6, 3.5, 15.7, 28.2, 3.3)

n <- c(24, 35, 23, 25, 30, 27, 25, 35, 23, 29, 25, 29,

       27, 25, 20, 25, 26, 30, 35, 35, 35, 32, 34, 24,

       36, 22, 32, 40, 32, 40, 30, 21, 27, 32, 29, 23, 

       31, 34, 31, 39)

y <- c(0, 0, 0, 0, 1, 0, 0, 0, 0, 2, 0, 0, 1, 1, 0, 1,

       0, 0, 0, 0, 20, 15, 25, 16, 20, 9, 27, 28, 29, 

       26, 20, 16, 27, 30, 25, 17, 23, 31, 31, 30)

g <- 4

n <- 10

group <- rep(LETTERS[1:g], each = n)

d <- data.frame(x = x, n = n, group = group, y = y, rate = y/n)

plotn(rate ~ x, data = d, mode = "m", group = "group")#図36の描画

------------------------------------------------------

全体的には大きくばらつくが、グループ内ではそれほどばらつかない。実際、GLMによる解析で、ピアソンのχ^2と自由度の比をとると、普通のGLMでは過分散である。

------------------------------------------------------

sum(residuals(glm(cbind(y,n-y) ~ x, family = binomial(link = "logit")), type = "pearson")^2)/length(y-2)

## [1] 17.89648

------------------------------------------------------

ゆえに、GLMMをしてみる価値があるだろう。GLMMを行って、固定効果の回帰曲線を描く。

------------------------------------------------------

fit1 <- glmmTMB(cbind(y,n-y) ~ x + (1|group), data = d, family = binomial(link = "logit")) #一般化線形混合モデル

fit2 <- summary(fit1)

fit2

##  Family: binomial  ( logit )

## Formula:          cbind(y, n - y) ~ x + (1 | group)

## Data: d

## 

##      AIC      BIC   logLik deviance df.resid 

##    139.5    144.6    -66.8    133.5       37 

## 

## Random effects:

## 

## Conditional model:

##  Groups Name        Variance Std.Dev.

##  group  (Intercept) 9.66     3.108   

## Number of obs: 40, groups:  group, 4

## 

## Conditional model:

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

## (Intercept) -3.09344    1.58463  -1.952   0.0509 .  

## x            0.08600    0.01326   6.484 8.96e-11 ***

## ---

## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

b0e <- coef(fit2)$cond[1,1]

b1e <- coef(fit2)$cond[2,1]

plotn(rate ~ x, data = d)#図37の描画

xx <- seq(min(x), max(x), length = 200)

yy <- 1/(1+exp(-b0e-b1e * xx))

overdraw(points(xx, yy, type = "l"))

------------------------------------------------------

固定効果の切片は-3.1、xの係数は0.086である。各グループごとの回帰曲線は下記の通り。

------------------------------------------------------

ranef(fit1)

## $group

##   (Intercept)

## A   -2.911579

## B   -2.899301

## C    2.357145

## D    3.723796

plotn(rate ~ x, data = d, mode = "m", group = "group")#図38の描画

for(i in 1:g){

  tmp <- d[d$group == unique(d$group)[i],]

  xx <- seq(min(tmp$x), max(tmp$x), length = 200)

  yy <- 1/(1+exp(-b0e-b1e * xx-ranef(fit1)$cond$group[i,1]))

  overdraw(points(xx, yy, type = "l", col = .default_col[i]))

}

------------------------------------------------------

実際、このデータのように、特定のブロックでは結実率が極端に低く、他では極端に高いみたいな状況は存在しうる。こんなデータであっても、GLMMは無理なく解析を行うことができる。

GLMMは過分散に限らず解析可能

　今回は、明らかな問題を抱えるデータ、過分散なデータを解析する体で、GLMMを紹介した。実際には、過分散なデータに限らず、データ内に階層構造を持っているのであればGLMMでランダム効果を考慮に入れた解析を行うのが望ましい。

　また、今回ではポアソン分布や二項分布を取り扱ったが、もちろん、分散を個別に推定可能な離散分布、例えば負の二項分布やベータ-二項分布を仮定してもよいし、データが連続値なら正規分布、ガンマ分布、そしてベータ分布を扱ってもよい、それはデータの性質をよく確認して、適切な確率分布を考慮すれば問題ないだろう。