# dmixnegbinom: The Negative Binomial Mixture Distribution In mederrRank: Bayesian Methods for Identifying the Most Harmful Medication Errors

## Description

Density function for a mixture of two negative binomial distributions.

## Usage

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

## 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).

## 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;α,β,E) = \frac{Γ(α + x)}{Γ(α) x!}\frac{1}{(1 + β/E)^x}\frac{1}{(1 + E/β)^α}

for x = 0,1,…,α, α, β, E > 0 and 0 < theta[5] ≤q 1. The mixture of two negative binomial distributions represents the marginal distribution of the counts N coming from Poisson data with parameter λ and a mixture of two gamma distributions as its prior. For details see the paper by Dumouchel (1999).

## Value

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

## Author(s)

Sergio Venturini [email protected],

Jessica A. Myers [email protected]

## 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.

dnbinom, rmixnegbinom.
  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 ## Not run: 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) ## End(Not run)