Dirichlet-Multinomial Models for Joint Cell Counts

Why Model the Full Count Vector?

A donor-level cell-count table records how many cells belong to each of $K$ cell types. The entries are not separate outcomes: their sum is the donor’s total number of observed cells, and increasing one type’s share necessarily affects the shares of others. A model for one type versus all remaining types can describe a marginal proportion, but it does not by itself describe the joint vector or its dependence structure.

The Dirichlet-multinomial model extends multinomial sampling by allowing donors to vary in their underlying compositions. It is useful when donor-to-donor heterogeneity makes the observed counts more variable than a multinomial model permits. As in the earlier material on donors as independent units, biological donors—not individual cells—provide the independent replicates for condition comparisons.

Definition: For donor $d$, let $\mathbf{Y}_d=(Y_{d1},\ldots,Y_{dK})$ be the counts for $K$ cell types, and let $N_d=\sum_{k=1}^K Y_{dk}$. A Dirichlet-multinomial model with parameters $\boldsymbol{\alpha}=(\alpha_1,\ldots,\alpha_K)$, where every $\alpha_k>0$, is the mixture $$ \boldsymbol{\theta}_d\sim\operatorname{Dirichlet}(\alpha_1,\ldots,\alpha_K),\qquad \mathbf{Y}_d\mid\boldsymbol{\theta}_d\sim\operatorname{Multinomial}(N_d,\boldsymbol{\theta}_d). $$ Here $\boldsymbol{\theta}_d$ is the donor’s latent composition. Write $\phi=\sum_{k=1}^K\alpha_k$ for the concentration and $\pi_k=\alpha_k/\phi$ for the mean proportion of type $k$, so $\alpha_k=\phi\pi_k$, $\sum_k\pi_k=1$, and $\sum_kY_{dk}=N_d$.

The model’s probability mass function, conditional on the total $N_d$, is

P(\mathbf{Y}_d=\mathbf{y}_d\mid N_d,\boldsymbol{\pi},\phi) = \frac{N_d!}{\prod_{k=1}^K y_{dk}!} \frac{\Gamma(\phi)}{\Gamma(N_d+\phi)} \prod_{k=1}^K \frac{\Gamma(y_{dk}+\phi\pi_k)}{\Gamma(\phi\pi_k)}.

This is a model for the whole vector, including its fixed total. It does not treat the $K$ counts as independent. In a condition comparison, one can allow the mean composition $\boldsymbol{\pi}$ to vary by condition while assuming a shared concentration $\phi$, or fit different concentrations if the scientific question and available replication support doing so.

What the Concentration Parameter Controls

The mean composition $\boldsymbol{\pi}$ describes the average donor-level proportions. The concentration $\phi$ determines how tightly donors’ latent compositions cluster around that mean. Large $\phi$ concentrates $\boldsymbol{\theta}_d$ near $\boldsymbol{\pi}$; small $\phi$ allows greater donor-to-donor variation. In the limit as $\phi$ grows without bound, the Dirichlet-multinomial approaches the multinomial model with composition $\boldsymbol{\pi}$.

The model makes the overdispersion precise. For a single donor with total $N$, the count and proportion moments are

E(Y_k)=N\pi_k,\qquad \operatorname{Var}(Y_k) =N\pi_k(1-\pi_k)\frac{N+\phi}{1+\phi},

\operatorname{Cov}(Y_j,Y_k) =-N\pi_j\pi_k\frac{N+\phi}{1+\phi}\quad(j\ne k).

For proportions $Y_k/N$, the variance is $ \pi_k(1-\pi_k)(N+\phi)/(N(1+\phi)) $. The covariance is negative because the total is fixed: a donor with a higher-than-average share of one type must have less share available for at least one other type. This does not imply that all biological processes linking cell types are negative; it is a consequence of the model’s fixed-total composition.

Worked Example: Interpreting a Concentration Parameter

Suppose a donor has $N=10$ cells, and the mean composition is $\boldsymbol{\pi}=(0.5,0.3,0.2)$. Let $\phi=10$. The Dirichlet parameters are $\boldsymbol{\alpha}=(5,3,2)$. For the third cell type, the expected count is $10(0.2)=2$. Its count variance is

\operatorname{Var}(Y_3) =10(0.2)(0.8)\frac{10+10}{1+10} =1.6\frac{20}{11} \approx 2.909.

Under a multinomial model with the same $N$ and mean composition, the variance would be $10(0.2)(0.8)=1.6$. Thus the Dirichlet-multinomial assigns more variation to this count, reflecting variation in latent donor compositions as well as multinomial sampling.

The model also gives the probability of one complete count vector. For $\mathbf{y}=(4,3,3)$, which sums to 10, the probability is

\frac{10!}{4!3!3!} \frac{(5)_4(3)_3(2)_3}{(10)_{10}} = \frac{4200(1680)(60)(24)}{335221286400} \approx 0.03031,

