分布の混合再来

　前回、とうとうGLMMの真の意味である、固定効果とランダム効果を考慮した解析を紹介した。GLMMは、その名前から混乱を招きやすいが、分布の混合が名前の由来ではない。しかし、分布の混合と密接なかかわりを持っている。さて、今回と次回、再度、分布の混合を取り扱う。解析の中身も、GLMに一度戻そうと思う。また、タイトル詐欺になってしまうが、やはり分布の混合はGLMMあたりで解説するのが収まりが良いと思うので、許してほしい。

0を含む正の連続値データ

　さて、今回扱うデータは、「0を含む」正の連続値をとるようなデータである。正の連続値なら、ガンマ分布とかベータ分布とか使えばいいのではないか、と考えるかもしれないが、実は全く一筋縄ではいかない。というのは、ガンマ分布もベータ分布も、0を台に含まない、つまり0の確率密度が定義されていないのである。言い換えれば、ガンマ分布やベータ分布に従うと仮定した時点で、0は出現しないものとして扱われる。仮にデータの中に0を含んでいて、ガンマ回帰などを行おうとすると、0はガンマ分布では定義されていないので、解析できません、と文句を言われてしまう。

　この問題、シンプルながらなかなか厄介で、しかもデータ解析を行っているとしばしば立ちはだかる問題である。例えば、ある処理に従って変化する植物の重さをモデリングすることを考えよう。この時、植物の重さが0に近づくと、この問題は発生しうる。そもそも植物の重さは、なにか計量計を使って量るのが普通だろう。そして、大抵は、量れる下限値が決まっている。植物としてはサンプリングされているので、重さがないわけではないが、下限値を下回ってしまうと0と表示されてしまう。

重さ（>0）のデータで、0付近をモデリングするのだからガンマ回帰が適切なはずなのに、0が出てしまうために、使えないのだ。これだけで、非常にいやらしい問題だとわかるだろう。解決方策はいくつかある。1～4つ目まで、順に高度な技術が求められる。

　1つ目、0データは欠損値として扱うことにする。せっかく取れたデータだが使えないのであきらめるということだ。けれども、0データが多かったら？　0が生じるメカニズムがわかっているのに、それを全部、使えないとするのはあまりにもったいないし、解析が失敗してしまうことを意味するだろう。

　2つ目、0にごく小さい値を足す。これで0を消すということだ。これで余すことなくデータは使える。しかし、実際解析したことがある立場から言うなら、推奨はされない。というのも、足す値によって結果が劇的に変わるからだ。10^-10と10^-9、どちらもごくごく小さな値だが、どちらを足すかで結果はかなり変わる。そして、どんな値を足すかは恣意的に決めざるを得ない。これは、ガンマ分布やベータ分布が0付近で劇的に確率密度が変わるため生じる問題で、結局、これらの確率分布を使う以上、避けられない問題なのだ。

　3つ目が、今回紹介するハードルモデルだ。詳細は後述するが、これを使えば、0を0のまま使っても無理なく解釈が可能だ。

　4つ目、これが最も高度で、0が生じるメカニズムを最も正確に記述できる。計量計の量れる重さの下限値以下の確率分布を積分値としてあらわすような尤度関数を定義し、それをもとにベイズ推定する。ベイズ推定は全く扱ってきていないので、立ち入りを避けるが、興味がある人は松浦健太郎氏の「StanとRでベイズ統計モデリング」で取り上げられる、打ち切りデータのモデリングの項を参考にしてほしい。

ハードルモデル：閾値を乗り越えれば正の値の確率分布に従う

　ハードルモデルhurdle modelは0が過剰に生じるようなデータで採用されるモデリング手法の一つである。その確率分布は、ベルヌーイ分布と正の値をとる確率分布の混合分布として表現される。今回は、ガンマ分布との混合を考えよう。ハードルモデルの発想は、何か見えない閾値＝ハードルが存在し、そのハードルを乗り越えた値は正の値として検出されるというものである。前述の、植物の重さなんかは良い例で、ここでのハードルは計量計の計量下限値である。逆に、そのハードルを越えられなければ、0として検出される。

　このハードルを超えられるか否かは、確率pに従うベルヌーイ分布に従うと考える。ベルヌーイ分布は0か1をとる確率分布で、今回は、1-pで0、pで1を生成するものとする。もし、0を生成すれば、それは0が観測された場合と同義である。一方、1を生成した場合は、ハードルを乗り越えたと考え、この値はガンマ分布に従う。

