VGAM (version 1.1-11)

# dirmultinomial: Fitting a Dirichlet-Multinomial Distribution

## Description

Fits a Dirichlet-multinomial distribution to a matrix response.

## Usage

dirmultinomial(lphi = "logitlink", iphi = 0.10, parallel = FALSE,
zero = "M")

## Value

An object of class "vglmff" (see

vglmff-class). The object is used by modelling functions such as vglm, rrvglm

and vgam.

If the model is an intercept-only model then @misc (which is a list) has a component called shape which is a vector with the

$$M$$ values $$\pi_j(1/\phi-1)$$.

## Arguments

lphi

Link function applied to the $$\phi$$ parameter, which lies in the open unit interval $$(0,1)$$. See Links for more choices.

iphi

Numeric. Initial value for $$\phi$$. Must be in the open unit interval $$(0,1)$$. If a failure to converge occurs then try assigning this argument a different value.

parallel

A logical (formula not allowed here) indicating whether the probabilities $$\pi_1,\ldots,\pi_{M-1}$$ are to be equal via equal coefficients. Note $$\pi_M$$ will generally be different from the other probabilities. Setting parallel = TRUE will only work if you also set zero = NULL because of interference between these arguments (with respect to the intercept term).

zero

An integer-valued vector specifying which linear/additive predictors are modelled as intercepts only. The values must be from the set $$\{1,2,\ldots,M\}$$. If the character "M" then this means the numerical value $$M$$, which corresponds to linear/additive predictor associated with $$\phi$$. Setting zero = NULL means none of the values from the set $$\{1,2,\ldots,M\}$$. See CommonVGAMffArguments for more information.

Thomas W. Yee

## Warning

This VGAM family function is prone to numerical problems, especially when there are covariates.

## Details

The Dirichlet-multinomial distribution arises from a multinomial distribution where the probability parameters are not constant but are generated from a multivariate distribution called the Dirichlet distribution. The Dirichlet-multinomial distribution has probability function $$P(Y_1=y_1,\ldots,Y_M=y_M) = {N_{*} \choose {y_1,\ldots,y_M}} \frac{ \prod_{j=1}^{M} \prod_{r=1}^{y_{j}} (\pi_j (1-\phi) + (r-1)\phi)}{ \prod_{r=1}^{N_{*}} (1-\phi + (r-1)\phi)}$$ where $$\phi$$ is the over-dispersion parameter and $$N_{*} = y_1+\cdots+y_M$$. Here, $$a \choose b$$ means $$a$$ choose $$b$$'' and refers to combinations (see choose). The above formula applies to each row of the matrix response. In this VGAM family function the first $$M-1$$ linear/additive predictors correspond to the first $$M-1$$ probabilities via $$\eta_j = \log(P[Y=j]/ P[Y=M]) = \log(\pi_j/\pi_M)$$ where $$\eta_j$$ is the $$j$$th linear/additive predictor ($$\eta_M=0$$ by definition for $$P[Y=M]$$ but not for $$\phi$$) and $$j=1,\ldots,M-1$$. The $$M$$th linear/additive predictor corresponds to lphi applied to $$\phi$$.

Note that $$E(Y_j) = N_* \pi_j$$ but the probabilities (returned as the fitted values) $$\pi_j$$ are bundled together as a $$M$$-column matrix. The quantities $$N_*$$ are returned as the prior weights.

The beta-binomial distribution is a special case of the Dirichlet-multinomial distribution when $$M=2$$; see betabinomial. It is easy to show that the first shape parameter of the beta distribution is $$shape1=\pi(1/\phi-1)$$ and the second shape parameter is $$shape2=(1-\pi)(1/\phi-1)$$. Also, $$\phi=1/(1+shape1+shape2)$$, which is known as the intra-cluster correlation coefficient.

## References

Paul, S. R., Balasooriya, U. and Banerjee, T. (2005). Fisher information matrix of the Dirichlet-multinomial distribution. Biometrical Journal, 47, 230--236.

Tvedebrink, T. (2010). Overdispersion in allelic counts and $$\theta$$-correction in forensic genetics. Theoretical Population Biology, 78, 200--210.

Yu, P. and Shaw, C. A. (2014). An Efficient Algorithm for Accurate Computation of the Dirichlet-Multinomial Log-Likelihood Function. Bioinformatics, 30, 1547--54.

dirmul.old, betabinomial, betabinomialff, dirichlet, multinomial.

## Examples

Run this code
nn <- 5; M <- 4; set.seed(1)
ydata <- data.frame(round(matrix(runif(nn * M, max = 100), nn, M)))
colnames(ydata) <- paste("y", 1:M, sep = "")  # Integer counts

fit <- vglm(cbind(y1, y2, y3, y4) ~ 1, dirmultinomial,
data = ydata, trace = TRUE)
depvar(fit)  # Sample proportions
weights(fit, type = "prior", matrix = FALSE)  # Total counts per row

if (FALSE) {
ydata <- transform(ydata, x2 = runif(nn))
fit <- vglm(cbind(y1, y2, y3, y4) ~ x2, dirmultinomial,
data = ydata, trace = TRUE)
Coef(fit)
coef(fit, matrix = TRUE)
(sfit <- summary(fit))
vcov(sfit)
}


Run the code above in your browser using DataLab