過分散な0以上n以下の離散値をとるデータのモデリング

　以前、負の二項分布を混合分布の一形態として紹介した。同じく、二項分布にも混合分布を導入し、過分散な割合データに適合できうることを今回紹介しよう。今回は、ベータ-二項分布と呼ばれるベータ分布と二項分布の混合分布を扱う。

　今回も前回と同様に、GLMでありながらも、分布に過分散な確率分布を指定することで解決を図る方法を紹介する。なので、タイトル詐欺気味ではあるのだが、負の二項分布とベータ-二項分布が分布の混合によって得られる恩恵を、わかりやすく示してくれているので、あえてここで紹介させてもらう。本当のGLMMにともかくチャレンジしたい人は次回に期待してほしい。

ベータ-二項分布：0以上n以下の離散値をとる確率分布

　二項分布でモデリングするには、過分散なデータに使える確率分布が、ベータ-二項分布beta binomial distributionである。この分布は0以上n以下の離散値（0, 1, 2……, n）をとる分布である。ベータ-二項分布は2つのパラメータ、試行回数パラメータnと形状パラメータα,βで特徴づけられる。GLMなどで解析するときは、平均パラメータμ = α/(α+β)、dispersion parameter φ = α+βをパラメータとして推定する。「ベータ-二項分布」は二項分布と名前も似ているし、下記のように確率質量関数Probability mass functionも実際類似している。

ベータ-二項分布はベータ分布と二項分布の混合分布

　ベータ二項分布の重要な性質として、この分布はベータ分布と二項分布の混合分布で表されるという点である。この性質があるために、あえて、一般化線形混合モデルの項でベータ-二項分布を紹介している。先に下記のように混合分布の定義から計算してみよう。ベータ-二項分布は二項分布の成功確率pがベータ分布に従ってばらつく場合の分布である。

　下記のB(α,β)はベータ関数を表し、Γをガンマ関数として、B(α,β) = ∫{p^(α-1)}{(1-p)^(β-1)}dp = Γ(α)Γ(β)/Γ(α+β)である、

　この定義から平均と分散を計算してみよう。今回の場合、有限和なので積分との順序交換が成立することに注目し、二項分布やベータ分布の時に行った計算を使いまわすことで、以下の結果を得られる。

かなり複雑なので見方が難しいが、nが大きくなるほど分散が大きくなり、μが0.5に近いほど分散が大きくなり、φが大きくなるほど分散は小さくなることをつかんでおこう。

　続いて、具体的に乱数を生成してデータの様子を確認しよう。Rのデフォルトではベータ-二項分布を扱えないので、自分で分布を混合する必要がある。

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

library(plotn)

library(viridis)

library(glmmTMB)

phi <- 6

i <- 3

n <- 100

p <- rbeta(n = n, shape1 = phi - i, shape2 = i)

rbinom(n = n, size = 30, prob = p)

##   [1] 17 15 12 14 20 12  9 20 22 12 19 18 10 18 19  3 10 18  8 18 20 16 15  8 10

##  [26] 10 11 20 23 18  5 22  7 12 23 21 12 16 24 18  7 19 12 16 24 14  6  8 19 16

##  [51]  4  5 21  8  9 17 12  8 10 16 20 25 20 15 17 22 19  9  6  4 11  7  8 10 15

##  [76] 13 11 15 13 14  7 16 18 12 14 15 20 22 15 18 19 22 12 10 18 20 29 25 12 24

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

このように乱数を生成してみるとわかるが、すべて0以上の整数であることがわかるだろう。いくつかのパターンで二項分布とともに乱数を生成して、ヒストグラムを描いてみよう。

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

phi <- 6

i <- 3

n <- 1000000

p <- rbeta(n = n, shape1 = phi - i, shape2 = i)

rbb <- rbinom(n = n, size = 30, prob = p)

rb <- rbinom(n = n, size = 30, prob = (phi-i)/phi)

d <- data.frame(value = c(rb, rbb), 

                distribution = c(rep("binomial", n), 

                                 rep("beta-binomial", n)))

histn(value ~ distribution, data = d, freq = F, legend = T, leg.sp = 8)#図1の描画

phi <- 6

i <- 2

n <- 1000000

p <- rbeta(n = n, shape1 = phi - i, shape2 = i)

rbb <- rbinom(n = n, size = 30, prob = p)

rb <- rbinom(n = n, size = 30, prob = (phi-i)/phi)