ハードルモデルの確率密度関数は下記のように定義できる。ここで、ベース分布とはベルヌーイ分布と混合する確率分布であり、θをパラメータとすると考えている。

式の形はいかついが、主張はシンプルである。x = 0を生成する確率は1-pであり、x > 0が生成する確率はpになるように調整がなされている。ベース分布が0を生成する確率g(0|θ)だけ差し引き、それで規格化したのちにpを掛け、x > 0の範囲で積分 or 和をとればx > 0の範囲で確率pとなる。このような、もとの確率分布から0を生成する確率を差し引いて規格化した分布をゼロ切断分布zero-truncated distributionと呼ぶ。連続分布の場合、あまり意味のない表現だが、離散分布では重要な概念である。

　ベース分布がガンマ分布の場合、g(0|θ) = 0だから（正確には不正確な表現だけど）、下記のようにあらわせる。上記のθはガンマ分布の具体的なパラメータ、kとλに置き換えた。

　さて、GLMやGLMMでは、様々な確率分布を扱えること以外に、線形予測子とリンク関数という重要な概念があった。この確率分布には、どのように線形予測子とリンク関数を組み込むのだろうか。実は、ハードルモデルは2つの線形予測子とリンク関数のセットを引数として持つことができる。ハードルモデルの線形予測子やリンク関数を考慮したデータの生成メカニズムは下記のような構造になっている。

　まず、上段、ベルヌーイ分布を支配するpに関して、線形予測子とリンク関数を設定できる。とはいえ、確率を扱うので、ここのリンク関数はロジットリンク一択である。y > 0の値をすべて1に変換すると、yは0と1のみになる。これに対して、二項（ベルヌーイ）回帰を行って、パラメータ群βhを推定する。

　続いて、下段ではガンマ分布を支配するλに関して、ベルヌーイ分布と独立にさらに線形予測子とリンク関数を設定できる。ここでは対数リンクとした。y = 0となるデータを除き、これに対してガンマ回帰を行うことでパラメータ群βを推定する。

　推定されるパラメータはβhとβの合計数であり、ベルヌーイ分布とガンマ分布に異なる線形予測子が組み込まれていてもよい。例えば、0は一定確率で生じているように見えるが、それ以外の値はxに対して増加するような傾向が見えるなら、ベルヌーイ分布には切片項だけを推定させ、ガンマ分布にはxを説明変数に含むように設定してもよい。このようにデータを2つの部分に分けて解析を行うので、ハードルモデルや次回紹介するゼロ過剰モデルはtwo-parts modelとも呼ばれる。

　ガンマ回帰とベルヌーイ回帰の両方を考慮した回帰曲線の推定値は、2つの推定値の積なので、以下のようにあらわすことができる。

　では、このような前提知識の下、データを解析してみよう。

具体的な解析

　例えば、ゼロ過剰なガンマ分布に従うデータは下記のように定義した関数から生成することができる。あらかじめ、線形予測子とリンク関数を想定した入力にしている。

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

library(plotn)

library(glmmTMB)

r_hgamma <- function(n, shape, rate = 1, prob){

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

  v <- mapply(function(x, shape, rate){

    if(x == 1) {

      rgamma(n = 1, shape = shape, rate = rate)

      } else {0}

    },

    x = tmp, shape = shape, rate = rate

    )

  unlist(v)

}

b0 <- 1

b0_h <- -0.5

shape <- 2

n <- 100

r_hgamma(n = n,  shape = shape, rate = shape/exp(b0), prob = 1/(1+exp(b0_h)))

##   [1] 3.1051700 1.0132823 4.7937589 3.3121354 0.0000000 2.4788693 1.6439880

##   [8] 5.7772205 1.6258768 0.0000000 3.0275223 2.0914953 2.2799799 3.6170907

##  [15] 0.0000000 5.6925916 2.1958523 1.6221110 4.9804616 0.5833789 0.0000000

##  [22] 0.5894626 0.0000000 2.3122172 1.0176533 7.3133324 0.0000000 2.2037698

##  [29] 0.0000000 0.0000000 0.0000000 3.4071564 0.0000000 1.0079533 0.0000000

