目次

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

・負の二項分布：0以上の離散値をとる確率分布

・負の二項分布はガンマ分布とポアソン分布の混合分布

・負の二項分布を使ったGLM（負の二項回帰）

・検出力

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

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

　ポアソン分布に対して過分散なデータを無理やりポアソン分布に当てはめると危険率が大きく上昇することを紹介した。それに対して、一般化線形混合モデル(GLMM)における分布の混合という考え方を導入することで、より分散が大きく、柔軟性の高い確率分布を記述し、解析できることを以前紹介した。その際、分布の混合の特別な場合として、ガンマ分布とポアソン分布の混合を紹介し、これらの混合分布が負の二項分布になることに軽く触れた。

　この性質は、過分散なデータの解析おいては、必ずしもランダム効果をもたらすデータの構造を考慮した解析をしなくても解決できうることを示している。つまり、一般化線形モデル(GLM)であっても、負の二項分布の仮定した解析をすれば、解決できるだろう、ということである。本項と、次項では、GLMでありながらも、分布に過分散な確率分布を指定することで解決を図る方法を紹介する。なので、タイトル詐欺気味ではあるのだが、負の二項分布とベータ-二項分布が分布の混合によって得られる恩恵を、わかりやすく示してくれているので、あえてここで紹介させてもらう。本当のGLMMにともかくチャレンジしたい人は次々回に期待してほしい。

負の二項分布：0以上の離散値をとる確率分布

　ポアソン分布でモデリングするには、過分散なデータに使える確率分布が、負の二項分布negative binomial distributionである。この分布は0以上の離散値を台とする分布である。負の二項分布は2つのパラメータ、サイズ（成功回数）パラメータrと成功確率パラメータpで特徴づけられる。「負の二項分布」は二項分布と名前も似ているし、下記のように確率質量関数Probability mass functionも実際類似している。

　負の二項分布の定義から、「成功確率pの事象をr回成功するまでにx回失敗する」とは、「(r+x-1)回の試行回数で(r-1)回成功し、 x回失敗する」かつ「(r+x)回目の試行で成功する」必要があるため、以下のように分布を計算、定義できる。

　二項分布は試行回数パラメータnと成功確率パラメータpを持つ。二項分布は「成功確率pの事象をn回繰り返したとき、1回成功する確率、2回成功する確率、……x回成功する確率、……、n回成功する確率」を記述している。二項分布は総試行回数nを固定している、ととらえることもできるだろう。

　対して、負の二項分布は「成功確率pの事象をr回成功するまでの失敗回数xの確率」を記述する。負の二項分布では総試行回数はr+x回になるのでxによって変動する、ととらえることもできる。

　具体的に、サイコロ投げを考えよう。1の出目が出ることを成功とすれば、成功確率は1/6である。ここで10回成功するまで、サイコロを投げ続けることにする。つまり、1の出目が合計で10回出ればそこで試行を止める。止まった時の失敗回数の分布が負の二項分布に従うということである。豪運なら10回投げて10回とも1の出目が出ることもある。そのときは失敗回数は0で、めったに起きないだろうな、と思えるだろう。このように11回投げる確率（1回失敗）、12回投げる確率（2回失敗）、を記述していくのである。中ほど位がたぶん確率が高くて、逆に失敗回数がかさみすぎる（例えば、1000回投げる羽目になる（失敗回数990回））ことは相当不運だろうな、確率が低いだろうな、という予測が立つ。

　二項分布ともども、我々の人間生活を記述する中で大変身近なテーマを扱える確率分布であり、実用性が高い。興味があれば、負の二項分布を使ったゲーム内での確率論の記事もあるので見てみてほしい。

　では、この定義から平均と分散を計算してみよう。

pが1に近づくほど、平均と分散はともに小さくなる。また、rが大きくなるほど、平均と分散はともに大きくなる。なお、一般化線形モデルにおいては、負の二項分布のサイズパラメータr、あるいはオッズ比(1-p)/pがdispersion parameterとして推定される。

