ハードルモデルともう一つのモデル

　前回、ゼロをとらない確率分布に対して、ハードルモデルを使って解析した。その時扱ったガンマ分布は、ゼロをとらないので、ゼロ切断分布はガンマ分布そのものであり、解釈も難しくなかった。今回は、ポアソン分布などのゼロも出現する確率分布を扱う。離散分布であっても同じようにハードルモデルを適用できるが、ゼロの出現確率に関して調整が必要である。さらに本項ではハードルモデルとともに、ゼロ過剰な離散値を扱えるモデルとして、ゼロ過剰モデルを扱う。よく似たモデルだが、ゼロの扱いが異なっている。

0が過剰な離散値データ：ハードルモデル的考え方

　ゼロが過剰なデータとは、その名の通り、期待されるよりもゼロの出現確率が高いデータである。このようなデータに関してゼロの出現メカニズムを無視して解析を行うと解析結果がゆがめられてしまう。そこで、ゼロの出現メカニズムについて考察しよう。

　前回と同様に、ハードルモデル的なゼロの出現が考えられる。つまり、ある見えないハードルがあって、そのハードルを越えられなければゼロ、越えたら正の確率分布に従う、と考える。例えば、ある（自動自家受粉しない）植物の種子生産数を調査したので、それをモデリングすることを考える。

個数という離散値のモデリングなので、ポアソン分布などを仮定するのが良いだろう。あるいはばらつくことを考えて慎重な人は負の二項分布を考慮するかもしれない。それでもたぶんこの解析は失敗する可能性がある。この植物は自動自家受粉しないので、結実するためには送粉者が来る必要がある（送粉者が来れば絶対結実するとする）。正確にデータを記述するなら送粉者の有無は観察しておくべきだった。

しかし、実際、それが難しいこともある。それに、得られているデータは種子数だけなので、ここからデータを読み解く必要がある。そこで見方を変えれば、結実していれば「送粉者が来るというハードル」を乗り越えており、結実していなければハードルを越えられていない、と考えることができる。こう捉えれば、ハードルモデルで解析することができる。

0が過剰な離散値データ：ゼロ過剰モデル的考え方

　対してゼロ過剰モデルでは、次のような例が適切だろう。あるチョウの訪花回数をモデリングすることを考える。

これも回数なのでポアソン分布や負の二項分布に従うと考えてモデリングすることを考えるだろう。ところで、チョウの状態に注目すると、そもそも訪花するためには花の存在を認識できている必要があるだろう。さらに、花の存在を認識できていたとして、訪花しないこともあり得る。例えば、満腹状態であるとか、そういった理由があるかもしれない。

つまり、ゼロに2つの異なる意味が存在する。この場合は、ゼロ過剰モデルで解析するのが適切である。

ゼロ過剰モデル：2つのゼロが混在するモデル

　ゼロ過剰モデルzero-inflated modelは0が過剰に生じるようなデータで採用されるモデリング手法の一つである。ハードルモデルと同様に、その確率分布は、ベルヌーイ分布とゼロ以上の離散確率分布の混合分布として表現される。今回は、ポアソン分布との混合を考えよう。ゼロ過剰モデルはある確率1-pでゼロが生成されるだけでなく、残りpの確率でポアソン分布に従い、そのポアソン分布からもゼロが生成される。

このように、元の離散分布（ベース分布）に対して、1-pだけ余分にゼロが生成されるわけである。後述するように、p = 1のとき、ベース分布に一致する。

　離散分布でのハードルモデルの模式図も掲載しておく。

ハードルモデルの場合、仮にp = 1となってもベース分布に一致せず、ベース分布からゼロが生じる確率を除去したゼロ切断分布zero-truncated distributionとなる。実はこの性質に従うと、ハードルモデルはゼロ過剰だけでなく、ゼロ過少もモデリングできる（この性質が後で重要になる）。

　ゼロ過剰モデルとハードルモデルの両方、確率密度関数を紹介する。

ハードルモデルは以前紹介したとおりである。ゼロ過剰モデルを紹介すると、x = 0となる確率のうちp g(0|θ)と、x > 0となる確率を合わせればちょうどベース分布をp倍したものである。ベース分布の確率質量の和は1となるから、ベース分布から生成される確率はpとなる。残りの1-pはゼロが生成されるということだ。このベース分布以外から生成されるゼロを構造的なゼロstructural zeroと呼ぶ。