d <- data.frame(value = c(rb, rbb), 

                distribution = c(rep("binomial", n), 

                                 rep("beta-binomial", n)))

histn(value ~ distribution, data = d, freq = F, legend = T, leg.sp = 8)#図2の描画

phi <- 6

i <- 0.5

n <- 1000000

p <- rbeta(n = n, shape1 = phi - i, shape2 = i)

rbb <- rbinom(n = n, size = 30, prob = p)

rb <- rbinom(n = n, size = 30, prob = (phi-i)/phi)

d <- data.frame(value = c(rb, rbb), 

                distribution = c(rep("binomial", n), 

                                 rep("beta-binomial", n)))

histn(value ~ distribution, data = d, freq = F, legend = T, leg.sp = 8)#図3の描画

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

平均的に同じ確率パラメータを持つ二項分布よりも分布がばらつく様子に注目しよう。

　R内にガンマ関数は実装されているので、下記のように、台、サイズパラメータ、平均パラメータ、dispersion parameterを引数とする確率分布を生成する関数を定義できる。これを使って、確率分布を眺めることも可能である。

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

dbeta_binom <- function(x, size, mu, phi){

  f <- function(x, size, mu, phi){

      choose(size, x)*(gamma(mu*phi + x)/gamma(mu*phi))*(gamma((1-mu)*phi + size - x)/gamma(((1-mu)*phi)))*(gamma(phi)/gamma(phi+size))

  }

  mapply(f, x = x, size = size, mu = mu, phi = phi)

}

size <- 30

dd <- data.frame(X = 0:size)

phi <- 10

i <- 9

mu <- (phi-i)/phi

dd$prob01 <- dbeta_binom(dd$X, size = size, mu = mu, phi = phi)

phi <- 10

i <- 5

mu <- (phi-i)/phi

dd$prob05 <- dbeta_binom(dd$X, size = size, mu = mu, phi = phi)

phi <- 10

i <- 3

mu <- (phi-i)/phi

dd$prob07 <- dbeta_binom(dd$X, size = size, mu = mu, phi = phi)

plotn(dd[,1], dd[,2:4], mode = "m", legend = T, leg.title = "prob",

      leg.lab = c(0.1, 0.5, 0.7), leg.sp = 5, ylab = "確率質量分布")#図4の描画

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

μが1に近づくにつれて右よりの分布になることが確認できる。同じ平均確率を持つ二項分布も描いてみよう。

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

dd <- data.frame(X = 0:size)

dd$prob01 <- dbinom(dd$X, size = size, prob = 0.1)

dd$prob05 <- dbinom(dd$X, size = size, prob = 0.5)

dd$prob07 <- dbinom(dd$X, size = size, prob = 0.7)

plotn(dd[,1], dd[,2:4], mode = "m", legend = T, leg.title = "prob",

      leg.lab = c(0.1, 0.5, 0.7), leg.sp = 5, ylab = "確率質量分布")#図5の描画

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

確かに、二項分布よりもベータ-二項分布はばらついている。

　今度は、φのほうの値を変えてみよう。

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

size <- 30

dd <- data.frame(X = 0:size)

phi <- 1

i <- 1/2

mu <- (phi-i)/phi

dd$disp1_2 <- dbeta_binom(dd$X, size = size, mu = mu, phi = phi)

phi <- 2

i <- 1

mu <- (phi-i)/phi

dd$disp2 <- dbeta_binom(dd$X, size = size, mu = mu, phi = phi)

phi <- 8

i <- 4

mu <- (phi-i)/phi

dd$disp8 <- dbeta_binom(dd$X, size = size, mu = mu, phi = phi)

phi <- 32

i <- 16

mu <- (phi-i)/phi

dd$disp32 <- dbeta_binom(dd$X, size = size, mu = mu, phi = phi)

plotn(dd[,1], dd[,2:5], mode = "m", legend = T, leg.title = "phi",

      leg.lab = c(1,2,8,32), leg.sp = 5, ylab = "確率質量分布")#図6の描画

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

ベータ-二項分布を使ったGLM（ベータ-二項回帰）

　では、ベータ-二項分布を使ったGLMについて考察してゆこう。ベータ-二項分布を使う場合、Rのプログラム上では、リンク関数はデフォルトでロジットリンクを使うことになっている。