　続いて、具体的に乱数を生成してデータの様子を確認しよう。rbinom関数を使うことで、負の二項分布に従うデータを生成できる。引数nは生成する個数、sizeはサイズパラメータ、probは成功確率パラメータである。なお、通常のglm関数ではベータ分布を使った解析をできないので、ここでglmmTMBパッケージをロードする（ほかには、MASSパッケージのglm.nb関数を利用する手もある）。

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

library(plotn)

library(viridis)

library(glmmTMB)

library(dplyr)

rnbinom(n = 30, size = 10, prob = 0.5)

##  [1]  4 12 10  8 13  8  8  3  7  5 13  3  4 13  2 13  7  8 12 10  5  8  9  8  6

## [26]  8 10  5 11 16

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

このように乱数を生成してみるとわかるが、すべて0以上の整数であることがわかるだろう。続いて、確率パラメータによって確率分布がどうなるかをいくつか確認してみる。試しに、size = 10で固定し、probには0.3, 0.5, 0.9を指定してみよう。

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

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

dd$prob01 <- dnbinom(dd$X, size = 10, prob = 0.3)

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

dd$prob07 <- dnbinom(dd$X, size = 10, prob = 0.9)

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

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

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

pが1に近づくにつれて左よりの分布になることが確認できる。次は、確率パラメータを固定して、さまざまなsizeを指定してみよう。

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

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

dd$size03 <- dnbinom(dd$X, size = 3, prob = 0.5)

dd$size09 <- dnbinom(dd$X, size = 9, prob = 0.5)

dd$size15 <- dnbinom(dd$X, size = 15, prob = 0.5)

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

      leg.lab = c(3, 9, 15), leg.sp = 5, ylab = "確率質量分布")#図2の描画

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

サイズパラメータが大きくなるほど、平均と分散が大きくなり、右寄りの分布になることがわかる。

負の二項分布はガンマ分布とポアソン分布の混合分布

　負の二項分布の重要な性質として、この分布はガンマ分布とポアソン分布の混合分布で表されるという点である。この性質があるために、あえて、一般化線形混合モデルの項で負の二項分布を紹介している。下記のような設定で、ポアソン分布、ガンマ分布、負の二項分布を生成し、分布の様子を見てみる。さらに、手作業でガンマ分布とポアソン分布を混合して生成した乱数も用意する。これらの分布を比較しよう。後述するが、負の二項分布のパラメータをsizeとprobとしたとき、ガンマ分布について形状母数shape = size、比母数rate = prob/(1-prob)と設定し、定義通りにポアソン分布と混合すれば、パラメータをsizeとprobをもつ負の二項分布になる。

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

size <- 6

prob <- 0.2

n <- 1000000

lambda <- size * (1 - prob)/prob

rnb <- rnbinom(n = n, size = size, prob = prob)#負の二項分布の乱数

rp <- rpois(n = n, lambda = lambda)#ポアソン分布の乱数

rg <- rgamma(n = n, shape = size, rate = prob/(1 - prob))#ガンマ分布の乱数

rbinom_by_mix <- function(n, size, prob){

  sapply(rgamma(n = n, shape = size, rate = prob/(1 - prob)), function(tmp) rpois(1, lambda = tmp))

}#手作業での分布の混合

rnb_h <- rbinom_by_mix(n = n, size = size, prob = prob)#混合分布

d <- data.frame(value = c(rnb, rp, rg, rnb_h), 

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

                                 rep("poisson", n),

                                 rep("gamma", n),

                                 rep("mixed", n)))

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

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

すると、もともと用意されていた負の二項分布を生成する関数と、自分で分布を混合した時、非常に近い分布になることがわかる。ポアソン分布は負の二項分布に対して分散が小さいこともわかる。負の二項分布はガンマ分布と非常に類似しているため、区別が難しい。そこで、もっと正確に分布の様子を測る方法で確認する。

　以下のように、手元の乱数から、ある値以下の分布の面積を求める関数を定義する。

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

cul_p <- function(x, min = NULL, max = NULL, continuous = F, length = 200){

  if(is.null(min)) min <- min(x)

  if(is.null(max)) max <- max(x)

  if(continuous == T){

    label <- seq(min, max, length = length)

  } else {

    label <- min:max

  }

  data.frame(label, p = sapply(label, function(tmp) sum(x <= tmp)/length(x)))

}

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