　例えば、ベース分布をポアソン分布とするようなハードルモデルやゼロ過剰モデルは以下のようになる。

　ゼロ過剰モデルはハードルモデルと同様に、2つの線形予測子とリンク関数を持つ。ベルヌーイ分布を支配するpに関してはロジットリンクを指定する。ベース分布に関しては、そのベース分布にあったものを選ぶが、ポアソン分布なら対数リンクを指定することが多いだろう。ゆえにゼロ過剰モデルもtwo-parts modelである。

具体的な解析

　例えば、ゼロ過剰モデルを想定した乱数の生成は下記のように行える。あらかじめ、線形予測子とリンク関数を想定した入力にしている。

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

library(plotn)

library(actuar)

library(glmmTMB)

r_zipois <- function(n, lambda, prob){

  tmp <- rbinom(n = n, size = 1, prob = prob)

  v <- mapply(function(x, lambda){

    if(x == 1) {

      rpois(n = 1, lambda)

      } else {0}

    },

    x = tmp, lambda = lambda

    )

  unlist(v)

}

b0 <- 0.1

b0_h <- -0.5

n <- 100

r_zipois(n = n,  lambda = exp(b0), prob = 1/(1+exp(b0_h)))

##   [1] 2 0 0 1 0 1 1 0 1 2 2 0 1 0 0 0 3 3 0 0 0 0 0 0 1 0 5 0 2 0 0 1 0 0 2 0 4

##  [38] 0 2 0 0 2 0 0 1 0 0 0 0 0 0 0 0 0 0 0 4 0 0 0 1 0 0 0 0 0 0 0 0 0 4 4 0 0

##  [75] 1 0 1 2 0 2 1 1 1 0 3 0 0 0 0 1 0 0 0 0 0 5 1 0 3 0

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

一方で、ハードルモデルを想定した乱数の生成はひと手間必要である。というのも、ゼロ切断分布による乱数生成が必要になるからである。そこでactuarパッケージの、rztpois関数を使うことにする（ztはzero truncatedの意味）。この関数はrpois関数と同様にlambda引数を持ち、平均がlambdaのポアソン分布に従う乱数を生成するが、ゼロは生成しない。ゆえにrztpoisから生成された乱数の平均はlambdaとは一致しないので注意である。

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

r_hpois <- function(n, lambda, prob){

  tmp <- rbinom(n = n, size = 1, prob = prob)

  v <- mapply(function(x, lambda){

    if(x == 1) {

      rztpois(n = 1, lambda)

      } else {0}

    },

    x = tmp, lambda = lambda

    )

  unlist(v)

}

b0 <- 0.1

b0_h <- -0.5

n <- 100

r_hpois(n = n,  lambda = exp(b0), prob = 1/(1+exp(b0_h)))

##   [1] 2 2 1 2 0 3 3 1 2 0 1 2 1 1 0 1 1 3 3 1 0 2 0 2 1 1 0 2 0 0 0 1 0 1 0 1 0

##  [38] 0 1 0 0 1 2 1 1 0 1 0 1 1 1 0 2 0 2 2 1 0 3 1 1 1 0 1 2 0 0 0 1 1 2 1 2 0

##  [75] 1 1 0 3 0 0 3 2 3 0 0 2 0 0 2 1 0 0 1 3 1 0 0 2 2 0

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

もっと乱数を生成してヒストグラムを描いてみよう。

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

n <- 10000

x <- r_hpois(n = n,  lambda = exp(b0), prob = 1/(1+exp(b0_h)))

histn(x, freq = F)#図1の描画

x <- r_zipois(n = n,  lambda = exp(b0), prob = 1/(1+exp(b0_h)))

histn(x, freq = F)#図2の描画

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

設定にもよるが、ゼロ過剰モデルのほうがより多くのゼロが生成される。

・ハードルモデルによる解析

　では、下記のようなデータを解析してみよう。とりあえず、ポアソン分布に従うものとして考える。

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