GLMではこのように確率分布とリンク関数を決めたのち、尤度を最大化するようなパラメータβを求める。なので、次に尤度関数を求めてみよう。データの生成過程を図示して様子をつかむことにする。 

ベータ分布が絡んでいるので、説明変数xがなくても最尤推定値が簡単にならないだろうと予測付くので、ここでは計算を行わない。

一般化線形モデルの検出力

　ではRを使って、GLMの検出力に関して検討しよう。以下のように、線形予測子の切片の係数は0ではないが、説明変数xの係数は0である場合をまず考える。リンク関数はロジットリンクを想定する。ベータ分布の形状パラメータα、βを指定する必要があるが、φを固定して、μi = αi/φ = 1/(1+exp(-β・xi))なので、αi = φ/(1+exp(-β・xi))、βi = φ(1-1/(1+exp(-β・xi)))と指定できる。データを生成した後で、描画してみよう。

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

b0 <- -1

b1 <- 0

phi <- 5

m <- 30

x <- runif(m, 0, 30)

p <- rbeta(n = m, shape1 = phi/(1+exp(-b0-b1*x)), shape2 = phi*(1-1/(1+exp(-b0-b1*x))))

n <- sample(20:30, size = m, replace = T)

y <- rbinom(n = m, size = n, prob = p)

plotn(y/n ~ x)#図7の描画

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

データの生成過程通り、xに対して特に関係のなさそうなデータが生成された。では、このデータをGLMで解析する。今回は、デフォルトのglm関数ではベータ-二項分布を扱えないので、glmmTMBパッケージ内の、glmmTMB関数とbetabinomial関数を使う。ここではベータ-二項分布とロジットリンクを指定するので、betabinomial(link="logit")と指定する。併せて、回帰曲線も描いてみよう。

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

fit1 <- glmmTMB(cbind(y, n-y) ~ x, family = betabinomial(link = "logit")) #一般化線形モデル

fit2 <- summary(fit1)

fit2

##  Family: betabinomial  ( logit )

## Formula:          cbind(y, n - y) ~ x

## 

##      AIC      BIC   logLik deviance df.resid 

##    183.4    187.6    -88.7    177.4       27 

## 

## 

## Dispersion parameter for betabinomial family (): 5.29 

## 

## Conditional model:

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

## (Intercept) -0.73086    0.32369  -2.258   0.0239 *

## x           -0.00164    0.01761  -0.093   0.9258  

## ---

## 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/n ~ x)#図8の描画

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

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

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

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

今回の場合、線形予測子の切片の推定値は-0.73、xの係数の推定値は-0.0016であった。

　普通の二項分布を仮定したGLMを使うことで実際にこのデータが過分散であることを示そう。ピアソンのχ^2と自由度の比を計算する。

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

fit3 <- glm(cbind(y, n-y) ~ x, family = binomial(link = "logit"))

fit4 <- summary(fit3)

fit4

## 

## Call:

## glm(formula = cbind(y, n - y) ~ x, family = binomial(link = "logit"))

## 

## Deviance Residuals: 

##     Min       1Q   Median       3Q      Max  

## -4.4861  -1.3762  -0.3787   1.2096   5.2668  

## 

## Coefficients:

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

## (Intercept) -0.683985   0.150471  -4.546 5.48e-06 ***

## x           -0.003639   0.008414  -0.433    0.665    

## ---

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

## 

## (Dispersion parameter for binomial family taken to be 1)

## 

##     Null deviance: 158.98  on 29  degrees of freedom

## Residual deviance: 158.79  on 28  degrees of freedom

## AIC: 257.94

## 

## Number of Fisher Scoring iterations: 4

sum(residuals(fit3, type = "pearson")^2)/(length(y)-2)

## [1] 5.159779

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

確かに、比は5.1となり二項分布にデータが従っているのであれば1になることが期待されるから、このデータは二項分布に対して過分散である。

　では、上記の設定でデータを1000回生成して、検出力を確認してみる。glmmTMBは計算がやや遅いので、1000回程度でとどめる。併せて、ベータ-二項分布から生成された過分散なデータで二項回帰も行ってみよう。

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

b0 <- -1

b1 <- 0

phi <- 5

m <- 30

p1bb <- NULL

p2bb <- NULL

b0ebb <- NULL

b1ebb <- NULL

p1b <- NULL

p2b <- NULL

b0eb <- NULL

b1eb <- NULL