##  [36] 2.7385534 0.0000000 0.0000000 0.9846316 0.0000000 0.0000000 2.7380465

##  [43] 2.1012453 2.7935054 3.5903841 0.0000000 1.6448305 0.0000000 3.3843574

##  [50] 3.1017506 4.3335586 0.0000000 2.8251490 0.0000000 4.7249172 1.4138173

##  [57] 1.1615968 0.0000000 4.0407703 6.8632643 2.3200155 2.2094785 0.0000000

##  [64] 0.6610596 0.5893749 0.0000000 0.0000000 0.0000000 2.9864599 1.3785199

##  [71] 1.8645623 2.6062540 4.2006592 0.0000000 0.3653368 3.1804218 0.0000000

##  [78] 0.6429347 0.0000000 0.0000000 1.0231629 1.5742456 5.6903338 0.0000000

##  [85] 0.0000000 1.0750316 0.0000000 0.0000000 4.5342978 2.6480068 0.0000000

##  [92] 0.0000000 2.9370420 2.7465387 1.1259302 0.0000000 0.0000000 3.2167045

##  [99] 3.1506699 0.0000000

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

ガンマ分布の時には生成されることがなかった、0が生成されている。上記の例では、1-1/(1+exp(-0.5)) = 0.378くらいの確率で0が生成され、残りがガンマ分布に従うように生成される。もっと大量に生成して、分布の様子を見てみよう。

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

n <- 10000

x <- r_hgamma(n = n,  shape = shape, rate = shape/exp(b0), prob = 1/(1+exp(b0_h)))

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

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

確かに、0がたくさん生成され、ガンマ分布とは異なる様子がわかる。

　では、下記のようなデータを解析してみよう。

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

x <- c(17.6, 10.5, 8.7, 27.3, 11.6, 16.1, 28.7, 7.6, 8.9, 

       21.2, 4.3, 5.5, 11.1, 13.2, 11.2, 8, 28.3, 25, 26.2,

       15.2, 11.1, 29.5, 13.9, 11.8, 17.4, 28.8, 1.3, 2.5,

       19.9, 17.9, 23.8, 26.7, 20.1, 10, 13.6, 5.4, 28.4, 26.1,

       21, 17.2, 8.6, 17.1, 1.5, 28.2, 4.5, 10.6, 6.8, 26.2,

       26.8, 16.3, 28.7, 23.9, 11.2, 21.1, 20.7, 1.8, 3.6, 25.5,

       20.9, 15.6, 23.5, 29.7, 19.6, 24.9, 15.2, 2.3, 19.2,

       17.3, 7.4, 12.2, 14.2, 8.2, 25.8, 2.5, 15.2, 4.4, 15.4,

       4.5, 26.9, 6.7, 0.3, 28.6, 4.7, 21.2, 28.6, 17.4, 25.5,

       16.8, 15.7, 7, 5.9, 16.1, 22, 23, 14.1, 15.5, 9.3, 17.6, 16, 1.3)

y <- c(0, 7.15, 0, 0, 2.92, 0, 0, 0, 0, 1.75, 1.2, 1.85, 1.94,

       0.77, 0, 0.3, 0, 0, 0, 0.58, 2.62, 4.65, 6.67, 0, 0.6, 

       0, 0.14, 1.85, 10.83, 0.1, 0, 2.84, 1.17, 4, 4.24, 0, 1.24, 

       1.55, 4.66, 0, 2.39, 3.11, 11.53, 0, 3.28, 0, 0, 7.15, 0.79,

       2.1, 0, 1.18, 3.49, 0, 4.39, 1.32, 6.52, 0, 1.77, 0, 1.11,

       4.29, 0, 1.26, 3.5, 0.89, 0.8, 0, 7.78, 1.8, 3.58, 14.48, 

       0, 0, 0, 4.86, 0.55, 6.91, 0, 2.42, 3.16, 4.11, 1.45, 2.67, 

       2.01, 3.41, 1.33, 1.37, 0.36, 0, 5.66, 1.79, 0, 1.78, 1.29, 

       3.34, 1.31, 4.7, 1.76, 0)

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

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

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

データはxとはあまり関係がなさそうな正の連続値と0を示すデータである。本当ならガンマ回帰したいところだが、0が存在するため、普通に解析しようとすると下記のようにエラーを出す。

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

