dirmultinomial {VGAM}

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")

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.

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 jth linear/additive predictor (\eta_M=0 by definition for P[Y=M] but not for \phi) and j=1,\ldots,M-1. The Mth 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.

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

Warning

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

Note

The response can be a matrix of non-negative integers, or else a matrix of sample proportions and the total number of counts in each row specified using the weights argument. This dual input option is similar to multinomial.

To fit a ‘parallel’ model with the \phi parameter being an intercept-only you will need to use the constraints argument.

Currently, Fisher scoring is implemented. To compute the expected information matrix a for loop is used; this may be very slow when the counts are large. Additionally, convergence may be slower than usual due to round-off error when computing the expected information matrices.

Author(s)

Thomas W. Yee

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

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

## End(Not run)