for (i in 1:1000){

  x <- runif(m, 0, 30)

  p <- rbeta(n = m, shape1 = phi/(1+exp(-b0-b1*x)), shape2 = phi*(1-1/(1+exp(-b0-b1*x))))

  n <- sample(20:30, size = m, replace = T)

  y <- rbinom(n = m, size = n, prob = p)

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

  fit1 <- glmmTMB(cbind(y, n-y) ~ x, data = d, family = betabinomial(link = "logit"))

  fit2 <- summary(fit1)

  fit3 <- glm(cbind(y, n-y) ~ x, data = d, family = binomial(link = "logit"))

  fit4 <- summary(fit3)

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

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

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

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

  p1b <- c(p1b, coef(fit4)[1,4])

  p2b <- c(p2b, coef(fit4)[2,4])

  b0eb <- c(b0eb, coef(fit4)[1,1])

  b1eb <- c(b1eb, coef(fit4)[2,1])

}

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

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

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

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

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

明らかに二項分布を仮定した時のxの係数のP値のヒストグラムは0付近に偏っており、危険率が高まってることがわかる。危険率を計算してみる。

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

sum(p1bb < 0.05)#ベータ-二項分布を仮定した時の切片の検出力

## [1] 809

sum(p2bb < 0.05)#ベータ-二項分布を仮定した時のx係数の危険率

## [1] 84

sum(p3 < 0.05)#二項分布を仮定した時の切片の検出力

## [1] 953

sum(p4 < 0.05)#二項分布を仮定した時のx係数の危険率

## [1] 407

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

すると、ベータ-二項分布を仮定した時は、xの係数の危険率は8%程度なのに対して、二項分布を仮定すると41%程度まで危険率が高まった。ベータ-二項分布でもやや危険率が高いが、こちらを仮定することがより適切であろう。

　次はxの係数が存在している場合を考えてみる。

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

b0 <- -3

b1 <- 0.2

phi <- 5

m <- 30

x <- runif(m, 0, 30)

p <- rbeta(n = m, shape1 = phi/(1+exp(-b0-b1*x)), shape2 = phi*(1-1/(1+exp(-b0-b1*x))))

n <- sample(20:30, size = m, replace = T)

y <- rbinom(n = m, size = n, prob = p)

plotn(y/n ~ x)#図13の描画

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

今回はxが増加するとともに、yも増加している傾向が見て取れる。ではGLMで解析しよう。

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

fit1 <- glmmTMB(cbind(y, n-y) ~ x, family = betabinomial(link = "logit")) #一般化線形モデル

fit2 <- summary(fit1)

fit2

##  Family: betabinomial  ( logit )

## Formula:          cbind(y, n - y) ~ x

## 

##      AIC      BIC   logLik deviance df.resid 

##    160.0    164.2    -77.0    154.0       27 

## 

## 

## Dispersion parameter for betabinomial family (): 7.45 

## 

## Conditional model:

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

## (Intercept) -2.66138    0.45629  -5.833 5.45e-09 ***

## x            0.19097    0.02698   7.078 1.47e-12 ***

## ---

## 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/n ~ x)#図14の描画

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

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

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

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

今回の場合、線形予測子の切片の推定値は-2.7、xの係数の推定値は0.19であった。

　普通の二項分布を仮定したGLMを使うことで実際にこのデータが過分散であることを示そう。ピアソンのχ^2と自由度の比を計算する。

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

fit3 <- glm(cbind(y, n-y) ~ x, family = binomial(link = "logit"))

fit4 <- summary(fit3)

fit4

## 

## Call:

## glm(formula = cbind(y, n - y) ~ x, family = binomial(link = "logit"))

## 

## Deviance Residuals: 

##     Min       1Q   Median       3Q      Max  

## -3.8686  -1.4726   0.5099   1.2722   4.2102  

## 

## Coefficients:

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

## (Intercept) -2.56153    0.22794  -11.24   <2e-16 ***

## x            0.18703    0.01395   13.41   <2e-16 ***

## ---

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

## 

## (Dispersion parameter for binomial family taken to be 1)

## 

##     Null deviance: 388.01  on 29  degrees of freedom

## Residual deviance: 124.17  on 28  degrees of freedom

## AIC: 210.95

## 

## Number of Fisher Scoring iterations: 5

sum(residuals(fit3, type = "pearson")^2)/(length(y)-2)

## [1] 4.408651

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

確かに、比は4.4となり二項分布にデータが従っているのであれば1になることが期待されるから、このデータは二項分布に対して過分散である。