つまり、図3のヒストグラムのある値以下の分布面積を求める関数であり、累積分布関数を求めることと等しい。この値は確率分布の面積であるから、いわゆるp値と同じで、0から始まり、1まで増加する。各乱数の累積分布関数を表示してみよう。

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

pnb <- cul_p(rnb, 0, 150)

pnb_h <- cul_p(rnb_h, 0, 150)

pp <- cul_p(rp, 0, 150)

pg <- cul_p(rg, 0, 150)

plotn(p ~ label, data = pnb, xlab = "x")#図4の描画

plotn(p ~ label, data = pnb_h, xlab = "x")#図5の描画

plotn(p ~ label, data = pp, xlab = "x")#図6の描画

plotn(p ~ label, data = pg, xlab = "x")#図7の描画

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

これを分布の最小値から最大値まで連続的に計算した値は、同じ分布であれば等しく、y=xの直線状に乗るはずである。2つの分布のp値をx軸、y軸にしてプロットしたものをP-Pプロットと呼び、後述するQ-Qプロットとともに分布の近さを測るのに使われる。

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

ps <- full_join(pnb, pnb_h, by = "label")

colnames(ps)[2:3] <- c("nb", "nb_h")

plotn(nb ~ nb_h, data = ps, 

      xlab = "Mixed poisson and gamma", 

      ylab = "Negative binomial")#図8の描画

#負の二項分布と混合分布のP-Pプロット

overdraw(abline(a = 0, b = 1, col = "red"))

ps <- full_join(pnb, pp, by = "label")

colnames(ps)[2:3] <- c("nb", "p")

plotn(nb ~ p, data = ps, 

      xlab = "Poisson", 

      ylab = "Negative binomial")#図9の描画

#負の二項分布とポアソン分布のP-Pプロット

overdraw(abline(a = 0, b = 1, col = "red"))

ps <- full_join(pnb, pg, by = "label")

colnames(ps)[2:3] <- c("nb", "g")

plotn(nb ~ g, data = ps, 

      xlab = "Gamma", 

      ylab = "Negative binomial")#図10の描画

#負の二項分布とガンマ分布のP-Pプロット

overdraw(abline(a = 0, b = 1, col = "red"))

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

すると、負の二項分布と混合分布では1直線に乗り、区別がつかないことを示している。負の二項分布とガンマ分布は近い分布だが、わずかに異なることもわかる。

　もう一つ、分布の近さを測るためにQ-Qプロットを紹介しよう。上記ではp値を使ったが、そのp値を与える確率変数の値をq値と呼ぶ（分位数quantileに由来する）。定義から、累積分布関数の逆関数であり、分位点関数とも呼ばれる。

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

by <- 0.001

qnb <- data.frame(label = (0:(1/by))/10, q = quantile(rnb, probs = seq(0,1, by = by)))

qnb_h <- data.frame(label =(0:(1/by))/10, q = quantile(rnb_h, probs = seq(0,1, by = by)))

qp <- data.frame(label =(0:(1/by))/10, q = quantile(rp, probs = seq(0,1, by = by)))

qg <- data.frame(label =(0:(1/by))/10, q = quantile(rg, probs = seq(0,1, by = by)))

plotn(q ~ label, data = qnb, xlab = "quantile (%)")#図11の描画

plotn(q ~ label, data = qnb_h, xlab = "quantile (%)")#図12の描画

plotn(q ~ label, data = qp, xlab = "quantile (%)")#図13の描画

plotn(q ~ label, data = qg, xlab = "quantile (%)")#図14の描画

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

P-Pプロットと同様に、q値も同じ分布ならy=xの直線上に乗る。Q-Qプロットは各種の線形モデルで、残差の分布を確認するときに使われる。

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

qs <- full_join(qnb, qnb_h, by = "label")

colnames(qs)[2:3] <- c("nb", "nb_h")

plotn(nb ~ nb_h, data = qs, 

      xlab = "Mixed poisson and gamma", 

      ylab = "Negative binomial")#図15の描画

#負の二項分布とポアソン分布のQ-Qプロット

overdraw(abline(a = 0, b = 1, col = "red"))

qs <- full_join(qnb, qp, by = "label")

colnames(qs)[2:3] <- c("nb", "p")

