rrisk.BayesZINB: Bayes estimation of a zero inflated negative binomial (ZINB)...

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

Description

Zero-inflated Negative-Binomial data are count data with an excess number of zeros. The ZINB (also referred to as 'gamma-poisson') model involves the prevalence parameter pi. The negative binomial distribution can be seen as a poisson(λ) distribution, where λ follows a gamma distribution.

Usage

rrisk.BayesZINB(data, 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.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 ZINB model applies to count data and can be interpreted as a mixture distribution with one component comprising the 'true' zeros and another component of negative-binomially distributed values with density parameter λ following a gamma distribution with shape and mean parameters modelled as dgamma(shape = 0.01, lambda = 0.01). 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{

   pi  ~ dbeta(prior.pi[1], prior.pi[2])
    
   dam ~ dgamma(0.01, 0.01)
   
   db ~ dgamma(0.01, 0.01)

   for (i in 1:n) {
           y[i] ~ dpois(mu[i])
                   
           mu[i] <- I[i] * lambda[i]
                   
           I[i] ~ dbern(pi)
                   
           lambda[i] ~ dgamma(dam, db)
                   }
 }

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

Lunn, D. et al. (2012). The BUGS book: A practical introduction to Bayesian analysis. CRC press. p.353, 117

Examples

#------------------------------------------
#Example of a ZINB model (compare with rrisk.BayesZIP)
#------------------------------------------
#generate ZINB data
pi <- 0.1
n <- 200
zinb.data <- rep(0,n)
zinb.data[sample(1:n, n*pi, replace = FALSE)] <- rpois(n*pi, lambda = 4)

# estimate using Bayes model for zero inflated data
resZINB <- rrisk.BayesZINB(data = zinb.data, prior.pi = c(1,1),
 burn = 100, update = 4000, max.time = '40seconds')
resZINB

#------------------------------------------
#Example of a ZINB model 
#------------------------------------------
set.seed(42)
n_true_neg <- 60
n_true_pos <- 33
n <- n_true_pos + n_true_neg

prev_true <- n_true_pos / n

a <- 6
b <- 2
lambda_true <- rgamma(n_true_pos, a, b)

y_neg <- rep(0, n_true_neg)
y_pos <- rpois(n_true_pos, lambda_true)
y <- c(y_pos, y_neg)

pi_prior     <- c(1, 1)

resZINB <- rrisk.BayesZINB(data = y,
                            prior.pi = pi_prior,
                            simulation = FALSE,
                            chains = 3,
                            burn = 4000,
                            thin = 1,
                            max.time = '60seconds',
                            update = 10000,
                            plots = TRUE
                            )

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