　では、上記の設定でデータを1000回生成して、検出力を確認してみる。

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

b0 <- -3

b1 <- 0.2

phi <- 5

m <- 30

p1bb <- NULL

p2bb <- NULL

b0ebb <- NULL

b1ebb <- NULL

for (i in 1:1000){

  x <- runif(m, 0, 30)

  p <- rbeta(n = m, shape1 = phi/(1+exp(-b0-b1*x)), shape2 = phi*(1-1/(1+exp(-b0-b1*x))))

  n <- sample(20:30, size = m, replace = T)

  y <- rbinom(n = m, size = n, prob = p)

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

  fit1 <- glmmTMB(cbind(y, n-y) ~ x, data = d, family = betabinomial(link = "logit"))

  fit2 <- summary(fit1)

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

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

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

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

}

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

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

sum(p1bb < 0.05)

## [1] 1000

sum(p2bb < 0.05)

## [1] 1000

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

すると切片の係数の検出力は100%、xの係数の検出力は100%程度である。xの係数の設定を変えたことで、xの係数が有意に検出されるようになった。

　最後、シミュレーションだけだが、切片だけが0の場合について、検討してみよう。

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

b0 <- 0

b1 <- 0.2

phi <- 5

m <- 30

p1bb <- NULL

p2bb <- NULL

b0ebb <- NULL

b1ebb <- NULL

p1b <- NULL

p2b <- NULL

b0eb <- NULL

b1eb <- NULL

for (i in 1:1000){

  x <- runif(m, 0, 30)

  p <- rbeta(n = m, shape1 = phi/(1+exp(-b0-b1*x)), shape2 = phi*(1-1/(1+exp(-b0-b1*x))))

  n <- sample(20:30, size = m, replace = T)

  y <- rbinom(n = m, size = n, prob = p)

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

  fit1 <- glmmTMB(cbind(y, n-y) ~ x, data = d, family = betabinomial(link = "logit"))

  fit2 <- summary(fit1)

  fit3 <- glm(cbind(y, n-y) ~ x, data = d, family = binomial(link = "logit"))

  fit4 <- summary(fit3)

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

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

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

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

  p1b <- c(p1b, coef(fit4)[1,4])

  p2b <- c(p2b, coef(fit4)[2,4])

  b0eb <- c(b0eb, coef(fit4)[1,1])

  b1eb <- c(b1eb, coef(fit4)[2,1])

}

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

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

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

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

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

やはり、二項分布を仮定すると切片の危険率が高まっていることがわかる。

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

sum(p1bb < 0.05)#ベータ-二項分布を仮定した時の切片の危険率

## [1] 73

sum(p2bb < 0.05)#ベータ-二項分布を仮定した時のx係数の検出力

## [1] 999

sum(p3 < 0.05)#二項分布を仮定した時の切片の危険率

## [1] 382

sum(p4 < 0.05)#二項分布を仮定した時のx係数の検出力

## [1] 1000

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

一般化線形モデルを使った解析

　では、今度はデータを与えられた時の解析の方針などを考えてゆこう。例えば、次のようなデータが与えられたとする。データの図示はまず最初にやろう。

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

n <- c(27, 29, 22, 23, 29, 21, 21, 20, 30, 30, 27, 26,

       21, 29, 29, 30, 29, 28, 20, 30, 21, 24, 27, 22, 

       22, 28, 20, 26, 22, 25, 22, 30, 26, 28, 26, 21, 

       20, 26, 28, 26, 30, 22, 30, 30, 28, 21, 20, 30, 

       23, 25, 23, 26, 22, 29, 20, 23, 22, 24, 24, 28,

       21, 27, 25, 24, 30, 28, 27, 25, 29, 20, 28, 30,

       26, 20, 30, 24, 24, 24, 22, 29, 27, 22, 20, 23, 

       27, 28, 27, 23, 25, 25, 28, 22, 27, 28, 30, 24, 

       22, 20, 23, 30)

y <- c(21, 20, 10, 15, 11, 12, 16, 15, 19, 12, 15, 23, 

       17, 22, 13, 17, 18, 20, 12, 16, 16, 20, 19, 20, 

       14, 26, 9, 15, 20, 9, 11, 14, 15, 24, 14, 14, 9,

       12, 19, 14, 16, 13, 19, 20, 21, 15, 18, 25, 9, 

       8, 12, 14, 11, 8, 16, 19, 14, 13, 5, 25, 6, 9, 

       13, 17, 20, 6, 22, 21, 11, 11, 20, 21, 21, 5, 18, 

       19, 21, 18, 10, 11, 18, 15, 15, 11, 11, 23, 15, 

       23, 16, 24, 26, 15, 14, 15, 12, 14, 17, 11, 19, 15)

