負の二項分布のパラメータ推定：R と WinBUGSの比較

過分散した整数データのモデリングによく使う負の二項分布、Rではglm.nb関数でできるが、WinBUGSの例をほとんど見ない。ちょうど使う用事があったので使い方をテストしてみた。

# サンプルデータは、平均 exp(2)、dispersion parameter 1.5 の負の二項分布から抽出した N = 1000 のデータ（後で対数リンク関数を使用するのでexp()で表記）
require(MASS)
Y <- rnegbin(1000, mu=exp(2), theta=1.5) 
# 注：乱数なので実行する度に結果は少しずつ変わる

# まずはRのglm.nb関数による推定
require(MASS)
summary(glm.nb(Y ~ 1))

Coefficients:

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

(Intercept)  1.97519    0.02771   71.28   <2e-16 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(1.5895) family taken to be 1)

# 平均 1.97519（expの中身に相当）、dispersion parameter 1.5895 と順当な推定

# 次にWinBUGS、負の二項分布だからdnegbin、と試してみるが…（thetaがdispersion parameterに相当しそう）

data <- list(N=1000, Y=Y)

parm <- list(theta=rlnorm(1), alpha=rnorm(1))

model <- function() {

theta ~ dgamma(1E-1, 1E-2)

alpha ~ dnorm(0, 1E-6)

for (i in 1:N) {

Y[i] ~ dnegbin(p[i], theta) # 負の二項分布では平均をそのまま使えない

p[i] <- theta/(theta+mu[i]) # pと平均muの関係（pは0~1の値を取る）

log(mu[i]) <- alpha # muは正の値

}

}

mcmc <- wbugs(data, parm, model, n.iter=10000, n.burnin=5000, debug=T) # 自作ラッパー関数wbugsを使用

# order of negative binomial *** must be an integer というエラー

# theta は負の二項分布の狭義の定義通り、整数でないといけないようだ

# 戸惑いつつも多項分布の事前分布で thetaが正の整数しか取らないようにし再計算

data <- list(N=1000, Y=Y, q=rep(1,20)/20)

parm <- list(theta=sample(c(1:5), 1), alpha=rnorm(1))

model <- function() {

theta ~ dcat(q[]) # 事前分布1~20の整数（多項分布使用）

alpha ~ dnorm(0, 1E-6)

for (i in 1:N) {

Y[i] ~ dnegbin(p[i], theta)

p[i] <- theta/(theta+mu[i]) 

log(mu[i]) <- alpha 

}

}

mcmc <- wbugs(data, parm, model, n.iter=10000, n.burnin=5000)

               mean        sd          2.5%      25%        50%        75%        97.5%   Rhat   n.eff

theta        2.000       0.000    2.000    2.000       2.000       2.000       2.000          1        1

alpha       1.976       0.026    1.927    1.958       1.976       1.993       2.026          1  1500

deviance 6032.409 1.446 6031.000 6032.000 6032.000 6033.000 6036.000    1  1500

# 平均 alphaの推定は1.976と正しく推定できているのだが、thetaは当然のように整数値となった…

# Rの?rnbinomなどを読み返してみると、"dispersion parameter (the shape parameter of the gamma mixing distribution). Must be strictly positive, need not be integer" とある。

# gamma mixing distributionのshapeパラメータに相当するということか。

# なので、ポアソン分布にガンマ分布を混ぜた階層ベイズとして負の二項分布の推定を試みる

data <- list(N=1000, Y=Y)

parm <- list(shape=rlnorm(1,log(1),log(2)), alpha=rnorm(1), lambda=rlnorm(1000,log(1),log(2)),F) 

model <- function() {

shape ~ dgamma(1E-1, 1E-2)

alpha ~ dnorm(0, 1E-6)

for (i in 1:N) {

Y[i] ~ dpois(lambda[i]) # 負の二項分布を階層ベイズで表現

lambda[i] ~ dgamma(shape, rate[i])

rate[i] <- shape / mu[i] # rate と shape、平均 mu の関係を表す

log(mu[i]) <- alpha 

}

}

mcmc <- wbugs(data, parm, model, n.iter=10000, n.burnin=5000)

                    mean     sd     2.5%      25%      50%      75%    97.5%  Rhat   n.eff

shape           1.587  0.087    1.422    1.522    1.587    1.647    1.755 1.002  1200

alpha           1.975  0.028    1.917    1.957    1.975    1.993    2.029 1.003  1200

# lambdaは省略

# 平均 1.975、shape 1.587 と 元の値やglm.nbの推定とほぼ一致、これが正解のようだ

# それにしてもややこしい…負の二項分布には異なる定義が多すぎる！