x <- c(2, 26, 29.2, 10.6, 12.1, 8.9, 4.5, 29, 22.6, 15.8, 25.4,

       25.1, 2.4, 4.3, 15, 20.2, 18.6, 6.6, 22.9, 22, 4.9, 24.6, 

       6.5, 23.3, 14, 11.7, 24, 29.4, 14.7, 9.9, 12.7, 13.5, 21.5, 

       18.3, 4.9, 12.4, 22.9, 4.1, 8.1, 11.4, 29.4, 21.6, 3.5,

       15.4, 12.1, 13.2, 14.5, 30, 0.6, 3.9, 4.4, 17.6, 17.9, 5.5,

       15.1, 12.2, 21.9, 5.3, 6.7, 25, 5.8, 15.2, 8.9, 14, 14.3, 

       26.1, 8.5, 16.5, 2.1, 13.1, 1.9, 4.2, 20.9, 16.7, 19.9, 10.3, 

       19.7, 15.1, 9.1, 12.7, 15.1, 20.9, 15.8, 15.7, 2.9, 19.4, 

       5.2, 6.2, 17.7, 21, 12.5, 10.9, 19.6, 23.4, 27, 25.9, 

       22.5, 22.6, 2.5, 28.6)

y <- c(6, 0, 0, 2, 12, 8, 6, 0, 0, 0, 0, 3, 11, 9, 0, 3, 0, 5,

       0, 5, 5, 0, 8, 0, 3, 6, 0, 0, 0, 12, 0, 2, 3, 0, 12, 11,

       0, 11, 8, 0, 0, 4, 5, 6, 2, 3, 9, 0, 13, 15, 10, 1, 1, 

       12, 7, 6, 0, 4, 4, 0, 5, 4, 4, 2, 0, 0, 5, 2, 18, 0, 10, 

       16, 0, 0, 0, 2, 0, 3, 11, 0, 3, 0, 0, 3, 17, 0, 8, 9, 0, 

       0, 2, 4, 0, 0, 0, 0, 1, 0, 14, 0)

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

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

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

データはxが増加するほどyは低い値を示す。

　y > 0の値を1に変換して、y2に格納して、0/1データとしてプロットし、0の出現について考えてみる。あとで使うので、y = 0をNAに変換してy3に格納もしておく。

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

d$y2 <- d$y

d$y2[d$y2 > 0] <- 1

d$y3 <- d$y

d$y3[d$y3 == 0] <- NA

plotn(y2 ~ x, data = d)#図4の描画

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

すると、0の出現はxが増えるほど増加する。

　データの生成に関して、何のメカニズムも想定していないので、このデータだけから通常のポアソン回帰にするか、ハードルモデルにするか、ゼロ過剰モデルにするかを決めるのは困難である。ただし、現状のデータは単純にポアソン回帰すると過分散であることがわかる。

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

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

## [1] 2.393472

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

そこで、天下り的であるが、今回はハードルモデルを想定して解析してみよう。前回同様にglmmTMBパッケージのglmmTMB関数を使う。ポアソン回帰するときのglm関数などと同じように表記してよい。そしてziformula引数には、ベルヌーイ分布を支配するパラメータpに関する線形予測子を入れる。今回は、ベルヌーイ分布のpはxと関係がありそうなので、~ xとして推定させよう。確率分布はポアソン分布とするが、ゼロ切断分布が必要になる。glmmTMBパッケージのtruncated_poisson分布がそれに相当する。ポアソン回帰の回帰曲線も描いてみよう。

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

fit1 <- glmmTMB(y ~ x, data = d, family = truncated_poisson(link = "log"), 

                ziformula = ~ x) #一般化線形モデル

fit2 <- summary(fit1)

fit2

##  Family: truncated_poisson  ( log )

## Formula:          y ~ x

## Zero inflation:     ~x

## Data: d

## 

##      AIC      BIC   logLik deviance df.resid 

##    380.7    391.1   -186.3    372.7       96 

## 

## 

## Conditional model:

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

## (Intercept)  2.62316    0.09035  29.033  < 2e-16 ***

## x           -0.08219    0.01006  -8.167 3.16e-16 ***

## ---

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

## 

## Zero-inflation model:

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

## (Intercept) -5.15921    1.00207  -5.149 2.62e-07 ***

## x            0.29943    0.05762   5.197 2.03e-07 ***

## ---

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

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

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

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

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

yy <- exp(b0e + b1e * xx)

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

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

そして、ベルヌーイ回帰の結果は以下のように描ける。

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

b0e_h <- coef(fit2)$zi[1,1]

b1e_h <- coef(fit2)$zi[2,1]

plotn(y2 ~ x, data = d)#図6の描画

xx <- seq(0, 30, length = 200)

yy <- 1/(1+exp(b0e_h + b1e_h * xx))

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

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