histn(y/n)#図21の描画

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

データは0以上の離散値データであり、カウントデータの割合を解析しようとしている。二項分布よりも分散は大きそうである。こういう時はベータ-二項分布を解析に使うとよい。一旦は練習なので、説明変数がないデータの解析をしてみることにする。GLMで解析してみると下記のようになる。

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

fit1 <- glmmTMB(cbind(y, n-y) ~ 1, family = betabinomial(link = "logit"))

fit2 <- summary(fit1)

fit2

##  Family: betabinomial  ( logit )

## Formula:          cbind(y, n - y) ~ 1

## 

##      AIC      BIC   logLik deviance df.resid 

##    586.7    591.9   -291.3    582.7       98 

## 

## 

## Dispersion parameter for betabinomial family ():  8.4 

## 

## Conditional model:

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

## (Intercept)   0.5055     0.0777   6.506 7.72e-11 ***

## ---

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

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

すると切片の係数は有意に0と異なっていることがわかる。さてここで尤度のことを思い出そう。今回推定された係数は尤度を最大にするように選ばれているはずである。対数尤度を計算するために下記のように関数を定義する。上記で作成した、ベータ-二項分布の確率質量を返す関数はここでも使うことができる。自作の確率分布で尤度を計算できるのは、少し驚きかもしれない。

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

f <- function(y, n, phi, b0t){

  sapply(b0t, function(tmp) sum(log(dbeta_binom(x = y, size = n, 

                                                mu = 1/(1+exp(-tmp)),

                                                phi = phi)

                                    )

                                )

         )

}#あるb0に対応する対数尤度を計算

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

この時、計算のためにdispersion parameterが必要で、それはfit2$sigmaに格納されている。この値は、glm関数でガンマ回帰をした時と異なり、最尤推定されている。今回得られた最尤推定値b0e = 0.51から±1の範囲の連続値を生成し、それらに対応する対数尤度を計算して、図示してみる。

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

phi_e <- fit2$sigma

db0 <- 1

b0min <- b0e - db0

b0max <- b0e + db0

logL <- f(y, n, phi_e, seq(b0min, b0max, length = 2000))

plotn(seq(b0min, b0max, length = 2000), logL, xlab = "b0", ylab = "logL", type = "l")#図22の描画

overdraw(abline(v = b0e),

         axis(side = 1, at = b0e, labels = "^b0"))

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

この図を見るとわかるように、確かに推定値b0eで最大の値を示していることがわかるだろう。さらに狭い範囲で図示すると次のようになる。

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

db0 <- 0.01

b0min <- b0e - db0

b0max <- b0e + db0

logL <- f(y, n, phi_e, seq(b0min, b0max, length = 2000)) #対数尤度の計算

plotn(seq(b0min, b0max, length = 2000), logL, xlab = "b0", ylab = "logL", type = "l")

overdraw(abline(v = b0e),

         axis(side = 1, at = b0e, labels = "^b0")) #図23の図示

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

glmmTMBの出力は、glmと異なり、対数尤度と逸脱度が表示される。下記のように上記の関数を使って計算できる。併せてAICも計算しよう。ガンマ回帰、ベータ回帰、負の二項分布の時と同様に、dispersion parameterが線形予測子以外に推定されているので、逸脱度への加算量は2×2である。

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

f(y, n, phi_e, b0e)

## [1] -291.325

-2*f(y, n, phi_e, b0e)

## [1] 582.6501

-2*f(y, n, phi_e, b0e) + 2*2

## [1] 586.6501

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

　続いて、説明変数xを含んだ時の解析を紹介する。例えば、次のようなデータである。図示も併せて行おう。

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

x <- c(5.6, 19.9, 3.1, 2.6, 0.3, 5, 13.7, 27.3, 19.8, 9.5, 8, 

       9.4, 22.8, 5.9, 29, 20.5, 22.3, 14.1, 3.8, 1.4, 2.9, 

       7.7, 23.5, 17.1, 7, 0.3, 28.6, 26.3, 8.9, 3.4, 3.2, 9.2, 

       5, 27.7, 5.8, 25.4, 1, 17.4, 28.4, 10.2, 8.9, 18.7, 1.8,

       7, 3.5, 24.5, 0.1, 18.1, 16.1, 4.5)