plotn(nb ~ p, data = qs, 

      xlab = "Poisson", 

      ylab = "Negative binomial")#図16の描画

#負の二項分布とポアソン分布のQ-Qプロット

overdraw(abline(a = 0, b = 1, col = "red"))

qs <- full_join(qnb, qg, by = "label")

colnames(qs)[2:3] <- c("nb", "g")

plotn(nb ~ g, data = qs, 

      xlab = "Gamma", 

      ylab = "Negative binomial")#図17の描画

#負の二項分布とポアソン分布のQ-Qプロット

overdraw(abline(a = 0, b = 1, col = "red"))

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

すると、負の二項分布と混合分布では1直線に乗り、区別がつかないことを示している。負の二項分布とガンマ分布は近い分布だが、わずかに異なることもわかる。以上のように、負の二項分布はガンマ分布とポアソン分布の混合分布で表せる。

　実際に、下記のように混合分布の定義から計算することもできる。

負の二項分布を使ったGLM（負の二項回帰）

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

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

例えば、一つの説明変数xが存在する時、上記のような状況でデータが生成されていると考える。そう考えると、xiが得られた時のyiが生じる条件付き確率は上記のように表現できる。尤度は下記の通りになる。説明変数も組み込んで、整理しきってしまう。

対数尤度が説明変数xを含まなければ下記のようにβ0を解析的に解くことができる。

すると、exp(β0)の推定値はyの平均を満たすことがわかる。

　一方で、オッズ比(1-p)/pをdispersion parameterととらえる、つまり(1-p)/p = 定数とする場合は下記のとおりである。こちらのパラメータ化は途中で計算が極めて煩雑になるので図に示す程度にとどめる。

検出力

　ではRを使って、GLMの検出力に関して検討しよう。以下のように、線形予測子の切片の係数は0ではないが、説明変数xの係数は0である場合をまず考える。今回はサイズパラメータをdispersion parameterとする場合を取り上げる（最後の解析で、オッズ比をdispersion parameterとする場合も取り上げる）。リンク関数は対数リンクを想定する。平均が(1-pi)r/pi = exp(β・xi)に従うことから、pi = r/(r+exp(β・xi))となり、rはサイズパラメータであることを思い出せば、下記のように指定できる。データを生成した後で、描画してみよう。

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

b0 <- 3

b1 <- 0

size <- 6

n <- 30

x <- runif(n, 0, 30)

y <- rnbinom(n = n, size = size, prob = size/(size + exp(b0 + b1*x)))

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

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

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

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

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

fit2 <- summary(fit1)

fit2

##  Family: nbinom2  ( log )

## Formula:          y ~ x

## 

##      AIC      BIC   logLik deviance df.resid 

##    230.4    234.6   -112.2    224.4       27 

## 

## 

## Dispersion parameter for nbinom2 family (): 6.26 

## 

## Conditional model:

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

## (Intercept)  3.240902   0.156924  20.653   <2e-16 ***

## x           -0.004368   0.009035  -0.483    0.629    

## ---

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

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

yy <- exp(b0e + b1e * xx)

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

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

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

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

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

b0 <- 3

b1 <- 0

size <- 6

n <- 30

p1 <- NULL

p2 <- NULL

p3 <- NULL

p4 <- NULL

for (i in 1:1000){

  x <- runif(n, 0, 30)

  y <- rnbinom(n = n, size = size, prob = size/(size + exp(b0 + b1*x)))

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

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

  fit2 <- summary(fit1)

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

  fit4 <- summary(fit3)

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

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

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

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

}

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

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

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

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

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

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

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

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

## [1] 1000

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

## [1] 59

sum(p3 < 0.05)#ポアソン分布を仮定した時の切片の検出力

## [1] 1000

sum(p4 < 0.05)#ポアソン分布を仮定した時のx係数の危険率

## [1] 352

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

すると、負の二項分布を仮定した時は、xの係数の危険率は5%程度なのに対して、ポアソン分布を仮定すると35%程度まで危険率が高まった。ゆえに、負の二項分布を仮定することがより適切であろう。

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

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

b0 <- 3.8

b1 <- -0.05

size <- 6

n <- 30

x <- runif(n, 0, 30)

