rrisk.BayesZIP: Bayes estimation of a zero inflated Poisson (ZIP) model

In BfRstats/rriskBayes2: Predefined Bayes models fitted with Markov chain Monte Carlo (MCMC) (related to the 'rrisk' project)

Description

Zero-inflated Poisson data are count data with an excess number of zeros. The ZIP model involves the Poisson parameter lambda and the prevalence parameter pi.

Usage

rrisk.BayesZIP(data, prior.lambda = c(1, 10), prior.pi = c(0.8, 1), 
 simulation = FALSE, chains = 3, burn = 4000, 
 thin = 1, update = 5000, 
 max.time = "15minutes", plots = FALSE)

Arguments

data	Matrix, data frame or data set with positive integers, including zeros and of the minimal length 10.
prior.lambda	Numeric vector containing minimum and maximum of a uniform distribution used as prior for the Poisson parameter lambda, e.g., lambda ~ prior.lambda(*,*)=unif(*,*).
prior.pi	Numeric vector containing parameters of a beta distribution describing prior knowledge about prevalence (proportion of contaminated samples), e.g., pi ~ prior.pi(*,*)=beta(*,*).
simulation	Not used any longer.
chains	Positive single numeric value, number of independent MCMC chains (default 3).
burn	Positive single numeric value, length of the burn-in period (default 4000).
thin	Positive single numeric value (default 1). The samples from every k-th iteration will be used for inference, where k is the value of thin. Setting thin > 1 can help to reduce the autocorrelation in the sample.
update	Positive single numeric value, length of update iterations for estimation (default 5000).
max.time	Maximum time for which the function is allowed to extend the chains. Acceptable units include 'seconds', 'minutes', 'hours', 'days', 'weeks' (default "15minutes") (see autorun.jags).
plots	Logical, if TRUE the diagnostic plots will be displayed in separate windows .

Details

The ZIP model applies to count data and can be interpreted as a mixture distribution with one component comprising the 'true' zeros and another component of Poisson distributed values with density parameter lambda. The prevalence parameter pi refers to the proportion of the second, true non-zero component.

The Bayesian model for estimation prevalence and lambda parameter has in rjags/JAGS (originally BRugs/Winbugs) syntax following form

model{

   lambda ~ dunif(prior.lambda[1], prior.lambda[2])

   pi ~ dbeta(prior.pi[1], prior.pi[2]) 

   for (i in 1:n) {  
                   y[i]  ~ dpois(mu[i])

                   mu[i] <- I[i] * lambda  

                   I[i] ~ dbern(pi)  
                   }   
 }

Value

The function rrisk.BayesZIP returns an instance of the bayesmodelClass class containing following information:

convergence	Logical, whether the model has converged (assessed by the user).
results	Data frame containing statitsics of the posterior distribution.
jointpost	Data frame giving the joint posterior probability distribution.
nodes	Names of the parameters jointly estimated by the Bayes model.
model	Model in rjags/JAGS syntax as a character string.
chains	Number of independent MCMC chains.
burn	Length of burn-in period.
update	Length of update iterations for estimation.

Note

The convergence of the model should be checked using the diagnostic plots.

References

Bohning, D., Dietz, E., Schlattman, P., Mendonca, L. and Kirchner, U. (1999). The zero-inflated Poisson model and the decayed, missing and filled teeth index in dental epidemiology. Journal of the Royal Statistical Society, Series A 162:195-209.

Examples

#------------------------------------------
# Example of ZIP model
#------------------------------------------
# generate ZIP data
pi <- 0.01
n <- 200
lambda <- 3.5
zip.data <- rep(0,n)
zip.data[sample(1:n,n*pi,replace=FALSE)]<-rpois(n*pi,lambda=lambda)

# estimate using Bayes model for zero inflated data
resZIP <- rrisk.BayesZIP(data = zip.data, 
prior.lambda = c(0,100), 
prior.pi = c(1,1),
burn = 4000, 
max.time = '40seconds',
update = 4000)
resZIP

# compare with naive results ignoring ZIP model
pi.crude <- sum(zip.data>0)/n
lambda.crude <- mean(zip.data)
print(pi.crude)
print(lambda.crude)
resZIP@results

BfRstats/rriskBayes2 documentation built on May 5, 2019, 2:42 p.m.