また、ポアソン回帰とベルヌーイ回帰の両方の結果を考慮した回帰曲線は、ガンマ-ハードルモデルの時と同様に以下のように描ける。

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

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

xx <- seq(0, 30, length = 200)

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

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

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

30付近でよりゼロに近づくような期待値となる回帰曲線が描けた。

　そしてこのモデルは過分散が解消されている。

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

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

## [1] 1.298825

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

　実際、今回の解析は下記のように設定してデータを生成したので、悪くない推定になっていることがわかるだろう。

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

b0 <- 2.5

b1 <- -0.065

b0_h <- -3.5

b1_h <- 0.2

n <- 100

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

y <- r_hpois(n = n,  lambda = exp(b1 * x + b0), prob = 1/(1+exp(b0_h + b1_h*x)))

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

以前紹介したように係数の推定は下記のように行った場合と同等のものが得られる。

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

summary(glmmTMB(y2 ~ x, data = d, family = binomial(link = "logit")))

##  Family: binomial  ( logit )

## Formula:          y2 ~ x

## Data: d

## 

##      AIC      BIC   logLik deviance df.resid 

##     79.9     85.1    -37.9     75.9       98 

## 

## 

## Conditional model:

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

## (Intercept)  5.15921    1.00207   5.149 2.62e-07 ***

## x           -0.29943    0.05762  -5.197 2.03e-07 ***

## ---

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

summary(glmmTMB(y3 ~ x, data = d, family = truncated_poisson(link = "log")))

##  Family: truncated_poisson  ( log )

## Formula:          y3 ~ x

## Data: d

## 

##      AIC      BIC   logLik deviance df.resid 

##    300.8    305.0   -148.4    296.8       57 

## 

## 

## Conditional model:

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

## (Intercept)  2.62316    0.09035  29.033  < 2e-16 ***

## x           -0.08219    0.01006  -8.167 3.16e-16 ***

## ---

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

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

　また、AICの計算は、上記のように個別に解析したAICの和で表すことができる。

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

AIC(fit1)

## [1] 380.679

AIC(glmmTMB(y2 ~ x, data = d, family = binomial(link = "logit"))) + 

AIC(glmmTMB(y3 ~ x, data = d, family = truncated_poisson(link = "log")))

## [1] 380.679

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

・ゼロ過剰モデルによる解析

　次は、下記のようなデータを解析してみよう。とりあえず、ポアソン分布に従うものとして考える。

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

x <- c(12.7, 5, 5, 29.9, 22.8, 1.9, 2.8, 27, 1, 27.9, 1.9,

       11.6, 15, 3.2, 14.6, 28.4, 5.9, 24.9, 24.6, 10.7, 6.6, 

       0.9, 20.1, 17.4, 14.3, 22.2, 2.7, 0.3, 11.8, 12.9, 18.8, 

       7.7, 24.7, 1.5, 23.2, 9.8, 24.7, 15.8, 27.1, 14.7, 0.4, 

       10.1, 26.4, 9.2, 13.5, 21.3, 9.6, 29.9, 23.4, 3.4, 7.8, 

       0.3, 20, 18.8, 27.7, 10, 23.5, 8.9, 14.3, 24.2, 13.6, 1.2,

       21.7, 25.9, 5.9, 19.8, 15.6, 23.8, 22.3, 14.9, 17.4, 26,

       3.1, 16.5, 27.7, 11.6, 27.6, 15.9, 10.3, 20.7, 17.1, 18.6, 

       3.7, 28.8, 24.3, 2.4, 20, 5, 17.8, 13.8, 15.4, 1.9, 24.4, 

       11.1, 11.9, 15.4, 6.4, 17, 25.6, 6.2)

y <- c(7, 8, 11, 0, 0, 14, 13, 0, 17, 0, 12, 8, 1, 11, 0, 0, 

       4, 3, 0, 0, 12, 8, 0, 0, 0, 2, 8, 10, 5, 5, 1, 8, 0, 5,

       0, 4, 0, 0, 0, 0, 14, 5, 3, 9, 2, 0, 2, 4, 0, 5, 0, 10,

       0, 0, 1, 14, 0, 8, 0, 1, 2, 14, 6, 3, 10, 0, 7, 0, 0, 0,

       2, 0, 16, 2, 0, 6, 1, 3, 0, 0, 0, 0, 12, 0, 0, 11, 1, 

       10, 0, 0, 7, 5, 0, 5, 7, 4, 10, 0, 0, 7)

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

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

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

