VGAM (version 1.0-4)

dirmultinomial: Fitting a Dirichlet-Multinomial Distribution

Description

Fits a Dirichlet-multinomial distribution to a matrix response.

Usage

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

Arguments

lphi

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

iphi

Numeric. Initial value for $\varphi$. 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,\dots ,{\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,\dots ,M\}$. If the character "M" then this means the numerical value $M$, which corresponds to linear/additive predictor associated with $\varphi$. Setting zero = NULL means none of the values from the set $\{1,2,\dots ,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/\varphi -1)$.

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,\dots ,Y_M=y_M)=(\genfrac{}{}{0ex}{}{N_{\ast }}{y_1,\dots ,y_M})\frac{\prod _{j=1}^M\prod _{r=1}^{y_j}({\pi }_j(1-\varphi )+(r-1)\varphi )}{\prod _{r=1}^{N_{\ast }}(1-\varphi +(r-1)\varphi )}$$ where $\varphi$ is the over-dispersion parameter and $N_{\ast }=y_1+\cdots +y_M$. Here, $(\genfrac{}{}{0ex}{}{a}{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 $\varphi$) and $j=1,\dots ,M-1$. The $M$th linear/additive predictor corresponds to lphi applied to $\varphi$.

Note that $E(Y_j)=N_{\ast }{\pi }_j$ but the probabilities (returned as the fitted values) ${\pi }_j$ are bundled together as a $M$-column matrix. The quantities $N_{\ast }$ 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/\varphi -1)$ and the second shape parameter is $shape2=(1-\pi )(1/\varphi -1)$. Also, $\varphi =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.

See Also

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

Examples

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

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

# }
# NOT RUN {
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)
# }