where $(a)_r=a(a+1)\cdots(a+r-1)$ is the rising factorial. The calculation assigns probability to the joint vector, not to three independent counts.

Deriving the Overdispersion

The variance formula follows from the mixture representation. Conditional on $\boldsymbol{\theta}$, a multinomial count has mean $N\theta_k$ and variance $N\theta_k(1-\theta_k)$. Under a Dirichlet distribution with total concentration $\phi$, $E(\theta_k)=\pi_k$ and $\operatorname{Var}(\theta_k)=\pi_k(1-\pi_k)/(1+\phi)$. Applying the law of total variance gives

\begin{aligned} \operatorname{Var}(Y_k) &=E\{\operatorname{Var}(Y_k\mid\boldsymbol{\theta})\} +\operatorname{Var}\{E(Y_k\mid\boldsymbol{\theta})\}\\ &=N E\{\theta_k(1-\theta_k)\}+N^2\operatorname{Var}(\theta_k)\\ &=N\pi_k(1-\pi_k)\frac{\phi}{1+\phi} +N^2\frac{\pi_k(1-\pi_k)}{1+\phi}\\ &=N\pi_k(1-\pi_k)\frac{N+\phi}{1+\phi}. \end{aligned}

The first term in the second line is the average conditional multinomial variance; the second captures variation in the latent donor composition. Combining them explains why ignoring donor heterogeneity can understate uncertainty. The same conditioning argument, using the negative conditional multinomial covariance and the Dirichlet covariance, yields the negative covariance formula above.

Worked Example: Comparing Variability With the Multinomial

Consider $N=100$, $\boldsymbol{\pi}=(0.5,0.3,0.2)$, and $\phi=20$. For cell type 3, the variance of its observed proportion is

\operatorname{Var}(Y_3/100) =0.2(0.8)\frac{100+20}{100(1+20)} =0.16\frac{120}{2100} \approx 0.009143.

The corresponding standard deviation is $\sqrt{0.009143}\approx0.09562$, or about 9.56 percentage points. A multinomial model gives variance $0.2(0.8)/100=0.001600$ and standard deviation $0.04000$, or 4.00 percentage points. The larger Dirichlet-multinomial value reflects between-donor composition variability. It does not mean that the expected proportion changes: in both models it remains 0.2.

For types 1 and 3, the covariance of their proportions is

\operatorname{Cov}(Y_1/100,Y_3/100) =-0.5(0.2)\frac{120}{100(21)} \approx -0.005714.

This covariance is part of the joint uncertainty. Estimating or testing each type as if its proportion were unrelated to the other components discards that structure.

Testing for Condition-Associated Shifts

Suppose there are two conditions, with independent donors and the same set of $K$ cell types measured for every donor. A global null hypothesis of no mean-composition shift is

H_0:\boldsymbol{\pi}_1=\boldsymbol{\pi}_2 \qquad\text{versus}\qquad H_a:\boldsymbol{\pi}_1\ne\boldsymbol{\pi}_2.

Each composition sums to one, so it has $K-1$ free components. When the Dirichlet-multinomial model is correctly specified, the likelihood-ratio statistic compares the model with a common composition against one with condition-specific compositions. With independent donors, an interior parameter value, and enough donors for the large-sample approximation, twice the log-likelihood difference is approximately chi-square with $K-1$ degrees of freedom if the concentration is shared. The test concerns a global shift in mean composition, not a causal effect of condition.

A Wald calculation gives a transparent illustration of how overdispersion affects a global comparison. Under $H_0$, the covariance matrix for the first $K-1$ donor proportions follows from the Dirichlet-multinomial moments. The estimated difference between condition means is compared with its covariance, including the negative covariances among components. In practice, $\phi$ is usually estimated from the donor data; uncertainty in that estimate should be incorporated by the fitted likelihood or an appropriate covariance estimate.

Worked Example: A Global Two-Condition Comparison

Suppose $K=3$, with 20 independent donors per condition and $N=100$ cells per donor. Assume for this illustrative Wald calculation that the common concentration is known to be $\phi=20$. The observed mean proportions are $\hat{\boldsymbol{\pi}}_1=(0.60,0.25,0.15)$ and $\hat{\boldsymbol{\pi}}_2=(0.40,0.35,0.25)$. Under the null, the pooled composition is $(0.50,0.30,0.20)$. Use the first two components; the third is determined by their sum.

For the first two components, the Dirichlet-multinomial covariance per donor under the null is $fA$, where

f=\frac{N+\phi}{N(1+\phi)}=\frac{120}{2100}=0.0571429,\qquad A= \begin{pmatrix} 0.50(0.50)&-0.50(0.30)\\ -0.50(0.30)&0.30(0.70) \end{pmatrix} = \begin{pmatrix} 0.25&-0.15\\ -0.15&0.21 \end{pmatrix}.