先ほどのよく似たデータである。

　y > 0の値を1に変換して、y2に格納して、0/1データとしてプロットし、0の出現について考えてみる。あとで使うので、y = 0をNAに変換してy3に格納もしておく。

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

d$y2 <- d$y

d$y2[d$y2 > 0] <- 1

d$y3 <- d$y

d$y3[d$y3 == 0] <- NA

plotn(y2 ~ x, data = d)#図9の描画

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

すると、0の出現はxが増えるほど増加する。

　このデータも単純にポアソン回帰すると過分散であることがわかる。

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

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

## [1] 2.474405

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

今回は、練習だと思ってゼロ過剰モデルを想定して解析してみよう。ziformula引数には、ベルヌーイ分布を支配するパラメータpに関する線形予測子として、~ xを入力する。確率分布は普通のポアソン分布とすれば、ゼロ過剰モデルとして計算していることになる。ポアソン回帰の回帰曲線も描いてみよう。

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

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

                ziformula = ~ x) #一般化線形モデル、familyがpoissonであることに注目

fit2 <- summary(fit1)

fit2

##  Family: poisson  ( log )

## Formula:          y ~ x

## Zero inflation:     ~x

## Data: d

## 

##      AIC      BIC   logLik deviance df.resid 

##    385.2    395.6   -188.6    377.2       96 

## 

## 

## Conditional model:

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

## (Intercept)  2.541867   0.073936   34.38   <2e-16 ***

## x           -0.071803   0.008144   -8.82   <2e-16 ***

## ---

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

## 

## Zero-inflation model:

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

## (Intercept) -3.04474    0.64735  -4.703 2.56e-06 ***

## x            0.16326    0.03581   4.560 5.12e-06 ***

## ---

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

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

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

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

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

yy <- exp(b0e + b1e * xx)

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

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

そして、ベルヌーイ回帰の結果は以下のように描ける。

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

b0e_h <- coef(fit2)$zi[1,1]

b1e_h <- coef(fit2)$zi[2,1]

plotn(y2 ~ x, data = d)#図11の描画

xx <- seq(0, 30, length = 200)

yy <- 1/(1+exp(b0e_h + b1e_h * xx))

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

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

また、ポアソン回帰とベルヌーイ回帰の両方の結果を考慮した回帰曲線は以下のように描ける。

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

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

xx <- seq(0, 30, length = 200)

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

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

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

30付近でよりゼロに近づくような期待値となる回帰曲線が描けた。

　そしてこのモデルは過分散が解消されている。

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

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

## [1] 1.238182

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

　実際、今回の解析は下記のように設定してデータを生成したので、悪くない推定になっていることがわかるだろう。

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

b0 <- 2.5

b1 <- -0.065

b0_h <- -3.5

b1_h <- 0.2

n <- 100

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

y <- r_zipois(n = n,  lambda = exp(b1 * x + b0), prob = 1/(1+exp(b0_h + b1_h*x)))

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

ところで今回は下記のように推定しても係数は一致しない。理由は簡単で、y2は構造的なゼロとポアソン分布に由来するゼロを区別していないためである。当然、AICも一致しない。

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

summary(glmmTMB(y2 ~ x, data = d, family = binomial(link = "logit")))

##  Family: binomial  ( logit )

## Formula:          y2 ~ x

## Data: d

## 

##      AIC      BIC   logLik deviance df.resid 

##    103.6    108.8    -49.8     99.6       98 

## 

## 

## Conditional model:

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

## (Intercept)  3.10660    0.64379   4.825 1.40e-06 ***

## x           -0.17174    0.03519  -4.881 1.06e-06 ***

## ---

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

summary(glmmTMB(y3 ~ x, data = d, family = poisson(link = "log")))

##  Family: poisson  ( log )

## Formula:          y3 ~ x

## Data: d

## 

##      AIC      BIC   logLik deviance df.resid 

##    284.0    288.2   -140.0    280.0       57 

## 

## 

## Conditional model:

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

## (Intercept)  2.529966   0.073234   34.55   <2e-16 ***

## x           -0.068764   0.007673   -8.96   <2e-16 ***

## ---

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

AIC(fit1)

## [1] 385.2024

AIC(glmmTMB(y2 ~ x, data = d, family = binomial(link = "logit"))) + 

AIC(glmmTMB(y3 ~ x, data = d, family = poisson(link = "log")))