glm(y ~ x, data = d, family = Gamma(link = "log"))

## Error in eval(family$initialize) : 

##   non-positive values not allowed for the 'Gamma' family

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

そこで、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)#図3の描画

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

すると、0の出現はxとはあまり関係がなく、一定の割合で出現しているように見える。そこで、ハードルモデルを実行して解析してみよう。glmmTMBパッケージのglmmTMB関数はハードルモデルを扱える。まず、ガンマ回帰するときのglm関数などと同じように表記する。そしてziformula引数には、ベルヌーイ分布を支配するパラメータpに関する線形予測子を入れる。今回は、ベルヌーイ分布のpはxとは関係がなさそうなので、~ 1として切片項だけ推定させよう。確率分布はガンマ分布とするが、通常のガンマ分布だとエラーを返されるので、0を無視するように設定して作られたglmmTMBパッケージのziGamma分布を使うことにする。

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

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

fit2 <- summary(fit1)

fit2

##  Family: Gamma  ( log )

## Formula:          y ~ x

## Zero inflation:     ~1

## Data: d

## 

##      AIC      BIC   logLik deviance df.resid 

##    415.9    426.3   -203.9    407.9       96 

## 

## 

## Dispersion estimate for Gamma family (sigma^2): 0.672 

## 

## Conditional model:

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

## (Intercept)  1.35774    0.19956   6.804 1.02e-11 ***

## x           -0.01604    0.01168  -1.373     0.17    

## ---

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

## 

## Zero-inflation model:

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

## (Intercept)  -0.7538     0.2144  -3.516 0.000438 ***

## ---

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

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

Conditional modelの項は、今までも解説したとおりである。目新しい項目はZero-inflation modelの項目で、ここでベルヌーイ回帰の結果が出ている。結果の読み取り方は、図を見たほうが早いだろう。回帰曲線を描いてみる。まずはガンマ回帰の結果である。以下のようにすることでConditional modelの項の係数を取得できる。

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

plotn(y ~ x, data = d)

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

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

xx <- seq(min(x), max(x), length = 200)#図4の描画

yy <- exp(b0e + b1e * xx)

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

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

確かに、xとはあまり関係がない様子がわかるだろう。そして、ベルヌーイ回帰の結果は以下のように描ける。下記のようにすればzero-inflation modelの係数を取得できる。

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

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

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

xx <- c(min(x), max(x))

yy <- rep(1/(1+exp(b0e_h)),2)

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

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

当然、切片項だけを推定したので、水平線が引けるだけである。

　そして、2つの推定値を考慮した回帰式は以下のようになる。

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

plotn(y ~ x, data = d)

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

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

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

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

ベルヌーイ回帰の線形予測子の結果から、ガンマ回帰の結果を一定割合で低下させた結果となっている。。

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

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

b0 <- 1

b1 <- 0

b0_h <- -0.5

b1_h <- 0

shape <- 2

n <- 100

x <- runif(n, 0, 30)

y <- r_hgamma(n = n,  shape = shape, rate = shape/exp(b1 * x + b0), prob = 1/(1+exp(b0_h + b1_h*x)))

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

実は係数の推定は下記のように行った場合と同等のものが得られる。まずzero-inflation modelのほうはベルヌーイ回帰だったので、0/1に変換したy2とglm関数を使うと下記の通り。

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

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

##  Family: binomial  ( logit )

## Formula:          y2 ~ 1

## Data: d

## 

##      AIC      BIC   logLik deviance df.resid 

##    127.4    130.0    -62.7    125.4       99 

## 

## 

## Conditional model:

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

## (Intercept)   0.7538     0.2144   3.516 0.000438 ***

## ---

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

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

正負は逆転しているが同じ値であることがわかる。これは、線形予測子を入れるときに普通のベルヌーイ回帰だと負の値にして入力するが、zero-inflation modelではそうしないために正負が逆転する。P値なども同じ値だ。

　続いて、ガンマ回帰は、0の値を除いたy3とglm関数を使うと下記の通り。

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

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

##  Family: Gamma  ( log )

## Formula:          y3 ~ x

## Data: d

## 

##      AIC      BIC   logLik deviance df.resid 

##    288.5    295.2   -141.2    282.5       65 

## 

## 

## Dispersion estimate for Gamma family (sigma^2): 0.672 

