mederrRank (version 0.1.0)

dmixnegbinom: The Negative Binomial Mixture Distribution

Description

Density function for a mixture of two negative binomial distributions.

Usage

dmixnegbinom(x, theta, E, log.p = FALSE)

Value

dmixnegbinom gives the density corresponding to the E and theta values provided.

Arguments

x

vector of (non-negative integer) quantiles.

theta

vector of parameters for the negative binomial distribution mixture.

E

vector of (non-negative integer) expected counts.

log.p

logical; if TRUE, probabilities p are given as log(p).

Author

Sergio Venturini ,

Jessica A. Myers

Details

The mixture of two negative binomial distributions has density $$P(N = x) = theta[5] f(x;theta[1],theta[2],E) + (1 - theta[5]) f(x;theta[3],theta[4],E),$$ where $$f(x;\alpha,\beta,E) = \frac{\Gamma(\alpha + x)}{\Gamma(\alpha) x!}\frac{1}{(1 + \beta/E)^x}\frac{1}{(1 + E/\beta)^\alpha}$$ for \(x = 0,1,\ldots,\alpha\), \(\alpha, \beta, E > 0\) and \(0 < theta[5] \leq 1\). The mixture of two negative binomial distributions represents the marginal distribution of the counts \(N\) coming from Poisson data with parameter \(\lambda\) and a mixture of two gamma distributions as its prior. For details see the paper by Dumouchel (1999).

References

DuMouchel W. (1999), "Bayesian Data Mining in Large Frequency Tables, with an Application to the FDA Spontaneous Reporting System". The American Statistician, 53, 177-190.

Myers, J. A., Venturini, S., Dominici, F. and Morlock, L. (2011), "Random Effects Models for Identifying the Most Harmful Medication Errors in a Large, Voluntary Reporting Database". Technical Report.

See Also

dnbinom, rmixnegbinom.

Examples

if (FALSE) {
data("simdata", package = "mederrRank")
ni <- simdata@numi
theta0 <- c(10, 6, 100, 100, .1)
ans <- mixnegbinom.em(simdata, theta0, 50000, 0.01,
	se = FALSE, stratified = TRUE)
theta <- ans$theta.hat
N.E <- cbind(ans$N[1:ni], ans$E[1:ni])[sort(ans$N[1:ni], index.return = TRUE)$ix, ]
N.ix <- match(unique(N.E[, 1]), N.E[, 1])
N <- N.E[N.ix, 1]
E <- N.E[N.ix, 2]
dens <- dmixnegbinom(N, theta, E)
hist(N.E[, 1], breaks = 40, freq = FALSE)
points(N, dens)
}