## [1] 387.6189

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

・ハードルモデルデータをゼロ過剰モデルで解析すると？

　これらの解析を通じて気になるのは、もし、ハードルモデルを仮定したデータをゼロ過剰モデルで解析する、あるいはその逆をするとどうなるかである。まず、下記のようなハードルモデルに従うデータを生成し、ハードルモデルとゼロ過剰モデルを使って解析をしてみる。

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

b0 <- 0.1

b1 <- 0

b0_h <- -0.5

b1_h <- 0

n <- 100

p1_1 <- NULL

p2_1 <- NULL

p1_1_h <- NULL

p2_1_h <- NULL

b0e1 <- NULL

b1e1 <- NULL

b0e1_h <- NULL

b1e1_h <- NULL

p1_2 <- NULL

p2_2 <- NULL

p1_2_z <- NULL

p2_2_z <- NULL

b0e2 <- NULL

b1e2 <- NULL

b0e2_z <- NULL

b1e2_z <- NULL

for (i in 1:1000){

  x <- runif(n, 0, 30)

  y <- r_hpois(n = n,  lambda = exp(b1 * x + b0), 

               prob = 1/(1+exp(b0_h + b1_h*x)))

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

  fit1 <- glmmTMB(y ~ x, data = d, family = truncated_poisson(link = "log"), ziformula = ~ x) #一般化線形モデル

  fit2 <- summary(fit1)

  fit3 <- glmmTMB(y ~ x, data = d, family = poisson(link = "log"), ziformula = ~ x) #一般化線形モデル

  fit4 <- summary(fit3)

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

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

  p1_1_h <- c(p1_1_h, coef(fit2)$zi[1,4])

  p2_1_h <- c(p2_1_h, coef(fit2)$zi[2,4])

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

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

  b0e1_h <- c(b0e1_h, coef(fit2)$zi[1,1])

  b1e1_h <- c(b1e1_h, coef(fit2)$zi[2,1])

  p1_2 <- c(p1_2, coef(fit4)$cond[1,4])

  p2_2 <- c(p2_2, coef(fit4)$cond[2,4])

  p1_2_z <- c(p1_2_z, coef(fit4)$zi[1,4])

  p2_2_z <- c(p2_2_z, coef(fit4)$zi[2,4])

  b0e2 <- c(b0e2, coef(fit4)$cond[1,1])

  b1e2 <- c(b1e2, coef(fit4)$cond[2,1])

  b0e2_z <- c(b0e2_z, coef(fit4)$zi[1,1])

  b1e2_z <- c(b1e2_z, coef(fit4)$zi[2,1])

}

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

実はこれを実行するとすごい量のエラーを吐く。後述するが、ゼロ過剰モデルの推定がうまくいかないことが多発するためである。さて、先に危険率を計算しよう。エラーを吐く場合、p値などがNAと表示される。こういった解析分は除去して計算が必要である。まずはハードルモデル。

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

#ハードルモデル

#エラーを吐いた解析を除いたシミュレーション総数

sum(!is.na(p1_1))

## [1] 1000

sum(!is.na(p2_1))

## [1] 1000

sum(!is.na(p1_1_h))

## [1] 1000

sum(!is.na(p2_1_h))

## [1] 1000

#有意なp値の数

sum(p1_1 < 0.05, na.rm = T)

## [1] 92

sum(p2_1 < 0.05, na.rm = T)

## [1] 50

sum(p1_1_h < 0.05, na.rm = T)

## [1] 234

sum(p2_1_h < 0.05, na.rm = T)

## [1] 43

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

幸い、どのシミュレーションもエラーを出していないようだ（理由は明白で後述する）。ポアソン回帰とベルヌーイ回帰のxの係数の危険率はともに5%程度で、妥当だろう。切片の検出力は9%、23%程度である。p値のヒストグラムを描く。

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

#ハードルモデルのp値

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

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

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

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

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

推定値のヒストグラムと平均値も出す。悪くない推定になっていることがわかるだろう。

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

#ハードルモデル

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

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

histn(b0e1_h, xlab = "b0_h", ylab = "頻度")#図19の描画

histn(b1e1_h, xlab = "b1_h", ylab = "頻度")#図20の描画

#推定値の平均値

mean(b0e1, na.rm = T)

## [1] 0.07212613

mean(b1e1, na.rm = T)

## [1] -0.0002599534