The difference in the first two observed mean proportions is $\mathbf{d}=(0.20,-0.10)^\mathsf{T}$. Since each group has 20 donors, the covariance of the difference under the null is $V=(1/20+1/20)fA=0.00571429A$. The Wald statistic is

W=\mathbf{d}^{\mathsf T}V^{-1}\mathbf{d} =\frac{1}{0.00571429} \frac{0.21(0.20)^2+2(0.15)(0.20)(-0.10)+0.25(-0.10)^2}{0.25(0.21)-(-0.15)^2} \approx 28.58.

The denominator in the fraction is $0.0525-0.0225=0.0300$; its numerator is $0.0084-0.0060+0.0025=0.0049$. Thus the fraction is $0.1633$, and $W=0.1633/0.00571429\approx28.58$. Compared with a chi-square distribution with $2$ degrees of freedom, the approximate p-value is $P(\chi^2_2\ge28.58)=e^{-28.58/2}\approx6.2\times10^{-7}$. Under the stated model and approximation, this is strong evidence against equal mean compositions. It is not the probability that the null hypothesis is true.

This calculation treats $\phi$ as known and uses a large-sample reference distribution. A real analysis would generally estimate the concentration and assess whether the number of donors and model fit support the approximation. The count total per donor does not turn cells into independent replicates: the covariance above is for independent donor-level compositions.

Assumptions, Diagnostics, and Alternatives

The model is most defensible when the following conditions are plausible:

Each donor contributes one correctly aggregated count vector, and donors are independent within the comparison. Dependence from repeated or paired measurements requires a model that represents that design.
The cell-type categories are defined consistently across donors, and each vector has the stated total $N_d$. The model conditions on these totals; it does not explain why donors have different numbers of captured cells.
Within a condition, donor compositions can reasonably be represented by a Dirichlet distribution around the mean composition. A common concentration across conditions is an additional assumption when it is imposed.
The observed counts are adequately described by multinomial sampling conditional on each donor’s latent composition. Technical effects, unmodelled covariates, or sampling biases can violate this assumption.

Basic checks include comparing observed donor-to-donor proportion variability with the fitted variance, inspecting whether a few donors dominate the apparent shift, and checking whether fitted probabilities reproduce the prevalence of zero counts and broad patterns of the count vectors. Systematic misfit, strong subgroups, or extreme donors may signal that a single Dirichlet distribution is too restrictive. A zero count is allowed by the model; however, repeated zeros or a rare type with very little information can make its mean and concentration-related uncertainty unstable.

The multinomial model is simpler and is appropriate when the extra donor heterogeneity is negligible; otherwise, it can produce overly narrow uncertainty estimates. The beta-binomial model handles overdispersion for one proportion, such as one cell type versus the rest, but does not jointly model all cell types. The Dirichlet-multinomial provides a compact joint model, but its covariance structure is constrained by a single concentration and the mean composition. More flexible models may be appropriate when that structure is inadequate; the choice should be driven by the scientific design and observable model misfit, not by the smallest p-value.

Common Mistake and Key Takeaway

A common mistake is to pool all cells within each condition and then analyze the pooled proportions as if the number of cells were the number of independent experimental units. That reasoning treats cells from the same donor as independent evidence about the condition and ignores donor-to-donor composition variability. The correction is to use donor-level count vectors as replicates and represent their extra variability, for example with a Dirichlet-multinomial model. More sequenced cells can improve measurement of a donor’s composition, but they do not replace independent donors.

Key takeaway: The Dirichlet-multinomial models a donor’s complete cell-type count vector conditional on its total. Its mean composition describes average proportions, while its concentration controls extra donor heterogeneity. Condition tests must use independent donors and account for the joint compositional covariance.

Check Your Understanding

Answer each question, then compare with the solutions.

For $K=4$, how many independent components does a mean composition have, and what degrees of freedom are used for a global two-condition equality test under the shared-concentration likelihood-ratio approximation?
What happens to donor-to-donor heterogeneity as $\phi$ increases while $\boldsymbol{\pi}$ stays fixed?
For $N=50$, $\pi_k=0.4$, and $\phi=9$, calculate the Dirichlet-multinomial variance of the proportion $Y_k/N$.
Why is a zero count not automatically incompatible with a Dirichlet-multinomial model?
A study has many cells but only one donor in each condition. Can those cells supply independent donor-level replication for the condition test? Explain.

Solutions. (1) There are $K-1=3$ independent components, so the global test has $3$ degrees of freedom under the stated approximation. (2) The latent compositions become more concentrated around $\boldsymbol{\pi}$, so between-donor heterogeneity decreases and the model approaches multinomial sampling. (3) The variance is $0.4(0.6)(50+9)/(50(1+9))=0.24(59/500)=0.02832$, in squared proportion units. (4) The model assigns positive probability to vectors containing zero components; a zero is a valid count, though patterns of zeros can still reveal poor fit. (5) No. Cells from one donor do not provide independent donor-level replicates, so the condition comparison cannot estimate between-donor variation from those data alone.