y <- rnbinom(n = n, size = size, prob = size/(size + exp(b0 + b1*x)))

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

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

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

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

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

fit2 <- summary(fit1)

fit2

##  Family: nbinom2  ( log )

## Formula:          y ~ x

## 

##      AIC      BIC   logLik deviance df.resid 

##    213.1    217.3   -103.6    207.1       27 

## 

## 

## Dispersion parameter for nbinom2 family (): 8.47 

## 

## Conditional model:

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

## (Intercept)  3.757021   0.143243  26.228  < 2e-16 ***

## x           -0.057973   0.009388  -6.175 6.62e-10 ***

## ---

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

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

yy <- exp(b0e + b1e * xx)

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

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

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

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

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

b0 <- 3.8

b1 <- -0.05

size <- 6

n <- 30

p1 <- NULL

p2 <- NULL

for (i in 1:1000){

  x <- runif(n, 0, 30)

  y <- rnbinom(n = n, size = size, prob = size/(size + exp(b0 + b1*x)))

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

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

  fit2 <- summary(fit1)  

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

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

}

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

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

sum(p1 < 0.05)

## [1] 1000

sum(p2 < 0.05)

## [1] 995

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

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

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

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

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

y <- c(1, 9, 3, 4, 6, 4, 3, 0, 3, 1, 3, 1, 3, 4, 4,

       4, 2, 1, 5, 7, 3, 2, 7, 6, 4, 2, 4, 3, 0, 3, 

       3, 2, 2, 0, 3, 3, 6, 3, 4, 0, 1, 2, 4, 3, 3, 

       7, 5, 5, 3, 0, 2, 0, 4, 1, 4, 5, 3, 8, 1, 2, 

       1, 11, 7, 5, 4, 1, 5, 2, 6, 4, 9, 3, 1, 4, 3, 

       8, 6, 8, 4, 0, 5, 5, 4, 4, 6, 2, 3, 4, 6, 5, 

       4, 5, 6, 3, 2, 2, 3, 1, 8, 1)

histn(y)#図28の描画

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

データは0より大きい離散値データである。ポアソン分布よりも分散は大きそうである。実際、データは平均よりも分散が大きい。

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

mean(y)

## [1] 3.62

var(y)

## [1] 5.207677

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

こういう時は負の二項分布を解析に使うとよい。一旦は練習なので、説明変数がないデータの解析をしてみることにする。まず、今回はsizeパラメータをdispersion parameterとする場合を考える。GLMで解析してみると下記のようになる。

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

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

fit2 <- summary(fit1)

fit2

##  Family: nbinom2  ( log )

## Formula:          y ~ 1

## 

##      AIC      BIC   logLik deviance df.resid 

##    442.2    447.4   -219.1    438.2       98 

## 

## 

## Dispersion parameter for nbinom2 family (): 7.89 

## 

## Conditional model:

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

## (Intercept)  1.28647    0.06348   20.27   <2e-16 ***

## ---

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

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

すると切片の係数は有意に0と異なっていることがわかる。今回の場合は負の二項分布の平均は下記のように3.62であり、確かにexp(切片の係数)と一致している。

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

mean(y)

## [1] 3.62

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

exp(b0e)

## [1] 3.62

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

さてここで尤度のことを思い出そう。今回推定された係数は尤度を最大にするように選ばれているはずである。対数尤度を計算するために下記のように関数を定義する。

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