## 

## Conditional model:

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

## (Intercept)  1.35774    0.19956   6.804 1.02e-11 ***

## x           -0.01604    0.01168  -1.373     0.17    

## ---

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

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

　こちらも同じ値になっていることがわかるだろう。

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

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

AIC(fit1)

## [1] 415.8693

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

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

## [1] 415.8693

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

AICを計算できるということは、説明変数の選択が可能ということであり、ガンマ回帰、ベルヌーイ回帰、両方に関して、より予測性の高いモデルを構築することが可能である。

　では、上記の設定でデータを1000回生成して、検出力を確認してみる。glmmTMBは計算がやや遅いので、1000回程度でとどめる。

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

b0 <- 1

b1 <- 0

b0_h <- -0.5

shape <- 2

n <- 100

p1 <- NULL

p2 <- NULL

p1_h <- NULL

b0e <- NULL

b1e <- NULL

b0e_h <- NULL

for (i in 1:1000){

  x <- runif(n, 0, 30)

  y <- r_hgamma(n = n,  shape = shape, 

                rate = shape/exp(b1 * x + b0), 

                prob = 1/(1+exp(b0_h)))

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

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

  fit2 <- summary(fit1)

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

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

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

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

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

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

}

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

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

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

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

危険率や検出力は下記のとおりである。

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

#ガンマ回帰

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

## [1] 999

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

## [1] 74

#ベルヌーイ回帰

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

## [1] 718

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

ガンマ回帰に関して、xの係数の危険率がやや高い感じがするが、これがシミュレーションの数が少ないことに起因するかは不明である。おそらくは問題ない範疇であろうと推測する。

　また、各係数については下記の通り。

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

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

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

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

#ガンマ回帰

mean(b0e, na.rm = T)

## [1] 0.9888244

mean(b1e, na.rm = T)

## [1] -0.00020859

#ベルヌーイ回帰

mean(b0e_h, na.rm = T)

## [1] -0.5063652

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

いずれも、推定値の平均値は設定どおりになっていることがわかるだろう。

　次は以下のようなデータについて解析してみよう。

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

x <- c(20.3, 26.1, 17.6, 6.6, 24.7, 7.9, 3.5, 11.6, 7.9, 21.3, 

       9.7, 12.9, 28.1, 16.3, 15.2, 25.3, 13.6, 2.1, 0.3, 26.5,

       0.7, 10.4, 16, 3.6, 29.3, 12, 20.9, 28.1, 23.1, 25.8, 9.1, 

       13.2, 18.9, 16.6, 0.3, 18, 9.8, 1.1, 11.3, 4.1, 3.6, 15, 

       8.2, 23.8, 7, 21.8, 7.4, 4.2, 12.3, 20.7, 16.5, 2.1, 5.9, 

       2.9, 5.6, 22.7, 26.5, 22.3, 9.2, 7.4, 1.7, 18.6, 9.1, 2.5,

       19.6, 11.9, 7.5, 25.9, 1.1, 20.5, 5, 14.1, 22.4, 10, 15.9, 

       15.2, 21.8, 28.2, 26.4, 4.6, 5.9, 21.6, 11.6, 21.9, 29.8, 

       6.4, 22, 1.1, 24.2, 12.6, 26.5, 22.1, 26.6, 23.8, 5.5, 2.6,

       4.2, 19.3, 8.5, 6.7)

y <- c(12.08, 24.79, 6.79, 0, 21.08, 4.6, 7.08, 3.08, 6.09, 

       2.83, 0, 0, 8.07, 2.14, 1.61, 6.74, 12.34, 0, 2.15,

       8.37, 2.12, 0, 9.55, 0, 6.14, 0, 3.98, 38.71, 0, 9.86,

       0, 1.69, 3.07, 2.96, 0, 0, 0, 0, 4.08, 5.39, 4.66, 14.13,

       3.68, 4.72, 0, 8.96, 12.16, 0, 0, 7.01, 11.48, 0, 0, 0, 

       0, 3.06, 15.42, 8.49, 0, 3.67, 0.91, 0, 0, 0, 13.79, 9.93,

       1.74, 13.25, 4.38, 0, 0, 12.58, 8.58, 0, 0, 7.27, 2.08, 

       13.02, 4.46, 0.87, 0, 8.13, 1.61, 2.26, 14.38, 1.33, 15.57,

       0, 5.52, 0, 12.63, 10.44, 12.28, 4.85, 0, 0, 2.48, 9.22, 8.68, 0)

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

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

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

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)#図14の描画

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