mean(b0e1_h, na.rm = T)

## [1] -0.5305668

mean(b1e1_h, na.rm = T)

## [1] 0.000519859

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

次に問題児のゼロ過剰モデルを確認しよう。

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

#ゼロ過剰モデル

#エラーを吐いた解析を除いたシミュレーション総数

sum(!is.na(p1_2))

## [1] 974

sum(!is.na(p2_2))

## [1] 974

sum(!is.na(p1_2_z))

## [1] 916

sum(!is.na(p2_2_z))

## [1] 902

#有意なp値の数

sum(p1_2 < 0.05, na.rm = T)

## [1] 99

sum(p2_2 < 0.05, na.rm = T)

## [1] 62

sum(p1_2_z < 0.05, na.rm = T)

## [1] 14

sum(p2_2_z < 0.05, na.rm = T)

## [1] 6

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

大体、30～100程度、エラーを出している。ポアソン回帰とベルヌーイ回帰のxの係数の危険率は62/974 × 100 = 6%、6/902 × 100 = 0.6%程度である。ベルヌーイ回帰の危険率は異様に低い。切片の検出力は99/974 × 100 = 10%、14/916 × 100 = 1.5%程度で、同様にベルヌーイ回帰の検出力は異様に低い。p値のヒストグラムを描く。

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

#ゼロ過剰モデルのp値

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

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

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

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

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

すると明らかにベルヌーイ回帰の挙動がおかしい。その理由は推定値のヒストグラムを見ればわかる。p値はNAになっても推定値は算出されていることが多い。

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

#ゼロ過剰モデル

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

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

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

histn(b1e1_h, xlab = "b1_h", ylab = "頻度")#図28の描画

#推定値の平均値

mean(b0e2, na.rm = T)

## [1] 0.1315804

mean(b1e2, na.rm = T)

## [1] -0.00188318

mean(b0e2_z, na.rm = T)

## [1] -93.56835

mean(b1e2_z, na.rm = T)

## [1] -6.281456

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

ポアソン回帰のほうは問題ないが、ベルヌーイ回帰のほうに問題が出ていることがわかる。ベルヌーイ回帰の切片項に関して、異常なほど低い値が検出されている。外れ値の影響で平均値もかなりずれている。切片項に関して異常に大きな負の値が検出されるということは、pの値としては1に極めて近い値となっているということである。ゼロ過剰分布においてp = 1となるときは、ポアソン分布になることを思い出すと、p = 1（実際には1にはなれない）に近づいても過剰に低い切片項を算出し、エラーを出すということは、この時生成された乱数は、ゼロ過少であった可能性が高い。ゼロ過少データは、ゼロ過剰モデルでカバーできないので、モデルを無理やりフィットさせるために、過剰に低い切片項が検出されたと推測できる。ハードルモデルはゼロ過少でも解析可能である。ここからわかるように、ゼロ過剰モデルは解析するデータをよく見ておかないと、解析に失敗する可能性が高い。ちなみにこういう時は中央値を確認すると、マシな値が出てくるのが定石である。

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

#推定値の中央値

median(b0e2, na.rm = T)

## [1] 0.1257326

median(b1e2, na.rm = T)

## [1] -0.001930435

median(b0e2_z, na.rm = T)

## [1] -2.174794

median(b1e2_z, na.rm = T)

## [1] -0.03313515

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

切片項の推定値はまだ過小評価だが、平均値よりはマシな値となっている。

・ゼロ過剰モデルデータをハードルモデルで解析すると？

　今度はゼロ過剰モデルに従ってデータを生成して、ハードルモデルで解析することを考える。ここではもうヒストグラムなどは表示せず、エラーの数や推定値の平均などを眺めるだけにとどめよう。興味がある人は各自ヒストグラムを確認してみるとよい。

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

b0 <- 0.1

b1 <- 0

b0_h <- -0.5

b1_h <- 0

n <- 100

p1_1 <- NULL

p2_1 <- NULL

p1_1_z <- NULL

p2_1_z <- NULL

b0e1 <- NULL

b1e1 <- NULL

b0e1_z <- NULL

b1e1_z <- NULL

p1_2 <- NULL

p2_2 <- NULL

p1_2_h <- NULL

p2_2_h <- NULL

b0e2 <- NULL

b1e2 <- NULL

b0e2_h <- NULL

b1e2_h <- NULL