n <- c(26, 27, 30, 23, 20, 27, 21, 30, 24, 21, 30, 28, 23, 28,

       27, 30, 30, 22, 26, 30, 21, 29, 29, 22, 26, 24, 28, 20,

       20, 26, 28, 28, 24, 27, 28, 26, 26, 22, 26, 27, 29, 20,

       21, 25, 22, 30, 25, 20, 29, 24)

y <- c(15, 5, 7, 16, 17, 10, 2, 0, 6, 7, 15, 6, 3, 8, 0, 5, 6,

       5, 12, 22, 20, 13, 1, 2, 9, 17, 0, 2, 15, 14, 18, 6, 18,

       0, 14, 3, 22, 5, 6, 10, 14, 3, 18, 16, 18, 0, 14, 2, 15, 19)

plotn(y/n ~ x)#図24の図示

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

では解析をしてみよう。解析ののち、回帰曲線を描く。

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

fit1 <- glmmTMB(cbind(y, n-y) ~ x, family = betabinomial(link = "logit")) #一般化線形モデル

fit2 <- summary(fit1)

fit2

##  Family: betabinomial  ( logit )

## Formula:          cbind(y, n - y) ~ x

## 

##      AIC      BIC   logLik deviance df.resid 

##    258.6    264.4   -126.3    252.6       47 

## 

## 

## Dispersion parameter for betabinomial family (): 11.6 

## 

## Conditional model:

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

## (Intercept)  0.99295    0.17604   5.641  1.7e-08 ***

## x           -0.13911    0.01528  -9.104  < 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]

phi_e <- fit2$sigma

plotn(y/n ~ x)#図25の図示

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

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

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

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

では、今回も対数尤度の定義から、2変数の場合の対数尤度を求めて図示してみよう。推定値b0e = 0.99、b1e = -0.14周りの対数尤度は下記のように計算できる。

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

f <- function(x, y, n, phi, b0t, b1t){

  mapply(function(tmp1, tmp2) {

    sum(log(dbeta_binom(x = y, size = n, mu = 1/(1+exp(-tmp1-tmp2*x)), phi = phi)

            )

    )

    }, b0t, b1t)

}#あるb0, b1に対応する対数尤度を計算

db0 <- 0.01

db1 <- 0.001

b0min <- b0e - db0

b0max <- b0e + db0

b1min <- b1e - db1

b1max <- b1e + db1

logL <- outer(seq(b0min, b0max, length = 200), seq(b1min, b1max, length = 200), f, x = x, y = y, n = n, phi = phi_e)

lev <- unique(c(seq(min(logL), max(logL)- 0.01, length = 10), 

                seq(max(logL)-0.01, max(logL)- 0.0001, length = 10), 

                seq(max(logL) - 0.0005, max(logL)-0.00001, length = 5)))

image(seq(b0min, b0max, length = 200), seq(b1min, b1max, length = 200), logL, col = magma(400),

      xlab = "b0", ylab = "b1", las = 1)#図26の図示

contour(seq(b0min, b0max, length = 200), seq(b1min, b1max, length = 200), logL,

        levels = lev, add = T)

abline(v = b0e)

abline(h = b1e)

axis(side = 1, at = b0e, labels = "^b0")

axis(side = 2, at = b1e, labels = "^b1")

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

すると、確かに最尤推定値のb0eとb1eがその周りで最も高い尤度となっていることがわかる。以下のような図示の仕方もある。

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

persp(seq(b0min, b0max, length = 200), seq(b1min, b1max, length = 200), logL, 

      xlab = "b0", ylab = "b1", zlab = "logL",

      theta=120, phi=20, expand=0.5, ticktype="detailed")#図27の図示

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

上記の関数を使って対数尤度、逸脱度を計算できる。併せてAICも計算しよう。ガンマ回帰の時と同様に、dispersion parameterが線形予測子以外に推定されているので、逸脱度への加算量は2×3である。

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

f(x, y, n, phi_e, b0e, b1e)

## [1] -126.3227

-2*f(x, y, n, phi_e, b0e, b1e)

## [1] 252.6453

-2*f(x, y, n, phi_e, b0e, b1e) + 2*3

## [1] 258.6453

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