f <- function(y, s, b0t){

  sapply(b0t, function(tmp) sum(log(dnbinom(x = y, size = s, 

                                            prob = s/(s + exp(tmp))

                                          )

                                    )

                                )

         )

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

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

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

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

size_e <- fit2$sigma

db0 <- 1

b0min <- b0e - db0

b0max <- b0e + db0

logL <- f(y, size_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")) #図29の図示

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

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

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

db0 <- 0.01

b0min <- b0e - db0

b0max <- b0e + db0

logL <- f(y, size_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")) #図30の図示

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

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

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

f(y, size_e, b0e)

## [1] -219.1098

-2*f(y, size_e, b0e)

## [1] 438.2196

-2*f(y, size_e, b0e) + 2*2

## [1] 442.2196

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

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

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

x <- c(22.6, 23.1, 26.4, 5.8, 20.8, 10.8, 14.1, 24.9, 26.6,

       11.6, 15.3, 1.9, 2.6, 13.6, 8.1, 27.7, 21.9, 10.7, 

       28.1, 4.8, 27.2, 1.2, 12.1, 19, 24.2, 22.3, 28.6, 

       21.4, 21, 20.4, 7.2, 27.6, 19.4, 17.5, 25.9, 6.9, 

       21.2, 27.6, 28.9, 22.7, 0.2, 15, 23.9, 6.3, 3.6, 

       1.8, 4.9, 10.5, 20, 21.9)

y <- c(12, 14, 5, 92, 20, 60, 28, 10, 14, 52, 26, 127, 93, 

       13, 91, 8, 14, 47, 13, 118, 12, 112, 53, 35, 5, 9, 

       8, 17, 31, 13, 80, 7, 13, 23, 8, 49, 12, 17, 2, 14,

       134, 57, 7, 120, 52, 141, 96, 62, 39, 12)

plotn(y ~ x)#図31の図示

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

このデータだけでポアソン回帰した時に過分散かどうかを判断するのは難しい。では解析をしてみよう。今回もsizeパラメータをdispersion parameterとする場合を考える。解析ののち、回帰曲線を描く。あわせてポアソン回帰して、ピアソンのχ^2と自由度の比を求めてみよう。

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

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

fit2 <- summary(fit1)

fit2

##  Family: nbinom2  ( log )

## Formula:          y ~ x

## 

##      AIC      BIC   logLik deviance df.resid 

##    379.1    384.8   -186.5    373.1       47 

## 

## 

## Dispersion parameter for nbinom2 family (): 10.8 

## 

## Conditional model:

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

## (Intercept)  5.051185   0.101263   49.88   <2e-16 ***

## x           -0.104092   0.005979  -17.41   <2e-16 ***

## ---

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

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

fit4 <- summary(fit3)

fit4

## 

## Call:

## glm(formula = y ~ x, family = poisson(link = "log"))

## 

## Deviance Residuals: 

##     Min       1Q   Median       3Q      Max  

## -5.6348  -1.4550  -0.3872   1.2170   4.2538  

## 

## Coefficients:

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

## (Intercept)  5.003511   0.033022  151.52   <2e-16 ***

## x           -0.100187   0.002766  -36.22   <2e-16 ***

## ---

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

## 

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

## 

##     Null deviance: 1760.20  on 49  degrees of freedom

## Residual deviance:  207.79  on 48  degrees of freedom

## AIC: 465.94

## 

## Number of Fisher Scoring iterations: 4

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

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

size_e <- fit2$sigma

E <- exp(b0e+b1e*x)

r <- 1/size_e

hr <- (y - E)/sqrt(E*(1+r*E)) #負の二項分布を仮定した時のpearson residuals

sum(hr^2)/(length(y)-2)

## [1] 1.069598

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

## [1] 1.069598

b0p <- coef(fit4)[1,1]

b1p <- coef(fit4)[2,1]

E <- exp(b0p+b1p*x)

hr <- (y - E)/sqrt(E) #ポアソン分布を仮定した時のpearson residuals

sum(hr^2)/(length(y)-2)

## [1] 4.236467

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

## [1] 4.236467

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

確かに、ポアソン分布を仮定した時、比が1より大きいのでポアソン回帰では過分散である。負の二項分布を仮定すれば、比は1に近いので、負の二項分布がより適切であることがわかる。

　では、回帰曲線を描いてみよう。

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

plotn(y ~ x)#図32の図示

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

yy <- exp(b0e + b1e * xx)

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

yy <- exp(b0p + b1p * xx)

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

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

係数自体は負の二項分布でも、ポアソン分布でも類似したものが推定される。

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

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

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

  mapply(function(tmp1, tmp2) {

    sum(log(dnbinom(x = y, size = s, prob = s/(s + exp(tmp1 + tmp2*x))

                    )

            )

    )

    }, 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, s = size_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)#図33の図示

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")#図34の図示

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

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

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

f(x, y, size_e, b0e, b1e)

## [1] -186.5314

-2*f(x, y, size_e, b0e, b1e)

## [1] 373.0629

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

## [1] 379.0629

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

　最後、次のようなデータである。図示も併せて行おう。

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

x <- c(1.4, 6.9, 0.7, 8.6, 9.4, 10.4, 13.1, 2.6, 

       20, 5, 23.6, 13.6, 22, 0.4, 22.1, 13.2, 27.1,

       9.9, 1.2, 11.8, 16.3, 24.5, 25.7, 20.7, 14.7,

       10.2, 1.1, 26.6, 24, 28.4, 10.9, 5.1, 26.6, 6, 

       20.7, 1.7, 23.8, 8.6, 12.8, 2.6, 7.7, 19.9, 1.8,

       9.4, 19.6, 23.5, 23.2, 19.8, 4.1, 17.1)

y <- c(180, 66, 112, 57, 57, 35, 24, 105, 34, 73, 

       8, 31, 28, 121, 11, 21, 11, 50, 159, 61, 57,

       6, 11, 15, 19, 68, 165, 9, 9, 8, 47, 80, 23,

       63, 5, 117, 15, 35, 24, 118, 37, 24, 149, 82, 

       26, 12, 15, 20, 85, 27)

plotn(y ~ x)#図35の図示

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

このデータだけでポアソン回帰した時に過分散かどうかを判断するのは難しい。では解析をしてみよう。今回はオッズ比をdispersion parameterとする場合を考える。その場合は、nbinom1関数をfamilyに指定するのであった。解析ののち、回帰曲線を描く。あわせてポアソン回帰して、ピアソンのχ^2と自由度の比を求めてみよう。

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

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

fit2 <- summary(fit1)

fit2

##  Family: nbinom1  ( log )

## Formula:          y ~ x

## 

##      AIC      BIC   logLik deviance df.resid 

##    399.8    405.6   -196.9    393.8       47 

## 

## 

## Dispersion parameter for nbinom1 family (): 3.41 

## 

## Conditional model:

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

## (Intercept)  4.987026   0.059833   83.35   <2e-16 ***

## x           -0.102008   0.006115  -16.68   <2e-16 ***

## ---

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

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

fit4 <- summary(fit3)

fit4

## 

## Call:

## glm(formula = y ~ x, family = poisson(link = "log"))

## 

## Deviance Residuals: 

##     Min       1Q   Median       3Q      Max  

## -3.9782  -1.6199  -0.3663   1.0015   4.9195  

## 

## Coefficients:

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

## (Intercept)  4.996239   0.028244  176.90   <2e-16 ***

## x           -0.103309   0.002889  -35.76   <2e-16 ***

## ---

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

## 

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

## 

##     Null deviance: 1847.26  on 49  degrees of freedom

## Residual deviance:  223.11  on 48  degrees of freedom

## AIC: 496.5

## 

## Number of Fisher Scoring iterations: 4

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

## [1] 1.09163

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

## [1] 4.847838

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

確かに、ポアソン分布を仮定した時、比が1より大きいのでポアソン回帰では過分散である。負の二項分布を仮定すれば、比は1に近いので、負の二項分布がより適切であることがわかる。

　では、回帰曲線を描いてみよう。

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

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

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

odds_e <- fit2$sigma

b0p <- coef(fit4)[1,1]

b1p <- coef(fit4)[2,1]

plotn(y ~ x)#図36の図示

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

yy <- exp(b0e + b1e * xx)

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

yy <- exp(b0p + b1p * xx)

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

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

係数自体は負の二項分布でも、ポアソン分布でも類似したものが推定される。

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

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

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

  mapply(function(tmp1, tmp2) {

    sum(log(dnbinom(x = y, size = exp(tmp1 + tmp2*x)/o, prob = 1/(1+o))

                    )

            )

    }, 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, o = odds_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)#図37の図示

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")#図38の図示

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

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

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

f(x, y, odds_e, b0e, b1e)

## [1] -196.9106

-2*f(x, y, odds_e, b0e, b1e)

## [1] 393.8213

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

## [1] 399.8213

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

以上のように、より大きな分散を持つ確率分布を使うことで、自然と過分散なデータに対応できた。次回も1つの確率分布を仮定して解析を行い、次々回から本格的なGLMMに挑戦する。