すると、xが増加するとともに、y > 0となる出現率が増加しているように見える。こういう時は、ziformula引数を ~ ｘとすることでベルヌーイ回帰で切片項とxの係数を推定してくれる。先にガンマ回帰の回帰曲線を描いてしまおう。

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

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

fit2 <- summary(fit1)

fit2

##  Family: Gamma  ( log )

## Formula:          y ~ x

## Zero inflation:     ~x

## Data: d

## 

##      AIC      BIC   logLik deviance df.resid 

##    482.8    495.8   -236.4    472.8       95 

## 

## 

## Dispersion estimate for Gamma family (sigma^2): 0.408 

## 

## Conditional model:

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

## (Intercept) 1.129165   0.175778   6.424 1.33e-10 ***

## x           0.049990   0.009299   5.376 7.61e-08 ***

## ---

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

## 

## Zero-inflation model:

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

## (Intercept)  1.07043    0.43109   2.483    0.013 *  

## x           -0.14031    0.03299  -4.253 2.11e-05 ***

## ---

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

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

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"))

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

そして下記のようにすればzero-inflation modelの係数を取得でき、ベルヌーイ回帰の結果も図示できる。

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

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

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

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

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)#図17の描画

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

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

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

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

ガンマ回帰単独と比べて、ゼロ付近でより低い値が生成されるような推定になっている。

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

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

b0 <- 1

b1 <- 0.06

b0_h <- 1

b1_h <- -0.15

shape <- 2

n <- 100

x <- runif(n, 0, 30)

y <- r_hgamma(n = n,  shape = shape, rate = shape/exp(b1 * x + b0), prob = 1/(1+exp(b0_h + b1_h*x)))

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

係数の推定は下記のように行った場合と同等のものが得られ、それらのAICの和がハードルモデルの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 

##    107.9    113.1    -52.0    103.9       98 

## 

## 

## Conditional model:

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

## (Intercept) -1.07043    0.43109  -2.483    0.013 *  

## x            0.14031    0.03299   4.253 2.11e-05 ***

## ---

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

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

##  Family: Gamma  ( log )

## Formula:          y3 ~ x

## Data: d

## 

##      AIC      BIC   logLik deviance df.resid 

##    374.9    381.4   -184.4    368.9       63 

## 

## 

## Dispersion estimate for Gamma family (sigma^2): 0.408 

## 

## Conditional model:

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

## (Intercept) 1.129165   0.175778   6.424 1.33e-10 ***

## x           0.049990   0.009299   5.376 7.61e-08 ***

## ---

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

AIC(fit1)

## [1] 482.8088

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

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

## [1] 482.8088

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

　では、上記の設定でデータを1000回生成して、検出力を確認してみる。glmmTMBは計算がやや遅いので、1000回程度でとどめる。

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

b0 <- 1

b1 <- 0.06

b0_h <- 1

b1_h <- -0.15

shape <- 2

n <- 100

p1 <- NULL

p2 <- NULL

p1_h <- NULL

p2_h <- NULL

b0e <- NULL

b1e <- NULL

b0e_h <- NULL

b1e_h <- NULL

for (i in 1:1000){

  x <- runif(n, 0, 30)

  y <- r_hgamma(n = n,  shape = shape, 

                rate = shape/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 = ziGamma(link = "log"), ziformula = ~ x) #一般化線形モデル

  fit2 <- summary(fit1)

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

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

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

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

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

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

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

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

}

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

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

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

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

#ガンマ回帰

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

## [1] 994

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

## [1] 1000

#ベルヌーイ回帰

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

## [1] 583

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

## [1] 1000

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

　また、各係数については下記の通り。

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

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

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

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

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

#ガンマ回帰

mean(b0e, na.rm = T)

## [1] 0.9811484

mean(b1e, na.rm = T)

## [1] 0.06047969

#ベルヌーイ回帰

mean(b0e_h, na.rm = T)

## [1] 1.050159

mean(b1e_h, na.rm = T)

## [1] -0.1570608

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

確かに、推定値の平均値は設定したものと同じになることが確認できる。