for (i in 1:1000){

  x <- runif(n, 0, 30)

  y <- r_zipois(n = n,  lambda = exp(b1 * x + b0), 

               prob = 1/(1+exp(b0_h + b1_h*x)))

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

  fit1 <- glmmTMB(y ~ x, data = d, family = poisson(link = "log"), ziformula = ~ x) #一般化線形モデル

  fit2 <- summary(fit1)

  fit3 <- glmmTMB(y ~ x, data = d, family = truncated_poisson(link = "log"), ziformula = ~ x) #一般化線形モデル

  fit4 <- summary(fit3)

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

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

  p1_1_z <- c(p1_1_z, coef(fit2)$zi[1,4])

  p2_1_z <- c(p2_1_z, coef(fit2)$zi[2,4])

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

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

  b0e1_z <- c(b0e1_z, coef(fit2)$zi[1,1])

  b1e1_z <- c(b1e1_z, coef(fit2)$zi[2,1])

  p1_2 <- c(p1_2, coef(fit4)$cond[1,4])

  p2_2 <- c(p2_2, coef(fit4)$cond[2,4])

  p1_2_h <- c(p1_2_h, coef(fit4)$zi[1,4])

  p2_2_h <- c(p2_2_h, coef(fit4)$zi[2,4])

  b0e2 <- c(b0e2, coef(fit4)$cond[1,1])

  b1e2 <- c(b1e2, coef(fit4)$cond[2,1])

  b0e2_h <- c(b0e2_h, coef(fit4)$zi[1,1])

  b1e2_h <- c(b1e2_h, coef(fit4)$zi[2,1])

}

#ゼロ過剰モデル

#エラーを吐いた解析を除いたシミュレーション総数

sum(!is.na(p1_1))

## [1] 996

sum(!is.na(p2_1))

## [1] 996

sum(!is.na(p1_1_z))

## [1] 996

sum(!is.na(p2_1_z))

## [1] 995

#有意なp値の数

sum(p1_1 < 0.05, na.rm = T)

## [1] 100

sum(p2_1 < 0.05, na.rm = T)

## [1] 98

sum(p1_1_z < 0.05, na.rm = T)

## [1] 19

sum(p2_1_z < 0.05, na.rm = T)

## [1] 21

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

ベルヌーイ回帰の検出力は低いが、エラーを吐く確率はずいぶん改善されている。

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

#推定値の平均値

mean(b0e1, na.rm = T)

## [1] 0.08019106

mean(b1e1, na.rm = T)

## [1] -0.00170989

mean(b0e1_z, na.rm = T)

## [1] -10.26869

mean(b1e1_z, na.rm = T)

## [1] -0.5542752

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

やはり、外れ値の影響を受けてベルヌーイ回帰の推定は怪しい部分がある。中央値も確認しよう。

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

median(b0e1, na.rm = T)

## [1] 0.1115595

median(b1e1, na.rm = T)

## [1] -0.00242938

median(b0e1_z, na.rm = T)

## [1] -0.521971

median(b1e1_z, na.rm = T)

## [1] -0.001466032

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

中央値で見れば推定は、「平均的に」うまくいっていることがわかる。ハードルモデルも確認しよう。ハードルモデルはエラーを吐いていない。

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

#ハードルモデル

sum(!is.na(p1_2))

## [1] 1000

sum(!is.na(p2_2))

## [1] 1000

sum(!is.na(p1_2_h))

## [1] 1000

sum(!is.na(p2_2_h))

## [1] 1000

sum(p1_2 < 0.05, na.rm = T)

## [1] 75

sum(p2_2 < 0.05, na.rm = T)

## [1] 51

sum(p1_2_h < 0.05, na.rm = T)

## [1] 113

sum(p2_2_h < 0.05, na.rm = T)

## [1] 40

mean(b0e2, na.rm = T)

## [1] 0.07243546

mean(b1e2, na.rm = T)

## [1] -0.001079898

mean(b0e2_h, na.rm = T)

## [1] 0.3236617

mean(b1e2_h, na.rm = T)

## [1] 0.001452809

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

基本的にハードルモデルとゼロ過剰モデルはよく類似しているし、パラメータの推定値も似た値が算出されることが多い。なので、データの構造や、生成のメカニズムに応じて選択することになるだろう。構造的なゼロが出現することが仮定できるときはゼロ過剰モデルを選ぶべきだし、ゼロ過少ならハードルモデルでもよいだろう。