Zero-inflated Poisson

ZIP (Zero-inflated Poisson) モデルのベイズ推定。Martin et al. (2005)を参考にしてみる。R2WinBUGSを使用。

Rのコード。n個のプロット中の割合pの部分に個体が分布し(分布範囲内だが個体数は0というプロットもあり)、その個体数は平均mのポアソン分布に従う。

set.seed(1234)

p <- 0.5
m <- 1.5

n <- 500
r <- rbinom(n, 1, p)
y <- r * rpois(n, m)

nz <- ifelse(y == 0, 0, 1)

library(R2WinBUGS)

winepath <- "/Applications/Darwine/Wine.bundle/Contents/bin/winepath"
wine <- "/Applications/Darwine/Wine.bundle/Contents/bin/wine"
init1 <- list(alpha = -10, beta = -10)
init2 <- list(alpha = 0, beta = 0)
init3 <- list(alpha = 10, beta = 10)
samples <- bugs(data = list(n = n, y = y, nz = nz),
                inits = list(init1, init2, init3),
                parameters = c("alpha", "beta"),
                model = "zip.bug",
                n.chains = 3, n.iter = 5000,
                useWINE = TRUE, WINEPATH = winepath, WINE = wine,
                bugs.directory = "c:/Program Files/WinBUGS14/",
                working.directory = getwd(),
                debug = FALSE)
samples.list <- mcmc.list(lapply(1:(samples$n.chain),
                                 function(chain)
                                 	as.mcmc(samples$sims.array[,chain,])))

BUGSのコード

var
  n,		# number of data
  y[n],		# data
  nz[n],	# is non-zero ?
  mu[n],	# 
  p, m,
  alpha, beta

model {
  for (i in 1:n) {
    zeros[i] <- 0;
    zeros[i] ~ dpois(mu[i]);
    mu[i] <- (1 - nz[i]) * -log((1 - p + p * exp(-m))) +
             nz[i] * -(log(p) - m + y[i] * log(m) - logfact(y[i]));
  }

  logit(p) <- alpha;
  log(m) <- beta;
  alpha ~ dnorm(0.0, 1.0E-4);
  beta ~ dnorm(0.0, 1.0E-4);
}

結果

> summary(samples.list)

Iterations = 1:358
Thinning interval = 1 
Number of chains = 3 
Sample size per chain = 358 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

               Mean      SD Naive SE Time-series SE
alpha      -0.01431 0.12576 0.003837       0.004535
beta        0.46924 0.06474 0.001975       0.002476
deviance 1203.11825 2.10666 0.064282       0.079939

2. Quantiles for each variable:

              2.5%       25%        50%       75%     97.5%
alpha      -0.2490   -0.1001   -0.01358 7.025e-02    0.2316
beta        0.3333    0.4270    0.47165 5.135e-01    0.5883
deviance 1201.0000 1202.0000 1203.00000 1.204e+03 1209.0000

> exp(-0.01431)/(1+exp(-0.01431))
[1] 0.4964226
> exp(0.46924)
[1] 1.598779

BUGSのコードはMartin (2005)のAppendix 1を参考にしたのだが、なかなかトリッキーだ。

zeros[i] <- 0;
zeros[i] ~ dpois(mu[i]);

このあたり。mu[i]が対数尤度に-1を掛けたものになるので、これを最小化している。

タグ：winbugs zero-inflated poisson MCMC