######  6.1.2 Clarification of Compositional Heterogeneity

The intended notion of “compositional heterogeneity” herein does not include
irregularities corresponding to e.g. strains preferentially grouping with or
distancing from each other.

First, note that the Dirichlet-Multinomial can be thought of as a hierarchical
distribution, where the $f^{(s)}$ parameters of the Multinomial distribution
(living on the unit simplex due to the $\sum_{s=1}^{S}f^{(s)}=1$ constraint)
are Dirichlet distributed. (Whence the name “Dirichlet-Multinomial”.) Second,
while the marginal distributions of the Dirichlet distribution are anti-
correlated with one another, this is in some sense due only to the
normalization condition, $\sum_{s=1}^{S}f^{(s)}=1$. It can be shown that the
joint distribution of mutually independent Gamma random variables, when
normalized by their sum, is Dirichlet. In particular, the Dirichlet
distribution can not model a situation where e.g. particular strains $s_{1}$
and $s_{2}$ tend to attach to one another within the sampling pool.

One way to get distributions on the unit simplex able to model such
preferential grouping or distancing of strains is to use correlated random
variables with distributions derived from Gamma distributions and again
normalize by their sum. Such an approach would most closely resemble the
intuitive idea of a “correlated Dirichlet distribution”. Another approach is
to use the logistic-normal distribution. See [BGF01] for a study applying
methodology based on the latter approach to the study of multivariate count
data from ecology.

######  6.1.3 Over-Disperson of Dirichlet-Multinomial with respect to
Multinomial

The over-dispersion of the Dirichlet-Multinomial distribution with
concentration $\zeta_{C}S$, total number $n$, and frequencies $\vec{f}$ with
respect to the corresponding Multinomial distribution with total number $n$
and frequencies $\vec{f}$ is

$\frac{n-1}{1+\zeta_{C}S}\,.$ (6.2)

Thus the concentration is approximately proportional to the reciprocal of the
over-dispersion. (They are asymptotically equivalent up to a constant, i.e.
the limit of their ratio as ${\zeta_{C}\to\infty}$ is the constant $n-1$.)

######  6.1.4 hPoDM Definition

Under the hPoDM distribution, for the $d$’th droplet: the probability of the
strain distribution vector given the number of cells

$\mathbb{P}\left(\vec{\mathbf{N}}_{d}(0)=\vec{\mathbf{n}}\nonscript\;\middle|\nonscript\;N_{d}(0)=n\right)$

is determined by a Dirichlet-Multinomial distribution, while the probability
of the number of cells ${\mathbb{P}\left(N_{d}(0)=n\right)}$ is determined by
a Poisson distribution.

Thus for the unconditional probability:

$\begin{split}&\mathbb{P}\left(\vec{\mathbf{N}}_{d}(0)=\vec{\mathbf{n}}\right)\\\
=&\quad\mathbb{P}\left(N_{d}(0)=n\right)\\!\cdot\\!\mathbb{P}\left(\vec{\mathbf{N}}_{d}(0)=\vec{\mathbf{n}}\nonscript\;\middle|\nonscript\;N_{d}(0)=n\right)\\\
=&\quad\frac{e^{-\lambda}\lambda^{n}}{n!}\\!\cdot\\!\frac{\Gamma(\zeta_{C}S)}{\Gamma(\zeta_{C}S+n)}\\!\cdot\\!\left[\prod_{s=1}^{S}\frac{\Gamma(\zeta_{C}Sf^{(s)}+n^{(s)})}{\Gamma\left(\zeta_{C}Sf^{(s)}\right)}\right]\\!\cdot\\!\binom{n}{n^{(1)}\\!\cdots\\!n^{(S)}}\,.\end{split}$
(6.3)

#####  6.2 hPoDM Generalizes hPoMu

Just like for hTPMH, hPoDM converges to the default hPoMu model under
appropriate limits. Therefore when using the hPoDM family of distributions to
model the distribution of the initial formation of droplets, we “retain the
same language” used by the hPoMu distribution for that task.

In sections 6.2.1 and 6.2.2 I give the context needed and then show that the
hNBDM family really is an extension of hPoMu. Because hPoDM and hPoMu have the
same count distribution, this reduces to showing how the Dirichlet-Multinomial
distribution extends the multinomial distribution.

######  6.2.1 Approximate Form of Likelihood Ratio of DM with respect to
Multinomial and of hPoDM with respect to hPoMu and Bounds

Using Lemma 3.4 we get that

$\begin{split}&\frac{\Gamma(\zeta_{C}S)}{\Gamma(\zeta_{C}S+n)}\\!\cdot\\!\left[\prod_{s=1}^{S}\frac{\Gamma(\zeta_{C}Sf^{(s)}+n^{(s)})}{\Gamma(\zeta_{C}Sf^{(s)})}\right]\\\
=&\exp\left(\mathscr{R}_{\mathrm{DM}}(\vec{\mathbf{n}},\zeta_{C}S\vec{f})\right)\\!\cdot\\!\left(\frac{\zeta_{C}S}{\zeta_{C}S+n}\right)^{\zeta_{C}S+n-\frac{1}{2}}\\\
&\\!\cdot\\!\left[\prod_{s=1}^{S}\left(\frac{\zeta_{C}Sf^{(s)}+n^{(s)}}{\zeta_{C}Sf^{(s)}}\right)^{\zeta_{C}Sf^{(s)}+n^{(s)}-\frac{1}{2}}\right]\\!\cdot\\!\prod_{s=1}^{S}(f^{(s)})^{n^{(s)}}\,,\end{split}$
(6.4)

where (with the upper $\overline{\mathscr{R}_{\mathrm{}}}$ and lower
$\underline{\mathscr{R}_{\mathrm{}}}$ functions as defined before)

$\underline{\mathscr{R}_{\mathrm{}}}(\zeta_{C}S\vec{f},\zeta_{C}S\vec{f}+\vec{\mathbf{n}})<\mathscr{R}_{\mathrm{DM}}(\vec{\mathbf{n}},\zeta_{C}S\vec{f})<\overline{\mathscr{R}_{\mathrm{}}}(\zeta_{C}S\vec{f},\zeta_{C}S\vec{f}+\vec{\mathbf{n}})\,.$
(6.5)

Unfortunately the bounds in (6.5) appear to be tight in general only for
larger values of $\zeta_{C}$. (They become tighter as $\zeta_{C}\to\infty$;
cf. section 6.2.2.)

Thus the likelihood ratio of the Dirichlet-Multinomial distribution with
respect to the Multinomial distribution, which is the same as the likelihood
ratio of the hPoDM distribution with respect to the hPoMu distribution, based
on equation 6.4 above equals

$\exp\left(\mathscr{R}_{\mathrm{DM}}(\vec{\mathbf{n}},\zeta_{C}S\vec{f})\right)\\!\cdot\\!\left(\frac{\zeta_{C}S}{\zeta_{C}S+n}\right)^{\zeta_{C}S+n-\frac{1}{2}}\\!\cdot\\!\left[\prod_{s=1}^{S}\left(\frac{\zeta_{C}Sf^{(s)}+n^{(s)}}{\zeta_{C}Sf^{(s)}}\right)^{\zeta_{C}Sf^{(s)}+n^{(s)}-\frac{1}{2}}\right]\,.$
(6.6)

######  6.2.2 Sketch of Proof that hPoDM Converges in Distribution to hPoMu

Again, it suffices to show that the probability mass function of hPoDM
converges to that of hPoMu. Cf. section 5.2.3.

That the expression in (6.6) approaches $1$ as $\zeta_{C}\to\infty$ again
follows from the Squeeze Theorem and Lemma 3.3. (As before, that the
expressions involving the upper and lower Robbins functions approach $0$ is an
elementary consequence of the Archimedean property of the real numbers.) Thus
the hPoDM distribution can be seen to converge to the hPoMu distribution as
$\zeta_{C}\to\infty$.

Because the likelihood ratio of the Dirichlet-Multinomial distribution with
respect to the Multinomial distribution is the same as the likelihood ratio of
the hPoDM distribution with respect to the hPoMu distribution, it follows that
the above also shows how the Dirichlet-Multinomial distribution converges to
the Multinomial distribution as $\zeta_{C}\to\infty$.

Instead of showing that the likelihood ratio converges to $1$, one can also
show (e.g. according to one’s personal preference) that the logarithm of the
likelihood ratio converges to $0$, using Lemma 3.5 (the “logarithmic version”
of Lemma 3.3). The “fudge function” $\Phi$ is defined for convenience in
equation (3.A.3). Then, starting from (6.6), tedious algebra shows that the
logarithm of the likelihood ratio of (hPo)DM with respect to (hPo)Mu is

$\sum_{s=1}^{S}\Phi(\zeta_{C}Sf^{(s)},n^{(s)},n^{(s)}-\frac{1}{2})-\Phi(\zeta_{C}S,n,n-\frac{1}{2})+\mathscr{R}_{\mathrm{DM}}(\vec{\mathbf{n}},\zeta_{C}S\vec{f})\,.$
(6.7)

####  7 Working Model for Density Heterogeneity

The hPoDM working model has the same count distribution as hPoMu, and hence is
unable to model over-dispersion of the count distribution caused by density
heterogeneity. Thus the main change we need to make is from using the Poisson
distribution as the count distribution to using a distribution over-dispersed
relative to the Poisson distribution. Section 7.1 introduces a working model
that fills this gap, and then section 7.2 shows that this working model can
also be considered an extension of hPoMu.

#####  7.1 Introduction to hNBDM

As suggested by its name, the “hierarchical Negative Binomial Dirichlet-
Multinomial” (abbreviated as hNBDM) working model fills this gap by using a
new distribution for the count distribution, the negative binomial. Section
7.1.1 introduces the negative binomial distribution. Section 7.1.2 explains
how the negative binomial distribution models over-dispersion relative to the
Poisson distribution. Section 7.1.3 defines the hNBDM working model.

######  7.1.1 Negative Binomial (NB) Distribution

One may write the PMF of the negative binomial distribution as

$\frac{\Gamma(\zeta_{D}S+n)}{\Gamma(\zeta_{D}S)}\\!\cdot\\!\frac{(\zeta_{D}S)^{\zeta_{D}S}}{(\zeta_{D}S+\lambda)^{\zeta_{D}S+n}}\\!\cdot\\!\frac{\lambda^{n}}{n!}$
(7.1)

where $\lambda$ denotes the expected value of the distribution.

More typically this is parameterized in terms of $\zeta_{D}S$ instead of
$\zeta_{D}$. (Indeed in principle it should not even be necessary to specify a
number of strains $S$ in order to define this distribution.) I parameterize in
terms of $\zeta_{D}$ to be compatible with the parameterization of the
Dirichlet-Multinomial I use in section 6.1.1. See section 7.1.3 for how this
pays off later by simplifying the math.

The parameter $\zeta_{D}$ is herein called the “density concentration”. Using
the negative binomial distribution for modelling density heterogeneity has
precedent in the literature. For macroecological studies, the parameter
$\lambda$ “has been defined as the density of organisms in the area of
interest” [WFY84]. Herein $\lambda$ will be interpreted as the average cell
density in the sampling pool. Previous work has also claimed that $\zeta_{D}$
should be able to capture variability in the density of organisms: “The
practical use of [$\zeta_{D}S$] … in ecological studies of aggregation
requires caution because [$\zeta_{D}S$] is usually density-dependent” [CP89].

If a Poisson distribution with rate $\lambda$ corresponds to a homogeneous
cell density $\lambda$ throughout the sampling pool, a negative binomial
distribution with rate $\lambda$ is interpreted to have the same average cell
density. Larger values of density concentration $\zeta_{D}$ correspond to
smaller heterogeneity around the average density, and thus more closely
resemble the Poisson case (which is the limit as $\zeta_{D}\to\infty$), cf.
figure 4(a). Similarly, smaller values of density concentration $\zeta_{D}$
correspond to larger heterogeneity around the average density, cf. figure
4(b).

######  7.1.2 Over-Disperson of Negative Binomial with respect to Poisson

The over-dispersion of the negative binomial distribution with concentration
$\zeta_{D}S$ and mean $\lambda$ with respect to the Poisson distribution with
mean $\lambda$ is

$\frac{\lambda}{\zeta_{D}S}\,.$ (7.2)

Thus the concentration is proportional to the reciprocal of the over-
dispersion.

######  7.1.3 hNBDM Definition

For hNBDM we set $\zeta_{D}=\zeta_{C}=:\zeta$. So for the $d$’th droplet, the
strain distribution vector
$\mathbb{P}\left(\vec{\mathbf{N}}_{d}(0)=\vec{\mathbf{n}}\nonscript\;\middle|\nonscript\;N_{d}(0)=n\right)$
given the number of cells is Dirichlet-Multinomial distributed with parameter
$\zeta$, and the number of cells ${\mathbb{P}\left(N_{d}(0)=n\right)}$ is
Negative Binomial distributed with parameter $\zeta$. Thus the PMF defining
the unconditional probability is

$\begin{split}&\mathbb{P}\left(\vec{\mathbf{N}}_{d}(0)=\vec{\mathbf{n}}\right)\\\
=&\quad\mathbb{P}\left(N_{d}(0)=n\right)\\!\cdot\\!\mathbb{P}\left(\vec{\mathbf{N}}_{d}(0)=\vec{\mathbf{n}}\nonscript\;\middle|\nonscript\;N_{d}(0)=n\right)\\\
&\frac{\Gamma(\zeta S+n)}{\Gamma(\zeta S)}\\!\cdot\\!\frac{(\zeta S)^{\zeta
S}}{(\zeta S+\lambda)^{\zeta
S+n}}\\!\cdot\\!\frac{\lambda^{n}}{n!}\frac{\Gamma(\zeta S)}{\Gamma(\zeta
S+n)}\\!\cdot\\!\left[\prod_{s=1}^{S}\frac{\Gamma(\zeta
Sf^{(s)}+n^{(s)})}{\Gamma\left(\zeta
Sf^{(s)}\right)}\right]\\!\cdot\\!\binom{n}{n^{(1)}\\!\cdots\\!n^{(S)}}\end{split}$
(7.3)

which via tedious algebraic manipulations can be seen to equal

$\prod_{s=1}^{S}\left[\frac{\Gamma(\zeta Sf^{(s)}+n^{(s)})}{\Gamma(\zeta
Sf^{(s)})}\\!\cdot\\!\frac{(\zeta Sf^{(s)})^{\zeta Sf^{(s)}}}{(\zeta
Sf^{(s)}+f^{(s)}\lambda)^{\zeta
Sf^{(s)}+n^{(s)}}}\\!\cdot\\!\frac{(f^{(s)}\lambda)^{n^{(s)}}}{(n^{(s)})!}\right]\,.$
(7.4)

The marginal distributions of the hNBDM distribution correspond to mutually
independent Negative Binomial random variables whose parameters have been
rescaled by the $f^{(s)}$.

This makes the relationship of the Negative Binomial distribution to the
Dirichlet-Multinomial a direct analogue of the relationship of the Poisson
distribution to the Multinomial Distribution. Unsurprisingly these
observations have precedent in the literature. Cf. for example Theorem 1 of
[Zho18].

#####  7.2 hNBDM Generalizes hPoMu

Just like for hTPMH and hPoDM, hNBDM converges to the default hPoMu model
under appropriate limits. Therefore when using the hNBDM family of
distributions to model the distribution of the initial formation of droplets,
we “retain the same language” used by the hPoMu distribution for that task.

In sections 7.2.1 and 7.2.2 I give the context needed and then show that, as a
count distribution, the negative binomial distribution truly is a
generalization of the Poisson. In sections 7.2.3 and 7.2.4 I give the context
needed and then show that the hNBDM family really is an extension of hPoMu.

######  7.2.1 Approximate Form of Likelihood Ratio of Negative Binomial with
respect to Poisson and Bounds

Applying Lemma 3.4, one has

$\frac{\Gamma(\zeta_{D}S+n)}{\Gamma(\zeta_{D}S)}=e^{-n}\\!\cdot\\!\left(\frac{\zeta_{D}S+n}{\zeta_{D}S}\right)^{\zeta_{D}S+n-\frac{1}{2}}\\!\cdot\\!(\zeta_{D}S)^{n}\\!\cdot\\!\exp\left(\mathscr{R}_{\mathrm{NB}}(\zeta_{D}S,n)\right)\,,$
(7.5)

where $\mathscr{R}_{\mathrm{NB}}(\zeta_{D}S,n)$ satisfies

$\frac{1}{12(\zeta_{D}S+n)+1}-\frac{1}{12\zeta_{D}S}<\mathscr{R}_{\mathrm{NB}}(\zeta_{D}S,n)<\frac{1}{12(\zeta_{D}S+n)}-\frac{1}{12\zeta_{D}S+1}\,.$
(7.6)

It follows from equations (7.1) and (7.5) that the likelihood ratio of the
negative binomial distribution (with mean $\lambda$) with respect to the
Poisson distribution (with mean $\lambda$) is

$e^{-n}\\!\cdot\\!\left(\frac{\zeta_{D}S+n}{\zeta_{D}S}\right)^{\zeta_{D}S+n-\frac{1}{2}}\\!\cdot\\!\left(\frac{\zeta_{D}S}{\zeta_{D}S+\lambda}\right)^{\zeta_{D}S+n}\\!\cdot\\!e^{\lambda}\\!\cdot\\!\exp\left(\mathscr{R}_{\mathrm{NB}}(\zeta_{D}S,n)\right)\,.$
(7.7)

######  7.2.2 Sketch of Proof that NB Converges in Distribution to Poisson

It follows nearly immediately from Lemma 3.3 that the above expression (7.7)
converges to $1$ as $\zeta_{D}\to\infty$. (Again, both of the bounds for
$\mathscr{R}_{\mathrm{NB}}(\zeta_{D}S,n)$ from (7.6) approach $0$ for trivial
reasons, allowing us to apply the Squeeze Theorem.) In particular, it follows
that the negative binomial distribution converges to the Poisson distribution
as $\zeta_{D}\to\infty$.

######  7.2.3 Approximate Form of Likelihood Ratio of hNBDM with respect to
hPoMu and Bounds

Using either of the expressions (7.3) or (7.4) for the PMF of the hNBDM
distribution given above, more tedious algebraic manipulations (including
application of Lemma 3.4) allow one to show that the likelihood ratio of the
hNBDM distribution relative to the corresponding hPoMu distribution is

$\displaystyle\begin{split}&\left[\prod_{s=1}^{S}e^{-n^{(s)}}\\!\cdot\\!\left(\frac{\zeta
Sf^{(s)}+n^{(s)}}{\zeta Sf^{(s)}}\right)^{\zeta
Sf^{(s)}+n^{(s)}-\frac{1}{2}}\\!\cdot\\!\left(\frac{\zeta Sf^{(s)}}{\zeta
Sf^{(s)}+f^{(s)}\lambda}\right)^{\zeta
Sf^{(s)}+n^{(s)}}\\!\cdot\\!e^{f^{(s)}\lambda}\right]\\!\cdot\\!\mathscr{R}_{\mathrm{hNBDM}}(\zeta
S\vec{f},\vec{\mathbf{n}})\,,\end{split}$ (7.8)

where

$\mathscr{R}_{\mathrm{hNBDM}}(\zeta
S\vec{f},\vec{\mathbf{n}})=\mathscr{R}_{\mathrm{NB}}(\zeta
S,n)+\mathscr{R}_{\mathrm{DM}}(\zeta
S\vec{f},\vec{\mathbf{n}})=\sum_{s=1}^{S}\mathscr{R}_{\mathrm{NB}}(\zeta
Sf^{(s)},n^{(s)})\,.$ (7.9)

Observe that adding the bounds from (6.5) and (7.6) gives the following bounds
for $\mathscr{R}_{\mathrm{hNBDM}}(\zeta S\vec{f},\vec{\mathbf{n}})$:

$\begin{split}&\scalebox{0.75}{
$\displaystyle\sum_{s=1}^{S}\left[\frac{1}{12(\zeta
Sf^{(s)}\\!+\\!n^{(s)})\\!+\\!1}\\!-\\!\frac{1}{12\zeta
Sf^{(s)}}\right]\\!+\\!\left(\frac{1}{12\zeta
S\\!+\\!1}\\!-\\!\frac{1}{12\zeta S}\right)\\!+\\!\left(\frac{1}{12(\zeta
S\\!+\\!n)\\!+\\!1}\\!-\\!\frac{1}{12(\zeta S\\!+\\!n)}\right)$ }\\\
=&\quad\underline{\mathscr{R}_{\mathrm{}}}(\zeta S\vec{f},\zeta
S\vec{f}+\vec{\mathbf{n}})+\frac{1}{12(\zeta S+n)+1}-\frac{1}{12\zeta S}\\\
<&\quad\mathscr{R}_{\mathrm{DM}}(\zeta
S\vec{f},\vec{\mathbf{n}})+\mathscr{R}_{\mathrm{NB}}(\zeta
S,n)=\mathscr{R}_{\mathrm{hNBDM}}(\zeta S\vec{f},\vec{\mathbf{n}})\\\
<&\quad\overline{\mathscr{R}_{\mathrm{}}}(\zeta S\vec{f},\zeta
S\vec{f}+\vec{\mathbf{n}})+\frac{1}{12(\zeta S+n)}-\frac{1}{12\zeta S+1}\\\
=&\scalebox{0.75}{$\displaystyle\sum_{s=1}^{S}\left[\frac{1}{12(\zeta
Sf^{(s)}\\!+\\!n^{(s)})}\\!-\\!\frac{1}{12\zeta
Sf^{(s)}\\!+\\!1}\right]\\!+\\!\left(\frac{1}{12\zeta
S}\\!-\\!\frac{1}{12\zeta S\\!+\\!1}\right)\\!+\\!\left(\frac{1}{12(\zeta
S\\!+\\!n)}\\!-\\!\frac{1}{12(\zeta
S\\!+\\!n)\\!+\\!1}\right)$}\,.\end{split}$ (7.10)

Compared to the above, observe that summing (over $s$) the bounds from (7.6)
based on (7.1) actually gives strictly tighter bounds:

$\begin{split}&\sum_{s=1}^{S}\left[\frac{1}{12(\zeta
Sf^{(s)}\\!+\\!n^{(s)})\\!+\\!1}\\!-\\!\frac{1}{12\zeta Sf^{(s)}}\right]\\\
<&\sum_{s=1}^{S}\mathscr{R}_{\mathrm{NB}}(\zeta
Sf^{(s)},n^{(s)})=\mathscr{R}_{\mathrm{hNBDM}}(\zeta
S\vec{f},\vec{\mathbf{n}})\\\ <&\sum_{s=1}^{S}\left[\frac{1}{12(\zeta
Sf^{(s)}\\!+\\!n^{(s)})}\\!-\\!\frac{1}{12\zeta
Sf^{(s)}\\!+\\!1}\right]\,.\end{split}$ (7.11)

The additional terms from the bounds in (7.10) not present in the bounds from
(7.11) are the result of twice applying Lemma 3.4 unnecessarily, once to
$\frac{\Gamma(\zeta S+n)}{\Gamma(\zeta S)}$ and once to $\frac{\Gamma(\zeta
S)}{\Gamma(\zeta S+n)}$. Because they cancel each other out, they obviously do
not need to be bounded. Therefore the bounds do not contradict one another.

######  7.2.4 Sketch of Proof that hNBDM Converges in Distribution to hPoMu

Again, it suffices to show that the probability mass function of hNBDM
converges to that of hPoMu. Cf. section 5.2.3.

That (7.8) converges to $1$ as $\zeta\to\infty$ quickly follows from Lemma
3.3. This convergence implies that the hNBDM distribution converges to the
hPoMu distribution as $\zeta\to\infty$.

Equivalently, the convergence of the hNBDM distribution to the hPoMu
distribution as $\zeta\to\infty$ follows from the convergence of (7.7) to $1$
as $\zeta\to\infty$ (replacing $\zeta$, $n$, and $\lambda$ with $\zeta
f^{(s)}$, $n^{(s)}$, and $f^{(s)}\lambda$ respectively). Because (7.8) is the
product of (7.7) and (6.6), it also equivalently follows from the convergence
of both (7.7) and (6.6) to $1$ as $\zeta\to\infty$.

Instead of showing that the likelihood ratio converges to $1$, one can also
show (e.g. according to one’s personal preference) that the logarithm of the
likelihood ratio converges to $0$, using Lemma 3.5 (the “logarithmic version”
of Lemma 3.3). The “fudge function” $\Phi$ is defined for convenience in
equation (3.A.3). Then, starting from (7.8), tedious algebra shows that the
logarithm of the likelihood ratio of hNBDM with respect to hPoMu is

$\lambda-n+\sum_{s=1}^{S}\Phi(\zeta
Sf^{(s)},n^{(s)},n^{(s)}\\!-\\!\frac{1}{2})-\sum_{s=1}^{S}\Phi(\zeta
Sf^{(s)},f^{(s)}\lambda,n^{(s)})+\\!\mathscr{R}_{\mathrm{hNBDM}}(\zeta
S\vec{f},\vec{\mathbf{n}})\,.$ (7.12)

####  8 Working Model for Arbitrary Combinations of Heterogeneities

While the hNBDM family of distributions is useful for modelling density
heterogeneity, it does have the limitation that the value of density
heterogeneity is assumed to be coupled with corresponding value of
compositional heterogeneity. By relaxing this restriction we can generalize
the hNBDM family to the “generalized hNBDM” or “ghNBDM” family of
distributions. This makes it possible to model arbitrary combinations of
density and compositional heterogeneities.

#####  8.1 ghNBDM Family Definition

Under a distribution in the ghNBDM family, for the $d$’th droplet: the
probability for the strain distribution vector given the number of cells
${\mathbb{P}\left(\vec{\mathbf{N}}_{d}(0)=\vec{\mathbf{n}}\nonscript\;\middle|\nonscript\;N_{d}(0)=n\right)}$
is described by a Dirichlet-Multinomial distribution with concentration
parameter $\zeta_{C}$, while the probability for the number of cells
${\mathbb{P}\left(N_{d}(0)=n\right)}$ is described by a Negative binomial
distribution with concentration parameter $\zeta_{D}$. Thus for the
unconditional probability:

$\begin{split}&\mathbb{P}\left(\vec{\mathbf{N}}_{d}(0)=\vec{\mathbf{n}}\right)\\\
=&\quad\mathbb{P}\left(N_{d}(0)=n\right)\\!\cdot\\!\mathbb{P}\left(\vec{\mathbf{N}}_{d}(0)=\vec{\mathbf{n}}\nonscript\;\middle|\nonscript\;N_{d}(0)=n\right)\\\
=&\quad\frac{\Gamma(\zeta_{D}S+n)}{\Gamma(\zeta_{D}S)}\\!\cdot\\!\frac{(\zeta_{D}S)^{\zeta_{D}S}}{(\zeta_{D}S+\lambda)^{\zeta_{D}S+n}}\\!\cdot\\!\frac{\lambda^{n}}{n!}\\!\cdot\\!\frac{\Gamma(\zeta_{C}S)}{\Gamma(\zeta_{C}S+n)}\\!\cdot\\!\left[\prod_{s=1}^{S}\frac{\Gamma(\zeta_{C}Sf^{(s)}+n^{(s)})}{\Gamma\left(\zeta_{C}Sf^{(s)}\right)}\right]\\!\cdot\\!\binom{n}{n^{(1)}\\!\cdots\\!n^{(S)}}\,.\end{split}$
(8.1)

Unlike for the hNBDM family, for arbitrary members of the ghNBDM family one
may have that $\zeta_{C}\not=\zeta_{D}$. In general the marginal distributions
of members of the ghNBDM family are not mutually independent (and in fact have
non-zero cross-covariance) and need not be negative binomial distributed.

The hNBDM family is the special case of the ghNBDM family when
$\zeta_{C}=\zeta_{D}$. The hPoDM and hPoMu distributions emerge as limiting
cases of the ghNBDM family, or are bona fide members if we allow “$\infty$” as
a valid concentration parameter value. Cf. figure 3.6. The case where
$\zeta_{C}$ is allowed to approach infinity while $\zeta_{D}$ remains finite
corresponds to the hierarchical Negative Binomial Multinomial (hNBMu)
distribution.

Figure 3.6: Schematic of the ghNBDM family of distributions, assuming a fixed
$\lambda$ and vector of population strain frequencies $\vec{f}$.

#####  8.2 Over-Dispersions for these Distributions

Since the non-diagonal covariances of the hPoMu distribution are $0$, it of
course does not make sense to speak of over-dispersion relative to them using
the definition (4.2). However, we may still sensibly speak of the over-
dispersion of their marginals with respect to the corresponding marginals of
the hPoMu distribution (which unsurprisingly we will now do).

The over-dispersion of the $s$’th marginal of the ghNBDM distribution with
respect to the $s$’th marginal of the hPoMu distribution is

$\frac{[(1+\zeta_{D}S)+f^{(s)}(\zeta_{C}S-\zeta_{D}S)]}{(1+\zeta_{C}S)(\zeta_{D}S)}\lambda\,.$
(8.2)

Observe how this value depends on the particular $s$. Moreover, the expression
(8.2) above:

  * •

reduces to (8.3) when $\zeta_{D}=\zeta_{C}$,

  * •

approaches (8.4) when $\zeta_{D}\to\infty$,

  * •

approaches (8.5) when $\zeta_{C}\to\infty$,

  * •

approaches $0$ when both $\zeta_{D}\to\infty$ and $\zeta_{C}\to\infty$.

The last corresponds to hPoMu, of course. Details about the remaining special
cases (hNBDM, hPoDM, hNBMu) are given below.

######  8.2.1 hNBDM Over-dispersion

The over-dispersion of the $s$’th marginal of the hNBDM distribution with
respect to the $s$’th marginal of the hPoMu distribution is

$\frac{\lambda}{\zeta S}\,,$ (8.3)

which interestingly does not depend on the particular strain $s$, and moreover
equals the over-dispersion of a Negative Binomial distribution with respect to
a Poisson distribution, cf. (7.2).

The latter is to be expected given that the marginals of the hPoMu
distribution are $\operatorname{Pois}(f^{(s)}\lambda)$ distributed and the
marginals of the hNBDM distribution are $\operatorname{NB}(\zeta
Sf^{(s)},f^{(s)}\lambda)$ distributed. Note that (8.3) also tends to $0$ as
$\zeta\to\infty$ as it should. Finally observe how the over-dispersion (8.3)
is inversely proportional to the concentration parameter $\zeta$, just like
(7.2) which it equals.

######  8.2.2 hPoDM Over-dispersion

The over-dispersion of the $s$’th marginal of the hPoDM distribution with
respect to the $s$’th marginal of the hPoMu distribution is

$\frac{1-f^{(s)}}{1+\zeta_{C}S}\lambda\,.$ (8.4)

Observe how this value depends on the particular $s$. And of course (8.4)
tends to $0$ as $\zeta_{C}\to\infty$. Further observe how the over-dispersion
(8.4) is (asymptotically) proportional to the inverse of the concentration
parameter $\zeta_{C}$, just like (6.2).

######  8.2.3 hNBMu Over-dispersion

The over-dispersion of the $s$’th marginal of the hNBMu distribution with
respect to the $s$’th marginal of the hPoMu distribution is

$\frac{f^{(s)}\lambda}{\zeta_{D}S}\,.$ (8.5)

Observe how this value depends on the particular $s$. And of course (8.5)
tends to $0$ as $\zeta_{D}\to\infty$. Further observe how the over-dispersion
(8.5) is proportional to the inverse of the concentration parameter
$\zeta_{D}$, just like (6.2).

####  9 Conclusion

###### Findings and Contributions

We have identified

  * •

a default working model we can use to for the initial droplet formation,

  * •

key implicit assumptions made by that default working model,

  * •

concrete ways those assumptions could fail to be true in practice, and

  * •

working models generalizing the default working model which can account for
any of those failed assumptions.

###### Practical Implications

We can now start to ask how realistic the default working model is and what
changes might need to be made to use a more realistic working model instead.
Being able to use a more realistic working model means that the data-
generating process for MOREI can be simulated more accurately. Being able to
simulate the data-generating process for MOREI more accurately will ultimately
help to identify best-performing methods for characterizing microbial
interactions based on MOREI data, cf. part III.

###### Next Steps and Open Questions

In chapter 4 we compare the default model with several of these alternative
models to see whether there are noticeable differences. In chapter 5 we look
at the effects of these alternative models on the throughput of data useful
for characterizing microbial interactions. Finally in chapter 6 we propose and
investigate the performance of estimators which will allow us to ascertain
which of these models is most realistic in practice. We want to know how much
the default working model differs from other, potentially more realistic
working models and how to quantify any such differences.

### Chapter 3. Appendices

###  Appendix 3.A Results for Bounding Approximate Likelihood Ratios

######  Lemma 3.1.

Using Stirling’s Approximation as given in [Rob55],

$\begin{split}&\frac{1}{12(n+m)+1}-\frac{1}{12n}\\\
<&\log\left(\left(\frac{(n+m)!}{n!}\right)\middle/\left(e^{-m}\left(\frac{n+m}{n}\right)^{n+m+\frac{1}{2}}n^{m}\right)\right)\\\
<&\frac{1}{12(n+m)}-\frac{1}{12n+1}\,.\end{split}$ (3.A.1)

Proof: This follows basically immediately from the bounds in [Rob55] via
simple but tedious algebraic manipulations. $\square$

Heuristically, this says that $\frac{(n+m)!}{n!}\approx n^{m}$ when $n$ is
large. This is because the expression
${\left(\frac{n+m}{n}\right)^{n+m+\frac{1}{2}}}$ approaches $e^{m}$ as
$n\to\infty$ by Lemma 3.3.

######  Lemma 3.2.

Given $\alpha,c,y_{1},y_{2}\in\mathbb{R}$ (in particular $c$ is constant
w.r.t. $x$):

$\lim_{x\to\infty}\left(\frac{\alpha x+y_{1}}{\alpha x+y_{2}}\right)^{c}=1\,.$

Proof: By properties of limits, since $f(x)=x^{c}$ is continuous, this equals

$\left[\lim_{x\to\infty}\left(\frac{\alpha x+y_{1}}{\alpha
x+y_{2}}\right)\right]^{c}\,,$

which by L’Hôspital’s Rule equals $1^{c}=1$. $\square$

######  Lemma 3.3.

Given $\alpha,c,y_{1},y_{2}\in\mathbb{R}$, $\alpha>0$, one has

$\lim_{x\to\infty}\left(\frac{\alpha x+y_{1}}{\alpha x+y_{2}}\right)^{\alpha
x+c}=e^{y_{1}}e^{-y_{2}}\,.$

Proof: By properties of limits, the LHS of the above equals

$\left[\lim_{x\to\infty}\left(\frac{\alpha x+y_{1}}{\alpha
x+y_{2}}\right)^{\alpha x}\right]\left[\lim_{x\to\infty}\left(\frac{\alpha
x+y_{1}}{\alpha x+y_{2}}\right)^{c}\right]=\lim_{x\to\infty}\left(\frac{\alpha
x+y_{1}}{\alpha x+y_{2}}\right)^{\alpha x}\,,$

where the second equality uses Lemma 3.2. Rearranging further:

$=\lim_{x\to\infty}\left[\left(\left(\frac{x+\frac{y_{1}}{\alpha}}{x}\right)^{x}\right)^{\alpha}\left(\left(\frac{x+\frac{y_{2}}{\alpha}}{x}\right)^{x}\right)^{-\alpha}\right]\,,$

from which the result more or less immediately follows from the properties of
limits and the definition of the exponential function. (In particular because
$f_{1}(x)=x^{\alpha}$ and $f_{2}(x)=x^{-\alpha}$ are continuous.) $\square$

######  Lemma 3.4.

Using the generalization of the bounds on Stirling’s approximation to the
entire $\Gamma$ function as given in Theorem 5 of [Gor94]

$\begin{split}&\frac{1}{12(x+\beta)+1}-\frac{1}{12x}\\\
<&\log\left(\left(\frac{\Gamma(x+\beta)}{\Gamma(x)}\right)\middle/\left(e^{-\beta}\left(\frac{x+\beta}{x}\right)^{x+\beta-\frac{1}{2}}x^{\beta}\right)\right)\\\
<&\frac{1}{12(x+\beta)}-\frac{1}{12x+1}\,.\end{split}$ (3.A.2)

Proof: Lemma 3.4 follows almost immediately from Theorem 5 of [Gor94] via
tedious algebraic manipulations. $\square$

Heuristically, Lemma 3.4 can be interpreted as saying that
$\frac{\Gamma(x+\beta)}{\Gamma(x)}\approx x^{\beta}$ for large $x$. Again,
Lemma 3.3 implies that
${\left(\frac{x+\beta}{x}\right)^{x+\beta-\frac{1}{2}}}$ approaches
$e^{\beta}$ as $x\to\infty$.

The “fudge function” $\Phi$ is defined as

$\Phi(x,y,z):=(x+z)\left[\log(x+y)-\log(x)\right]\,.$ (3.A.3)

######  Lemma 3.5.

Given $y,z\in\mathbb{R}$ which are constant with respect to $x\in\mathbb{R}$,

$\lim_{x\to\infty}\Phi(x,y,z)=y\,.$

Proof: This is a corollary of Lemma 3.3. Specifically, because the logarithm
is continuous and continuous functions commute with limits,

$\begin{split}&\lim_{x\to\infty}\Phi(x,y,z)=\log\left(\lim_{x\to\infty}\exp(\Phi(x,y,z))\right)=\log\left(\lim_{x\to\infty}\left(\frac{x+y}{x}\right)^{x+z}\right)\\\
&=\log(e^{y})=y\,,\end{split}$

where the penultimate equality was an application of 3.3. $\square$

#####  3.A.1 Comments about Robbins-like Bounds

The bounds on $n!$ given in [Rob55] follow from those for the $\Gamma$
function given in Theorem 5 of [Gor94] using the identity
$n!=n\Gamma(n)=n(n-1)!$ (for all $n\in\mathbb{N}$). They are still slightly
different though, because for $n\in\mathbb{N}$ the expression from [Gor94]
gives ${\Gamma(n)=\frac{n!}{n}=(n-1)!\sim\sqrt{2\pi}n^{n-\frac{1}{2}}e^{-n}}$,
whereas directly121212As opposed to indirectly, by dividing the expression for
$n!$ by $n$. applying the expression from [Rob55] to $(n-1)$ gives
${(n-1)!\sim\sqrt{2\pi}(n-1)^{n-\frac{1}{2}}e^{-(n-1)}}$. Although these
expressions are technically distinct, it follows from Lemma 3.3 that they are
asymptotically equivalent ($\sim$, the limit of their ratio approaches $1$ as
$n\to\infty$). Therefore they are consistent with one another.

###  Appendix 3.B Cumulants of Hierarchical Count-Categorical Distributions

The following are consequences of the law of total covariance (of which the
law of total variance is a special case) and in turn from the law of total
expectation.

######  Lemma 3.6.

If $\vec{\mathbf{X}}$ is such that
$\vec{\mathbf{X}}|N\sim\operatorname{Mult}(N,\vec{f})$ or
$\vec{\mathbf{X}}|N\sim\operatorname{DirMult}(\zeta_{C}S,N,\vec{f})$, then for
all $s\in[S]$:

$\mathbb{E}[X_{s}]=f^{(s)}\mathbb{E}[N]\,.$

######  Lemma 3.7.

If $\vec{\mathbf{X}}$ is such that
$\vec{\mathbf{X}}|N\sim\operatorname{Mult}(N,\vec{f})$, then for all
${s\in[S]}$:

$\operatorname{Var}[X_{s}]=(f^{(s)})^{2}\operatorname{Var}[N]+f^{(s)}(1-f^{(s)})\mathbb{E}[N]\,.$
(3.B.1)

######  Lemma 3.8.

If $\vec{\mathbf{X}}$ is such that
$\vec{\mathbf{X}}|N\sim\operatorname{DirMult}(\zeta_{C}S,N,\vec{f})$, then for
all ${s\in[S]}$:

$\operatorname{Var}[X_{s}]=\left[\frac{f^{(s)}(1-f^{(s)})}{1+\zeta_{C}S}+(f^{(s)})^{2}\right]\operatorname{Var}[N]+\frac{f{s}(1-f^{(s)})}{1+\zeta_{C}S}\mathbb{E}[N](\mathbb{E}[N]+\zeta_{C}S)\,.$
(3.B.2)

Observe how (using L’Hôspital’s rule) equation (3.B.2) approaches (3.B.1) as
${\zeta_{C}\to\infty}$.

######  Lemma 3.9.

If $\vec{\mathbf{X}}$ is such that
$\vec{\mathbf{X}}|N\sim\operatorname{Mult}(N,\vec{f})$, then for all
${s_{1},s_{2}\in[S]}$ ($s_{1}\not=s_{2}$):

$\operatorname{Cov}[X_{s_{1}},X_{s_{2}}]=f^{(s_{1})}f^{(s_{2})}\left(\operatorname{Var}[N]-\mathbb{E}[N]\right)\,.$
(3.B.3)

######  Lemma 3.10.

If $\vec{\mathbf{X}}$ is such that
$\vec{\mathbf{X}}|N\sim\operatorname{DirMult}(\zeta_{C}S,N,\vec{f})$, then for
all ${s_{1},s_{2}\in[S]}$ (${s_{1}\not=s_{2}}$):

$\operatorname{Cov}[X_{s_{1}},X_{s_{2}}]=f^{(s_{1})}f^{(s_{2})}\left(\operatorname{Var}[N]-\frac{1}{1+\zeta_{C}S}\left[\zeta_{C}S\mathbb{E}[N]+\mathbb{E}[N^{2}]\right]\right)\,.$
(3.B.4)

Observe how (using L’Hôspital’s rule) equation (3.B.4) approaches (3.B.3) as
${\zeta_{C}\to\infty}$.

#####  3.B.1 Variance and Covariance of ghNBDM

From Lemma 3.8, for the ghNBDM distributions, for each ${s\in[S]}$:

$\operatorname{Var}[N^{(s)}(0)]=f^{(s)}\lambda+\frac{[f^{(s)}(1+\zeta_{D}S)+(f^{(s)})^{2}(\zeta_{C}S-\zeta_{D}S)]}{(1+\zeta_{C}S)(\zeta_{D}S)}\lambda^{2}\,.$
(3.B.5)

Note how (3.B.5) reduces to (3.B.9) when $\zeta_{D}=\zeta_{C}$. Similarly,
(3.B.5) approaches (3.B.11) as $\zeta_{D}\to\infty$, and (3.B.5) approaches
(3.B.13) as $\zeta_{C}\to\infty$, both conclusions following from L’Hôspital’s
rule.

From Lemma 3.10, for the ghNBDM distributions, for each ${s_{1},s_{2}\in[S]}$
(${s_{1}\not=s_{2}}$):

$\operatorname{Cov}[N^{(s_{1})}(0),N^{(s_{2})}(0)]=f^{(s_{1})}f^{(s_{2})}\left[\frac{\zeta_{C}S-\zeta_{D}S}{(\zeta_{D}S)(1+\zeta_{C}S)}\right]\lambda^{2}\,.$
(3.B.6)

Note how (3.B.6) reduces to (3.B.10) when $\zeta_{D}=\zeta_{C}$. Similarly,
(3.B.6) approaches (3.B.12) as $\zeta_{D}\to\infty$, and (3.B.6) approaches
(3.B.14) as $\zeta_{C}\to\infty$, both conclusions again following from
L’Hôspital’s rule. Information about special cases (hPoMu, hNBDM, hPoDM,
hNBMu) is given below.

######  3.B.1.1 Variance and Covariance of hPoMu

For the hPoMu distribution, for each ${s\in[S]}$:

$\operatorname{Var}[N^{(s)}(0)]=f^{(s)}\lambda\,,$ (3.B.7)

and for all ${s_{1},s_{2}\in[S]}$ (${s_{1}\not=s_{2}}$):

$\operatorname{Cov}[N^{(s_{1})}(0),N^{(s_{2})}(0)]=0\,.$ (3.B.8)

This follows either from applying Lemmas 3.7 and 3.9 directly, or recalling
that the marginals of the hPoMu distribution are mutually independent
$\operatorname{Pois}(f^{(s)}\lambda)$ random variables.

######  3.B.1.2 Variance and Covariance of hNBDM

For the hNBDM distribution, for each ${s\in[S]}$:

$\operatorname{Var}[N^{(s)}(0)]=f^{(s)}\lambda+\frac{f^{(s)}}{\zeta
S}\lambda^{2}\,,$ (3.B.9)

and for all ${s_{1},s_{2}\in[S]}$ (${s_{1}\not=s_{2}}$):

$\operatorname{Cov}[N^{(s_{1})}(0),N^{(s_{2})}(0)]=0.$ (3.B.10)

These follow from applying Lemmas 3.8 and 3.10 directly, or recalling that the
marginals of the hNBDM distribution are mutually independent
$\operatorname{NB}(\zeta Sf^{(s)},f^{(s)}\lambda)$ random variables. Notice
also how (3.B.9) approaches (3.B.7), as $\zeta\to\infty$, as it should. Of
course (3.B.10) already equals (3.B.8).

######  3.B.1.3 Variance and Covariance of hPoDM

For the hPoDM distribution, for each ${s\in[S]}$:

$\operatorname{Var}[N^{(s)}(0)]=f^{(s)}\lambda+\frac{f^{(s)}(1-f^{(s)})}{1+\zeta_{C}S}\lambda^{2}\,,$
(3.B.11)

and for all ${s_{1},s_{2}\in[S]}$ (${s_{1}\not=s_{2}}$):

$\operatorname{Cov}[N^{(s_{1})}(0),N^{(s_{2})}(0)]=-\frac{f^{(s_{1})}f^{(s_{2})}}{1+\zeta_{C}S}\lambda^{2}\,.$
(3.B.12)

These follow from applying Lemmas 3.8 and 3.10. Notice also how (3.B.11)
approaches (3.B.7), and (3.B.12) approaches (3.B.8), as $\zeta_{C}\to\infty$,
as it should.

######  3.B.1.4 Variance and Covariance of hNBMu

From Lemma 3.7, for the hNBMu distribution, for each ${s\in[S]}$:

$\operatorname{Var}[N^{(s)}(0)]=f^{(s)}\lambda+\frac{(f^{(s)})^{2}}{\zeta_{D}S}\lambda^{2}\,.$
(3.B.13)

As expected, (3.B.13) approaches (3.B.7) as $\zeta_{D}\to\infty$.

From Lemma 3.9, for the hNBMu distribution, for each ${s_{1},s_{2}\in[S]}$
(${s_{1}\not=s_{2}}$):

$\operatorname{Cov}[N^{(s_{1})}(0),N^{(s_{2})}(0)]=\frac{f^{(s_{1})}f^{(s_{2})}}{\zeta_{D}S}\lambda^{2}\,.$
(3.B.14)

Again, as expected, (3.B.14) approaches (3.B.8) as $\zeta_{D}\to\infty$.

[chapters]

###  Chapter 4 Log-Likelihood Ratio Statistics for Goodness of Fit

Herein I show that log likelihood ratios can be used to determine whether
failure of assumptions from section 3.3 causes substantially new behavior. See
sections 2, 3.2, and 4. Failure of the sampling without replacement assumption
has negligible effects in practice, but failures of the other assumptions
could be important. See sections 4 and 5.

Section 1 provides context for the chapter. Section 2 explains some of the
technical details of the comparisons. Section 3 explains the implementation
details of the comparisons. Section 4 describes and explains the comparisons.
Section 5 interprets the results from section 4, explaining what they mean for
modelling the initial formation of droplets. Section 6 summarizes the findings
and concludes with implications and next steps.

[chapters] [chapters]1

##### Chapter 4 Contents

####  1 Background and Significance

The overall scientific problem remains how to design the experiment so that
enough data is collected from a single experimental run to adequately power
statistical inference, without additional costly and time-consuming
experimental runs being required. Cf. section 8.2. As argued in section 8,
this amounts to knowing the joint distribution of the strain counts well
enough to predict the probabilities associated with any possible combination
of strains. In the previous chapter 3, I identified a simple working model
(hPoMu) for the initial formation of droplets which practitioners are likely
to assume is unquestionably accurate. I explained the assumptions required to
justify such a simple working model, and presented working models resulting
from relaxing some of the most unreasonable of those assumptions.

###### Broader field

The problem considered in this chapter falls under the general field of “model
selection”. See [HTF09, chapter 7] for a more widely accessible introduction
to the general theory of model selection, or [BBM99] for a highly
sophisticated discussion thereof. Cf. the discussion from the introduction to
Part II. Herein we only need to consider the subfield of model selection
corresponding to so-called “goodness of fit” tests, in particular those using
the classic likelihood ratio statistic. For a basic introduction, tailored for
non-statisticians, on the use of the likelihood ratio statistic for model
selection with ecological count data, see [LBG10]. For a more sophisticated
discussion of the likelihood ratio statistic, see [Kee10, chapter 17] or
[Vaa98, chapter 16] for information such as a general definition, or the
statement and proof of the result about its asymptotic $\chi^{2}$
distribution.

###### Specific problem

Herein I investigate whether the failure of any of the implicit hPoMu
assumptions from section 3.3 causes substantially new behavior. As explained
previously, failure of any of these assumptions can be associated with a new
working model that generalizes hPoMu. Therefore, if data generated by a given
such generalized working model fits hPoMu well, it seems reasonable to say
that the generalized working model, and thus also the corresponding assumption
failure, is not causing substantially new behavior.

On the other hand, given two working models which are both adequate for
describing the initial droplet formation process, the simpler working model
will be easier to understand. Therefore, given two working models with the
same descriptive power, we should choose the simpler working model. It follows
that justifying the use of a more complicated working model over a simpler
working model requires showing that the more complicated working model has
substantially more descriptive power than the simpler working model.

Therefore, to justify potentially using any of the other, more complicated
working models in place of hPoMu, we must first ascertain whether any of the
other working models have substantially more descriptive power than hPoMu.
Thus we want to ascertain how well hPoMu can describe data generated by the
other working models. Only those working models which do not always produce
data that can be described well by hPoMu have substantially more descriptive
power than hPoMu. Any other working models are not as useful.

In statistical terms, we know that a working model as simple as hPoMu will
inevitably lead to biased estimators. So we want to investigate whether the
failure of any of these assumptions can cause estimates derived from hPoMu to
become biased enough that they are no longer useful. Thinking of the estimates
produced by each of these working models as an algorithm, in the targeted
learning roadmap [vdLR11] we want to investigate whether we can rely on a
single “default” algorithm based on questionable assumptions.

###### Particular approach

In this chapter I describe the information needed to implement goodness of fit
tests based on the log-likelihood ratio statistic that compare any one of
these more complicated working models against the “default” working model
(hPoMu). After generating data from any of these working models whose
assumptions are more relaxed than those of hPoMu, if these tests demonstrate
that the resulting data are significantly different than what would be
expected under the null model (hPoMu), then this suggests that the predictions
made by hPoMu are inadequate. We could conclude that we may need to pay
attention to whether the corresponding assumption of hPoMu is violated in
practice. In contrast, if the tests demonstrate no significant difference
between the resulting data and what would be expected under the null model
(hPoMu), then this suggests we do not need to care whether the corresponding
assumption of hPoMu is violated in practice.

Comparing the fit of two models for a given dataset by comparing the values of
their likelihoods is a straightforward idea. In turn the (arguably) simplest
way to make that comparison is to consider (the logarithm of) the ratio of
those likelihoods. Because this is exploratory and preliminary work, there
does not seem to be any compelling reason to use approaches that are more
sophisticated than the simplest possible.

For problems similar to these, either the Wald statistic or the Rao score
statistic is often used in place of the likelihood ratio statistic. However,
in addition to being less straightforward, these statistics are also
“equivalent” (in a sense that can be made precise) to the likelihood ratio
statistic. See [Vaa98, section 16.2] or [Kee10, section 17.4] for discussion
of the relationship of the Wald, Rao score, and likelihood ratio statistics.

Similarly, the AIC (Akaike Information Criterion) and BIC (Bayesian
Information Criterion) methods can both (roughly) be thought of as more
sophisticated versions of the likelihood ratio statistic. See, for example,
[BA02], in particular [BA02, section 6.9.3]. Thus the same preference for
simplicity described above is also a reason for not using AIC or BIC.

The “null model” in this case, hPoMu, is (effectively by design) nested within
all of these generalized working models against which it is being compared.
Nevertheless, as explained in [LBG10], the assumption that the null model is
nested within the alternative model is not necessary in general to make
comparisons using the likelihood ratio statistic. So failure of nesting in
other contexts is not a reason to avoid using the likelihood ratio statistic.

In comparing the performances of “Negative Binomial Factor Analysis” (related
to the hNBDM model) and “Poisson Factor Analysis” (related to the hPoMu
model), it can be argued that [Zho18, section 6] effectively addresses a
problem similar to the one considered in this chapter. However, the underlying
questions being asked therein are substantially different from those questions
being asked herein. Hence, there seems to be no reason to expect that the
comparison from [Zho18, section 6] would map (well) to this problem.

####  2 Preliminaries

Below I give the formulae for the log likelihood ratios comparing more complex
working models to hPoMu. These exact formulae for the log likelihood ratios do
not invoke any version of Stirling’s approximation. I used these formulae as
described in section 3.2.

#####  2.1 hTPMH

The likelihood ratio of the hTPMH distribution with respect to the
corresponding hPoMu distribution is

$\mathbb{P}\left(\operatorname{Pois}(\lambda)\leq K_{\scriptscriptstyle
d\\!-\\!1}\right)^{-1}\\!\cdot\\!\frac{(K_{\scriptscriptstyle
d\\!-\\!1}-n)!}{(K_{\scriptscriptstyle
d\\!-\\!1})!}\\!\cdot\\!\left[\prod_{s=1}^{S}\frac{(K^{(s)}_{\scriptscriptstyle
d\\!-\\!1})!}{(K^{(s)}_{\scriptscriptstyle
d\\!-\\!1}-n^{(s)})!}\right]\\!\cdot\\!\left[\prod_{s=1}^{S}(f^{(s)})^{n^{(s)}}\right]^{-1}\,.$
(2.1)

Therefore the log likelihood ratio is

$\begin{split}-\log\left(\mathbb{P}\left(\operatorname{Pois}(\lambda)\leq
K_{\scriptscriptstyle
d\\!-\\!1}\right)\right)\\!+\\!\log\left((K_{\scriptscriptstyle
d\\!-\\!1}-n)!\right)\\!-\\!\log\left((K_{\scriptscriptstyle
d\\!-\\!1})!\right)\\!+\\!\\\
\sum_{s=1}^{S}\left[\log\left((K^{(s)}_{\scriptscriptstyle
d\\!-\\!1})!\right)\\!-\\!\log\left((K^{(s)}_{\scriptscriptstyle
d\\!-\\!1}-n^{(s)})!\right)\right]\\!-\\!\sum_{s=1}^{S}n^{(s)}\log(f^{(s)})\,.\end{split}$
(2.2)

Using the identity $\Gamma(n+1)=n!$, the above can be rewritten using the log-
gamma function131313Which is most likely more numerically stable and
computationally tractable than directly computing the logarithm of the
factorial. (which is directly implemented in SciPy [VGO+20]):

$\begin{split}-\log\left(\mathbb{P}\left(\operatorname{Pois}(\lambda)\leq
K_{\scriptscriptstyle
d\\!-\\!1}\right)\right)\\!+\\!\log\left(\Gamma(K_{\scriptscriptstyle
d\\!-\\!1}-n+1)\right)\\!-\\!\log\left(\Gamma(K_{\scriptscriptstyle
d\\!-\\!1}+1)\right)\\!+\\!\\\
\sum_{s=1}^{S}\left[\log\left(\Gamma(K^{(s)}_{\scriptscriptstyle
d\\!-\\!1}+1)\right)\\!-\\!\log\left(\Gamma(K^{(s)}_{\scriptscriptstyle
d\\!-\\!1}-n^{(s)}+1)\right)\right]\\!-\\!\sum_{s=1}^{S}n^{(s)}\log(f^{(s)})\,.\end{split}$
(2.3)

As an aside, because the above expression equals (5.12) from section 5.2.3,
one can derive an explicit formula for
$\mathscr{R}_{\mathrm{MHg}}(\vec{\mathbf{n}},\vec{\mathbf{K}}_{\scriptscriptstyle
d\\!-\\!1})$. It is not very elegant.

For all observed values of $K_{\scriptscriptstyle d\\!-\\!1}$,
$\mathbb{P}\left(\operatorname{Pois}(\lambda)\leq K_{\scriptscriptstyle
d\\!-\\!1}\right)^{-1}$ was indistinguishable from $1$ up to numerical
precision. (Indeed for $\lambda=2$ that probability appears to be
indistinguishable from $1$ up to numerical precision already for
$K_{\scriptscriptstyle d\\!-\\!1}\geq 20$, values many orders of magnitude
smaller than any observed value.) Therefore the contribution of the logarithm
of this term was treated as if it was exactly $0$ in the code implementation
and not explicitly included.

#####  2.2 hPoDM

The likelihood ratio of the hPoDM distribution with respect to the
corresponding hPoMu distribution is

$\frac{\Gamma(\zeta_{C}S)}{\Gamma(\zeta_{C}S+n)}\\!\cdot\\!\left[\prod_{s=1}^{S}\frac{\Gamma(\zeta_{C}Sf^{(s)}+n^{(s)})}{\Gamma\left(\zeta_{C}Sf^{(s)}\right)}\right]\\!\cdot\\!\left[\prod_{s=1}^{S}(f^{(s)})^{n^{(s)}}\right]^{-1}\,.$
(2.4)

Therefore the log likelihood ratio is

$\displaystyle\log\left(\Gamma(\zeta_{C}S)\right)-\log\left(\Gamma(\zeta_{C}S\\!+\\!n)\right)+\sum_{s=1}^{S}\left[\log\left(\Gamma(\zeta_{C}Sf^{(s)}\\!+\\!n^{(s)})\right)\\!-\\!\log\left(\Gamma\left(\zeta_{C}Sf^{(s)}\right)\right)\right]\\!-\\!\sum_{s=1}^{S}n^{(s)}\log(f^{(s)})\,.$
(2.5)

As an aside, because the above expression equals (6.7) from section 6.2.2, one
can derive an explicit formula for
$\mathscr{R}_{\mathrm{DM}}(\vec{\mathbf{n}},\zeta_{C}S\vec{f})$. Again, it is
not elegant.

#####  2.3 hNBDM

The likelihood ratio of the hNBDM distribution with respect to the
corresponding hPoMu distribution is

$\prod_{s=1}^{S}\left[\frac{\Gamma(\zeta Sf^{(s)}+n^{(s)})}{\Gamma(\zeta
Sf^{(s)})}\\!\cdot\\!\frac{(\zeta Sf^{(s)})^{\zeta Sf^{(s)}}}{(\zeta
Sf^{(s)}+f^{(s)}\lambda)^{\zeta
Sf^{(s)}+n^{(s)}}}\right]\\!\cdot\\!\left[\prod_{s=1}^{S}e^{-f^{(s)}\lambda}\right]^{-1}\,.$
(2.6)

Therefore the log likelihood ratio is

$\begin{split}\sum_{s=1}^{S}\left[\log\left(\Gamma(\zeta
Sf^{(s)}\\!+\\!n^{(s)})\right)-\log\left(\Gamma(\zeta Sf^{(s)})\right)+\zeta
Sf^{(s)}\log\left(\zeta Sf^{(s)}\right)\right.\\\ \left.-(\zeta
Sf^{(s)}+n^{(s)})\log\left(\zeta
Sf^{(s)}+f^{(s)}\lambda\right)\right]+\sum_{s=1}^{S}f^{(s)}\lambda\,.\end{split}$
(2.7)

As an aside, because the above expression equals (7.12) from section 7.2.4,
one can derive an explicit formula for $\mathscr{R}_{\mathrm{hNBDM}}(\zeta
S\vec{f},\vec{\mathbf{n}})$. Again, it is not elegant.

####  3 Methods

Section 3.1 explains which distributions were simulated, why, and how. Section
3.2 explains how logarithms of the likelihood ratios were computed for five of
the distributions. I used NumPy [HMvdW+20] version 1.20.2 and SciPy [VGO+20]
version 1.6.2 for computations, and Matplotlib [Hun07] version 3.4.1 and/or
Seaborn [Was21] version 0.11.1 for plots. Complete implementation details can
be found in the code at the relevant GitLab repository. See
https://gitlab.com/krinsman/droplets.

#####  3.1 Simulated Distributions

To test how much the assumptions of hPoMu may be violated with hPoMu remaining
an adequate model, I simulated five distributions. The first distribution, the
hierarchical truncated Poisson Multivariate Hypergeometric (hTPMH) accounts
for sampling without replacement (section 3.1.1), removing the sampling with
replacement approximation of hPoMu. The remaining four distributions removed
the homogeneity assumption of hPoMu by modelling density and compositional
heterogeneities via over-dispersion (sections 3.1.2 and 3.1.3). Details of the
simulations common to all distributions are explained in section 3.1.4.

######  3.1.1 Simulating Sampling without Replacement

I used an initial total population size of $500,000,000$ ($500$ million) cells
for defining the parameter $K$ of hTPMH. This is substantially less than the
number of cells ($\sim 10^{10}$) which would be sampled from in practice. The
NumPy version of the multivariate hypergeometric distribution does not scale
beyond a total population size of $10^{9}$.

######  3.1.2 Simulating Compositional Heterogeneity Only

The first two distributions model only compositional heterogeneity and have no
density heterogeneity. hPoDM (cf. section 6) with compositional concentration
$\zeta_{C}=100$ has the lowest compositional heterogeneity, followed by hPoDM
with $\zeta_{C}=1$.

######  3.1.3 Simulating Compositional and Density Heterogeneities

The final two distributions model both compositional and density
heterogeneity, with both assumed to correspond to equal over-dispersion. hNBDM
(cf. section 7) with concentration $\zeta:=\zeta_{C}=\zeta_{D}=100$ has the
lowest heterogeneities, followed by hNBDM with $\zeta=1$.

######  3.1.4 Simulation Implementation Details

For all simulations and distributions, the same “simulated community” (or
“SimCom”) of $91$ strains141414Herein I use “strains” to refer equally to
strains belonging to the same species(/genus/family/etc.) as well as to
strains belonging to different species(/genera/families/etc.), because the
distinction is irrelevant for setting up the abstract problem. It may matter
for the implementation of a specific experiment. was used. These correspond to
$90$ strains distributed across $9$ distinct relative abundances or population
frequencies: $10$ strains for each of $.01\%$, $.02\%$, $.05\%$, $.1\%$,
$.2\%$, $.5\%$, $1\%$, $2\%$, and $5\%$. The $91$st strain was a “remainder”
strain with abundance $\approx 12\%$. All distributions also used the same
value of the rate parameter, $\lambda=2$. This equals the expected number of
cells for each droplet for all distributions except151515The expectation of a
truncated Poisson distribution will be less than that of the corresponding
Poisson distribution. In this context any deficit compared to $\lambda=2$ is
negligible in practice, because the support was always truncated to a value no
smaller than $\approx 400$ million. technically for hTPMH.

Each of the distributions was statistically independently simulated $500$
times, with each simulation having $15$ million statistically independently
simulated droplets, except for hTPMH, for which the statistical dependence
structure of successive droplets is inherently Markov. (This corresponds to
$5$ batches, each with $3$ million droplets.)

#####  3.2 Log Likelihood Ratios

I derived explicit formulae for the logarithms of the likelihood ratios
relative to hPoMu for the hTPMH, hPoDM, and hNBDM families of distributions
(cf. section 2). I did not attempt to derive analogous formulae for the
hExhPoDM and hExhNBDM families of distributions to avoid making mistakes when
attempting to compute anti-derivatives of expressions involving the (logarithm
of the) Gamma function. (Cf. section 3.1.)

Then, for each of the five distributions (hTPMH, hPoDM $\zeta_{C}=100$, hNBDM
$\zeta=100$, hPoDM $\zeta_{C}=1$, hNBDM $\zeta=1$) with explicit formulae for
the logarithm of their likelihood ratios with respect to hPoMu, of the $500$
simulations I arbitrarily chose one (the $18$th) and evaluated the respective
formula on all $15$ million droplets. For each of the five distributions, I
also collected the results of doing so into gluttonous groups161616See section
3.4 for clarification of the specific meaning of the term “gluttonous groups”.
See also footnote 22 of section 3.4 regarding the groups defining the
diagonals of the heatmaps. (which by definition overlap) and computed the
geometric mean of the observed likelihood ratios for all groups.

####  4 Results

Section 4.1 discusses evidence indicating that sampling with replacement is
the most innocuous of the assumptions behind hPoMu. Section 4.2 discusses
evidence indicating that ghNBDM distributions with low (but nonzero)
heterogeneity may be adequately modeled by hPoMu. Finally, section 4.3
discusses evidence indicating that ghNBDM distributions with only moderate
heterogeneity can be easily distinguished from hPoMu.

#####  4.1 Sampling without Replacement

All available evidence suggests that, even after $15,000,000$ droplets have
been formed, the hTPMH distribution for the chosen value of $K$ is extremely
similar to the corresponding hPoMu distribution.

(a) Likelihood ratio histogram for hTPMH.

(b) Grouped likelihood ratios for hTPMH.

Figure 4.1: The likelihood ratio with respect to hPoMu never substantially
differs from $1$.

######  4.1.1 Distribution of Likelihood Ratios w.r.t. hPoMu

When plotted on the same scale as that used for the other distributions, the
histogram for hTPMH of observed likelihood ratios relative to hPoMu seen in
figure 1(a) has no right tail whatsoever. This indicates few, if any, of the
droplets were substantially “better explained” using the hTPMH distribution
compared to using hPoMu.

######  4.1.2 Typical Likelihood Ratios as a Function of Strain Frequency

The heat map of typical values for hTPMH of the likelihood ratio across
gluttonous groups showin in figure 1(b) also shows, even at worst, negligible
differences from hPoMu. The values are slightly smaller for the least abundant
strains, reflecting how the sampling with approximation is weakest for those
strains. Even so, the values are still very close to $1$ and well within the
range observed for other distributions.

#####  4.2 Very Similar Distributions

All available evidence suggests that the hPoDM with $\zeta_{C}=100$ and hNBDM
with $\zeta=100$ distributions are very similar to the corresponding hPoMu
distribution. Nevertheless, they also appear to differ slightly more than the
hTPMH distribution.

######  4.2.1 Distributions of Likelihood Ratios w.r.t. hPoMu

Unlike what occurs for hTPMH, figures 2(a) and 2(b) respectively do show the
distributions of observed likelihood ratios for hPoDM with $\zeta_{C}=100$ and
hNBDM with $\zeta=100$ having right tails. However, these right tails are
small.

(a) Compare with figure 1(a).

(b) Compare with figure 1(a).

Figure 4.2: Likelihood ratio histograms for very similar distributions.

######  4.2.2 Typical Likelihood Ratios as a Function of Strain Frequency

Similarly, while figures 3(a) and 3(b) also show greater differences relative
to hPoMu compared to those observed for the hTPMH distribution, the
differences still trend very small. (That the geometric means across the
gluttonous groups tend to be slightly less than $1$, rather than slightly
greater, most likely reflects, in addition to the overlapping nature of the
gluttonous groups, the depletion of droplets with $3$ or more strains relative
to hPoMu that will be observed in these distributions in chapter 5.)

(a) Compare with figure 1(b).

(b) Compare with figure 1(b).

Figure 4.3: Grouped likelihood ratios for very similar distributions.

#####  4.3 Slightly Different Distributions

All available evidence suggests that the hPoDM with $\zeta_{C}=1$ and hNBDM
with $\zeta=1$ distributions are (at least) slightly different from the
corresponding hPoMu distribution.

######  4.3.1 Distributions of Likelihood Ratios w.r.t. hPoMu

Figures 4(a) and 4(b) show very large right tails for the distributions of
observed likelihood ratios for both the hPoDM ($\zeta_{C}=1$) and the hNBDM
($\zeta=1$). This shows that many droplets were substantially “better
explained” using the true simulation distribution rather than hPoMu, something
which did not occur for any of the less heterogeneous distributions. It is
unclear whether either large right tail is “larger” than the other.

(a) Compare with figure 1(a) as well as with figure 2(a).

(b) Compare with figure 1(a) as well as with figure 2(b).

Figure 4.4: Likelihood ratio histograms for slightly similar distributions.

######  4.3.2 Typical Likelihood Ratios as a Function of Strain Frequency

Finally, figures 5(a) and 5(b) show typical values for the likelihood ratios
in each gluttonous group for hPoDM $\zeta_{C}=1$ and hNBDM $\zeta=1$ that
consistently “float above” $1$, in clear contrast to the “hovering around $1$”
trend observed in figures 1(b), 3(a), and 3(b) that correspond to less
heterogeneous distributions. One also sees that the values in figure 5(b)
trend slightly higher than those in figure 5(a), reflecting again how the
hNBDM family of distributions incorporates both kind of heterogeneities
whereas the hPoDM family does not.

(a) Compare with figure 1(b) as well as with figure 3(a).

(b) Compare with figure 1(b) as well as with figure 3(b).

Figure 4.5: Grouped likelihood ratios for slightly similar distributions.

####  5 Discussion

The results show the following. The sampling with replacement assumption of
hPoMu is completely harmless for adequately modelling the data. Depending on
our target parameters, and which values of heterogeneities are realistic, the
homogeneity assumptions of hPoMu could potentially inadequately model the
data. Density and compositional heterogeneities are the only failure of the
hPoMu assumptions with the potential to be important in practice.

Section 5.1 discusses the taxonomy of distributions identified in section 4.
Section 5.2 discusses the two most similar groups, which are either very
difficult or impossible to distinguish from hPoMu. Section 5.3 discusses the
distributions which can be easily distinguished from hPoMu, but whose
differences from hPoMu may not be practically important. Section 6 reviews the
next steps to ascertain which differences from hPoMu might affect our
practical goal of being able to infer microbial interactions.

#####  5.1 Overview of Groups of Distributions

The simulated distributions can be grouped into three levels of assumption
violation.

For the first two levels, the first consisting of hTPMH and the second
consisting of hPoDM $\zeta_{C}=100$ and hNBDM $\zeta=100$, the data produced
is difficult to distinguish from data produced by the hPoMu distribution.
Therefore hPoMu is most likely an adequate model and approximation for data
generated by distributions from these two levels.

The third level consists of hPoDM $\zeta_{C}=1$ and hNBDM $\zeta=1$, which
produce data that can be relatively easily distinguished from data produced by
hPoMu.

#####  5.2 Very Similar Distributions

While hTPMH, hPoDM $\zeta_{C}=100$, and hNBDM $\zeta=100$ are all very similar
to hPoMu, the above results suggest that hTPMH is even more similar to hPoMu
than the other two.

Moreover, the simulations actually overstate how different data generated by
hTPMH should be in practice from data generated by hPoMu. I showed in section
5.2.3 that hTPMH converges in distribution to hPoMu as the size of the total
population of cells being sampled from increases. While the simulation assumed
a population size of $500$ million, or $5\times 10^{8}$, cells, a real MOREI
experiment is likely to sample from a population containing upwards of $10$
billion, or $1\times 10^{10}$, cells, more than twenty times as many.

Thus, to the extent we consider that data generated by hPoDM $\zeta_{C}=100$
and hNBDM $\zeta=100$ are difficult to distinguish from data generated by
hPoMu, data that would be generated by hTPMH in practice should be considered
effectively indistinguishable from data generated by hPoMu. In particular,
even for extremely low abundance strains, these results show that the sampling
with replacement approximation is completely irrelevant in practice to
accurate modelling of the data generated by MOREI.

#####  5.3 Slightly Different Distributions

The results show that there are evident differences between data generated by
distributions in the third level, hPoDM $\zeta_{C}=1$ and hNBDM $\zeta=1$, and
data generated by hPoMu. Yet this still is not enough to be certain that there
are practically significant differences. (Setting aside how to be more precise
about what is a “practically significant” difference.)

For example, the long right tails from figures 4(a) and 4(b) are only easily
visible when plotting the number of counts on a log scale. The total fraction
of droplets represented by those long right tails is actually relatively
small. It’s unclear whether the existence of a small faction of droplets which
are much better explained by the true models than by hPoMu is an important
demerit against hPoMu when most droplets are nearly equally well explained by
both. Focusing on “typical” values of the likelihood ratios as depicted in
figures 5(a) and 5(b), the advantage of the true models over hPoMu is usually
at most around $6\%$.

On the other hand, the “typical” likelihood ratio values in figures figures
5(a) and 5(b) also show a distinct upwards trend as the abundance of one of
the involved strains decreases. As explained in section 5, quantifying the
extent of MOREI’s ability to advance the state of the art requires focusing on
the least abundant strains. Therefore, inasmuch as the discrepancies of these
models compared to hPoMu may be particularly important for the least abundant
strains, the discrepancies may also be important overall.

Thus the results in this chapter seem to point both ways about whether hPoMu
is an adequate working model for data generated by these distributions.

####  6 Conclusion

###### Findings and Contributions

I demonstrated that log likelihood ratios can be used to determine whether the
failure of assumptions from section 3.3 causes substantially new behavior.
Failure of the sampling without replacement assumption turns out to have
negligible effects in practice, but failures of the other assumptions could be
important.

###### Practical Implications

The results from this chapter clarify that we can safely ignore concerns
regarding the sampling without replacement assumption of hPoMu when choosing a
(working) model for the initial formation of droplets. The practical
implications of the results from this chapter are otherwise limited to the
extent that they do not unambiguously indicate whether or not we can ignore
concerns regarding the homogeneity assumptions of hPoMu. These results are
useful, however, for at least raising the possibility that violations of the
homogeneity assumptions could have important effects on our ability to predict
the numbers of droplets available for making inferences about given microbial
interactions.

###### Next Steps and Open Questions

Although these results do suggest that data generated by certain distributions
is unlikely to be adequately modeled by hPoMu, we need to be more precise to
substantiate such a claim. We need to decide precisely which aspects of the
data we want to be model accurately, to target, or otherwise we have no way to
measure whether hPoMu is an adequate model for the data. Only experiments
measuring those targeted aspects of the data will give results which could
definitively determine whether or when hPoMu is an adequate model. This is
chapter 5.

Moreover, even if hPoMu is inaccurate for those targeted aspects of the data,
it could still be an adequate working model as long as its predictions are
conservative compared to the truth, rather than overly optimistic. To make
that assessment, we need to know what values of heterogeneity might be
realistic in practice. That in turn requires having estimators for the
heterogeneities to use on empirical data. We fill that gap later in chapter 6.

[chapters]

###  Chapter 5 [

Goodness of Fit Using Contingency Tables of the Targeted Estimands]Goodness of
Fit Using Contingency Tables of the Targeted Estimands
Effects of Heterogeneities on Data Throughput

Herein I confirm that more severe failures of the assumptions from section 3.3
cause more severe discrepancies with the predictions derived from hPoMu. See
for example sections 4.1.3, 4.1.4, 4.2.1. The nature of the effect depends on
the chosen grouping (picky or gluttonous) of droplets defining the targeted
estimands. See sections 4.2.3 and 5.3. Effects of the failure of the sampling
with replacement assumption are further confirmed to be negligible in
practice. See for example sections 4.1.1 or 4.2.1.

Section 2 explains some of the technical details of the analyses. Section 3
explains implementation details of the analyses. Section 4 explains the
results of the analyses. Section 5 explains how the results are relevant to
modelling the initial formation of droplets.

[chapters] [chapters]1

##### Chapter 5 Contents

####  1 Background and Significance

When planning the previous chapter 4, I had anticipated the effect of the
homogeneity assumptions to be negligible. This was true of the effect of the
sampling with replacement assumption, which was investigated earlier.
Similarity of underlying distributions $\mathcal{P}_{1}$ and $\mathcal{P}_{2}$
should be enough to imply similarity of target parameters
$\boldsymbol{\Psi}(\mathcal{P}_{1})$ and $\boldsymbol{\Psi}(\mathcal{P}_{2})$.
The log-likelihood ratio statistics provide evidence of similarity of the
underlying distributions $\mathcal{P}_{1}$ and $\mathcal{P}_{2}$. Hence, if
the anticipated findings of the previous chapter 4 had been valid, then the
log-likelihood ratio statistics alone would have been sufficient to
demonstrate them.

The questions here are the same as those of the previous chapter 4. In
practice we only care about the relevant portion $\mathcal{Q}_{O}$ of the
probability distribution of the observed data that affects our target
parameters $\boldsymbol{\Psi}(\mathcal{P}_{O})$. (Cf. the overview of these
ideas in [vdLR11, section 1.5].) Again, in the previous chapter 4 we found
evidence that violations of the homogeneity assumptions could cause the entire
underlying probability distribution of the observed data $\mathcal{P}_{O}$ to
substantially differ from that under hPoMu. Yet even if two probability
distributions $\mathcal{P}_{1}$ and $\mathcal{P}_{2}$ substantially differ,
that does not necessarily mean that the relevant portions $\mathcal{Q}_{1}$
and $\mathcal{Q}_{2}$ also substantially differ. Hence, as mentioned before in
section 6, the evidence from chapter 4 is insufficient to say whether
violations of the homogeneity assumptions lead to practically significant
differences, ones that affect $\mathcal{Q}_{O}$ /
$\boldsymbol{\Psi}(\mathcal{P}_{O})$. Thus the goal for this chapter is to
look for and implement goodness of fit procedures that speak more directly to
comparisons of target parameters $\boldsymbol{\Psi}(\mathcal{P})$ and relevant
portions $\mathcal{Q}$ of probability distributions.

###### Broader field

As was the case for the previous chapter 4, the problem of this chapter also
belongs to the general field of model selection, and to the subfield of
goodness of fit tests in particular. Again, the more accessible [HTF09,
chapter 7] or the highly sophisticated [BBM99] are recommended as references
for the general theory of model selection. Cf. the discussion earlier from the
introduction to Part II or from section 1. While the previous chapter 4
discusses goodness of fit tests using the likelihood ratio statistic, this
chapter focuses on goodness of fit tests for contingency tables that use the
$\chi^{2}$-divergence [CZPA09, p. 57] statistic. (Cf. section 2.2 for
terminology.) See the first two chapters of [HCBNM02], or the introduction of
[Joh04], for reviews of the classical procedure for using the
$\chi^{2}$-divergence statistic (“Pearson” $\chi^{2}$ statistic) in goodness
of fit tests for contingency tables. The article [WSCR03] discusses “non-
classical” use of the $\chi^{2}$-divergence statistic for goodness of fit via
Monte Carlo approximations of null distributions, and applications thereof to
ecology. The review [Fie11] discusses the use of classical $\chi^{2}$
procedure as well as other, more sophisticated methods for analyzing
contingency tables (for goodness of fit and for other reasons).

###### Specific problem

Herein “practically significant” differences correspond to the target
parameters defined in section 8.4, the “data throughput” of the experiment.
These target parameters allow us to predict how much data (measured in numbers
of droplets) will be produced that can be used to infer given microbial
interactions, cf. again section 8.2. Thus, in this chapter I investigate how
failures of the implicit hPoMu assumptions, in particular failures of the
homogeneity assumptions, affect these target parameters.

As long as the average number of cells per droplet is not too large, the
number of droplets available to serve as controls will usually be much larger
than the number of droplets available to serve as treatments (i.e. droplets
where strains co-occur). Cf. the heuristic argument in section 2.C. Hence, if
an estimation method is ever too starved of data to make inferences about a
given microbial interaction, most likely the cause will be too few treatment
droplets. In terms of “data throughput”, bottlenecks that occur will usually
be caused by insufficiently many treatment droplets. Therefore the main goal
of this chapter is to understand effects on those target parameters
corresponding specifically to the treatment droplets.

These target parameters, and estimators thereof, can be grouped meaningfully
into contingency tables. (Cf. sections 2.3 and 2.4.) Therefore we can use
goodness of fit procedures designed for contingency tables to accomplish our
goal, to find and implement goodness of fit procedures that speak more
directly to comparisons of these target parameters.

The paper [WSCR03] discusses goodness of fit procedures for ecological
simulations. This includes ecological models that produce count data that can
be grouped into a contingency table, but the authors also show how the same
methodology can be applied more generally. The approach is particularly useful
because it can be applied in situations that violate the assumptions of the
“classical” Pearson goodness of fit procedure for contingency tables.

The paper [Zho18] studies mixed-membership models involving distributions
similar to those in this work. The analysis of mixed-membership models is
known to be related to the analysis of contingency tables, see for example the
review article [Fie11]. The key idea relevant here appears to be that
multivariate count data is closely related to multivariate categorical data.
For finite samples, observed multivariate count data can always be converted
into multivariate categorical data via truncation. In fact, the procedures in
e.g. section 2.3.1 can be thought of as “collapsing” or “marginalizing” a
sparse “multi-dimensional” contingency table with $2^{S}$ cells (“exclusive
and nonempty subsets” [Zho18]) into a “one-dimensional” contingency table
whose number of cells is at most $O(S^{2})$. So although the connection with
the problems and ideas in [Zho18] is certainly currently somewhat hazy, it
also seems reasonable to expect that the connection could be clarified in
future work.

###### Particular approach

Under the assumptions of the classical goodness of fit test for contingency
tables developed by Pearson [HCBNM02], the asymptotic distribution of the
$\chi^{2}$-divergence [CZPA09, p. 57] statistic is the Pearson $\chi^{2}$
distribution [WSCR03, p. 53]. Because this is only an asymptotic result, for
finite samples the sampling distribution of the $\chi^{2}$-divergence
statistic may actually substantially differ from the $\chi^{2}$-distribution,
even under the null hypothesis. Cf. section 3.2. Therefore I computed a Monte
Carlo approximation of the sampling distribution of the $\chi^{2}$-divergence
statistic and used that to compute $p$-values. Cf. [WSCR03, section 3] for a
fairly detailed explanation of Monte Carlo hypothesis testing.

Details of the classical approach, which assumes that the asymptotic
distribution is a valid approximation, again can be found in [HCBNM02]. That
assumption seemed unrealistic here, cf. again section 3.2. Bayesian versions
of this procedure also exist, cf. [Joh04]. However such an approach seemed
unnecessarily complicated for exploratory work like this.

####  2 Preliminaries

In this section I clarify details needed to analyze the data throughput
predicted by each model of initial droplet formation. Section 2.1 gives
precise definitions of the statistics that are relevant for studying data
throughput (following the argument given in section 5.C). Section 2.2 defines
the divergence used to compare distributions. Section 2.3 clarifies the
details of the hypothesis test comparing distributions to hPoMu using this
divergence. In addition to the “picky” and “global” way of grouping of
droplets by strain171717Herein I use “strains” to refer equally to strains
belonging to the same species(/genus/family/etc.) as well as to strains
belonging to different species(/genera/families/etc.), because the distinction
is irrelevant for setting up the abstract problem. It may matter for the
implementation of a specific experiment. for hypothesis testing defined
already in section 2.3, section 2.4 defines a “gluttonous”, “local”, or
“pairwise” way of grouping droplets by strain for hypothesis testing.

#####  2.1 Definition of Throughput

Previously we have been looking at distributions of numbers of cells in a
given droplet (of any strain or of a given strain). As argued in section 5.C,
I believe that understanding the distributions of numbers of strains in a
given droplet is the most relevant for predicting data throughput. In section
2.1.1 I clarify how to define these as statistics computed from the numbers of
cells. Then sections 2.1.2 and 2.1.3 give two different possible definitions
of throughput. This section is a more detailed version of section 8.4.

######  2.1.1 Strain Count Distribution

The notation $\mathbb{I}\\!\left\\{A\right\\}$ denotes the indicator function
for the event $A$.

For a given strain $s$, the strain presence RV (random variable):

$\mathbb{S}^{(s)}(0):=\mathbb{I}\\!\left\\{N^{(s)}(0)\geq
1\right\\}=\mathbb{I}\\!\left\\{N^{(s)}(0)>0\right\\}\,.$ (2.1)

indicates whether any cells of strain $s$ belong to the droplet. It takes
values in $\\{0,1\\}$ like any indicator random variable. Specifically, it
equals $1$ if strain $s$ is present in the droplet and $0$ if strain $s$ is
absent.

Similarly, for a given strain $s$, the strain absence RV:

$\overline{\mathbb{S}}^{(s)}(0):=1-\mathbb{S}^{(s)}(0)=\mathbb{I}\\!\left\\{N^{(s)}(0)=0\right\\}\,,$
(2.2)

takes values in $\\{0,1\\}$ like any indicator random variable, equalling $1$
if strain $s$ is absent in the droplet and $0$ if strain $s$ is present.

For a given droplet, the joint distribution of the strain presence RVs for all
strains corresponds to a binary random vector:

$\vec{\mathbf{\mathbb{S}}}(0):=(\mathbb{S}^{(1)}(0),\dots,\mathbb{S}^{(S)}(0))\,.$
(2.3)

Its values are in $\bigtimes_{s=1}^{S}\\{0,1\\}$ (the $S$-fold Cartesian
product of $\\{0,1\\}$ with itself).

For a given droplet, its strain count $\mathbb{S}(0)$ is the number of strains
present in the droplet. The strain count equals the sum of the strain presence
RVs:

$\mathbb{S}(0):=\sum_{s\in[S]}\mathbb{S}^{(s)}(0)=\sum_{s\in[S]}\mathbb{I}\\!\left\\{N^{(s)}(0)\geq
1\right\\}\,.$ (2.4)

It is an RV taking values in $\\{0\\}\cup[S]$.

Another useful notion is the support of a (non-negative) vector
$\vec{\mathbf{v}}\in\mathbb{R}_{\geq 0}^{S}$:

$\operatorname{supp}(\vec{\mathbf{v}}):=\\{s\in[S]:[\vec{\mathbf{v}}]_{s}>0\\}\subseteq[S]\,.$
(2.5)

It follows directly from the definitions that
$\operatorname{supp}(\vec{\mathbf{N}}(0))=\operatorname{supp}(\vec{\mathbf{\mathbb{S}}}(0))$
always, corresponding to the set of $s\in[S]$ such that
$\mathbb{S}^{(s)}(0)=1$. Moreover we also always have181818 The cardinality
$|\mathcal{S}|$ of a finite set $\mathcal{S}$ is the number of elements of the
set. E.g. ${|[S]|=S}$.  that
$\mathbb{S}(0)=|\operatorname{supp}(\vec{\mathbf{N}}(0))|=|\operatorname{supp}(\vec{\mathbf{\mathbb{S}}}(0))|$.

As argued in section 5.C, I believe that the notion of “throughput” that is
important for this problem corresponds to $\vec{\mathbf{\mathbb{S}}}(0)$ and
not (directly191919Only indirectly as mediated via
$\vec{\mathbf{\mathbb{S}}}(0)$, which of course can be computed from
$\vec{\mathbf{N}}(0)$.) to $\vec{\mathbf{N}}(0)$. While choosing
$\lambda\approx 2$ is intended to make $\vec{\mathbf{N}}(0)$ as similar to
$\vec{\mathbf{\mathbb{S}}}(0)$ as possible, they are definitely not the same,
something which e.g. the multiple representatives problem from section 5.C
makes clear.

Even when the first moments of $\vec{\mathbf{N}}(0)$ remain the same, changes
to the second moments of $\vec{\mathbf{N}}(0)$ can still cause changes in the
expected numbers of droplets serving as treatments or controls. A heuristic
way to understand this is to observe that the expected data throughput
corresponds to the second (or higher order) moments of
$\vec{\mathbf{\mathbb{S}}}(0)$; see sections 2.1.2 and 2.1.3. Although for any
distribution belonging to the ghNBDM working model (see section 8.1) the
expected cell counts for each strain are the same as those under hPoMu, the
choice of distribution from the ghNBDM working model still affects the
predicted data throughput.

Figure 5.1: Gluttonous vs. Picky Definitions of Data Throughput. Cf. text.

######  2.1.2 Gluttonous Definition of Data Throughput

Recall from section 3, or the introduction to Part I, that the treatment group
for inferring the effect of strain $s_{1}$ on the growth of strain $s_{2}$
corresponds to droplets where $s_{1}$ and $s_{2}$ co-occur. The set of all
droplets satisfying this condition is

$\begin{split}\mathcal{T}^{(s_{1},s_{2})}_{g}:=&\left\\{d\in[{D}]:\mathbb{S}^{(s_{1})}_{d}(0)\cdot\mathbb{S}^{(s_{2})}_{d}(0)=1\right\\}\\\
=&\left\\{d\in[{D}]:\mathbb{I}\\!\left\\{N^{(s_{1})}_{d}(0)\geq 1\text{ \&
}N^{(s_{2})}_{d}(0)\geq 1\right\\}=1\right\\}\\\
=&\left\\{d\in[{D}]:\operatorname{supp}(\vec{\mathbf{\mathbb{S}}}_{d}(0))\supseteq\\{s_{1},s_{2}\\}\right\\}\,.\end{split}$
(2.6)

This is called the “gluttonous” treatment group because all droplets where
$s_{1}$ and $s_{2}$ co-occur are included. Cf. figure 5.1.

Recall also from section 3, or the introduction to Part I, that the control
group for inferring the effect of strain $s_{1}$ on the growth of strain
$s_{2}$ corresponds to droplets where $s_{2}$ occurs but $s_{1}$ does not. The
set of all droplets satisfying this condition is

$\begin{split}\mathcal{C}^{(s_{1},s_{2})}_{g}:=&\left\\{d\in[{D}]:\overline{\mathbb{S}}^{(s_{1})}_{d}(0)\cdot\mathbb{S}^{(s_{2})}_{d}(0)=1\right\\}\\\
=&\left\\{d\in[{D}]:\mathbb{I}\\!\left\\{N^{(s_{1})}_{d}(0)=0\text{ \&
}N^{(s_{2})}_{d}(0)\geq 1\right\\}=1\right\\}\\\
=&\left\\{d\in[{D}]:\operatorname{supp}(\vec{\mathbf{\mathbb{S}}}_{d}(0))\supseteq\\{s_{2}\\}\,,\,s_{1}\not\in\operatorname{supp}(\vec{\mathbf{\mathbb{S}}}_{d}(0))\right\\}\,.\end{split}$
(2.7)

This is called the “gluttonous” control group because all droplets where
$s_{2}$ occurs but $s_{1}$ does not are included. Cf. figure 5.1.

These definitions place no constraints on the presence or absence of strains
$s\in[S]$ besides $s_{1}$ and $s_{2}$. This can be favorable e.g. when
studying very rare strains, for which there may be no or very few droplets
containing $s_{1}$ and $s_{2}$ only, but some containing $s_{1}$ and $s_{2}$
along with other strains. These definitions give the best chance of avoiding
“data starvation”.

From equation (2.6) we get (assuming the droplets are identically distributed)
that the expected number of droplets for the gluttonous treatment group for
inferring the effect of strain $s_{2}$ on the growth of $s_{1}$ is

$\mathbb{E}\left[\left\nonscript\;\middle|\nonscript\;\mathcal{T}^{(s_{1},s_{2})}_{g}\right\nonscript\;\middle|\nonscript\;\right]={D}\cdot\mathbb{E}\left[\mathbb{S}^{(s_{1})}(0)\
\mathbb{S}^{(s_{2})}(0)\right]\,.$ (2.8)

Similarly, from equation (2.7) we get (assuming the droplets are identically
distributed) that the expected number of droplets for the corresponding
gluttonous control group for inferring the effect of strain $s_{2}$ on the
growth of $s_{1}$ is

$\mathbb{E}\left[\left\nonscript\;\middle|\nonscript\;\mathcal{C}^{(s_{1},s_{2})}_{g}\right\nonscript\;\middle|\nonscript\;\right]={D}\cdot\mathbb{E}\left[\overline{\mathbb{S}}^{(s_{1})}(0)\cdot\mathbb{S}^{(s_{2})}(0)\right]={D}\cdot\mathbb{E}\left[\mathbb{S}^{(s_{2})}(0)-\mathbb{S}^{(s_{1})}(0)\cdot\mathbb{S}^{(s_{2})}(0)\right]\,.$
(2.9)

As claimed in section 2.1.1, these expressions for the expected gluttonous
data throughput involve the second moments of $\vec{\mathbf{\mathbb{S}}}(0)$.

######  2.1.3 Picky Definition of Data Throughput

Comparing figures 5.1 and 2, we see that the above “gluttonous” definitions
correspond to more possible combinations of strains than we probably would
have included when studying the effect of strain $s_{1}$ on the growth of
strain $s_{2}$ by manually plating cells. This motivates the following
definitions.

The set of droplets where $s_{1}$ and $s_{2}$ co-occur in a way corresponding
to the combination of strains we would use if we were manually plating cells
is

$\begin{split}\mathcal{T}^{(s_{1},s_{2})}_{p}:=&\left\\{d\in[{D}]:\mathbb{S}^{(s_{1})}_{d}(0)\cdot\mathbb{S}^{(s_{2})}_{d}(0)\prod_{\sigma\not\in\\{s_{1},s_{2}\\}}\overline{\mathbb{S}}^{(\sigma)}_{d}(0)=1\right\\}\\\
=&\left\\{d\in[{D}]:\mathbb{I}\\!\left\\{N^{(s_{1})}_{d}(0)\geq 1\text{ \&
}N^{(s_{2})}_{d}(0)\geq 1\text{ \& }\forall\sigma\not\in\\{s_{1},s_{2}\\},\
N^{(\sigma)}_{d}(0)=0\right\\}=1\right\\}\\\
=&\left\\{d\in[{D}]:\operatorname{supp}(\vec{\mathbf{\mathbb{S}}}_{d}(0))=\\{s_{1},s_{2}\\}\right\\}\,.\end{split}$
(2.10)

This is called the “picky” treatment group because droplets where $s_{1}$ and
$s_{2}$ co-occur are included only when all other strains are absent. Cf.
figure 5.1.

The set of droplets where only $s_{2}$ occurs, corresponding to the
combination of strains we would use if we were manually plating cells, is

$\begin{split}\mathcal{C}^{(s_{1},s_{2})}_{p}:=&\left\\{d\in[{D}]:\overline{\mathbb{S}}^{(s_{1})}_{d}(0)\cdot\mathbb{S}^{(s_{2})}_{d}(0)\prod_{\sigma\not\in\\{s_{1},s_{2}\\}}\overline{\mathbb{S}}^{(\sigma)}_{d}(0)=1\right\\}\\\
=&\left\\{d\in[{D}]:\mathbb{I}\\!\left\\{N^{(s_{1})}_{d}(0)=0\text{ \&
}N^{(s_{2})}_{d}(0)\geq 1\text{ \& }\forall\sigma\not\in\\{s_{1},s_{2}\\},\
N^{(\sigma)}_{d}(0)=0\right\\}=1\right\\}\\\
=&\left\\{d\in[{D}]:\operatorname{supp}(\vec{\mathbf{\mathbb{S}}}_{d}(0))=\\{s_{2}\\}\right\\}\,.\end{split}$
(2.11)

This is called the “picky” control group because droplets where $s_{2}$ occurs
but $s_{1}$ does not occur belong only when all other strains are absent. It
is not enough for $s_{1}$ alone to be absent. Cf. figure 5.1.

Unlike the gluttonous groups, the picky groups do place constraints on the
presence or absence of strains $s\in[S]$ besides $s_{1}$ and $s_{2}$. For both
the treatments and the controls, for “picky” groups all strains besides
$s_{1}$ or $s_{2}$ must have zero counts. Every picky group is by definition a
subset of its gluttonous counterpart.

The picky definitions can be favorable when we have plenty of droplets from
which to make estimates. They reduce the possibility of confounding effects on
the growth of strain $s_{2}$ that could be caused by strains that are not
strain $s_{1}$. On the other hand, in instances where there are only very few
or no droplets without “extra” strains, the picky definitions could lead to
“data starvation”.

From equation (2.10) we get (assuming the droplets are identically
distributed) that the expected number of droplets for the picky treatment
group for inferring the effect of strain $s_{2}$ on the growth of $s_{1}$ is

$\mathbb{E}\left[\left\nonscript\;\middle|\nonscript\;\mathcal{T}^{(s_{1},s_{2})}_{p}\right\nonscript\;\middle|\nonscript\;\right]={D}\cdot\mathbb{E}\left[\mathbb{S}^{(s_{1})}(0)\cdot\mathbb{S}^{(s_{2})}(0)\prod_{\sigma\not\in\\{s_{1},s_{2}\\}}\left(1-\mathbb{S}^{(\sigma)}(0)\right)\right]\,.$
(2.12)

Similarly, from equation (2.11) we get (assuming the droplets are identically
distributed) that the expected number of droplets for the corresponding picky
control group for inferring the effect of strain $s_{2}$ on the growth of
$s_{1}$ is

$\mathbb{E}\left[\left\nonscript\;\middle|\nonscript\;\mathcal{C}^{(s_{1},s_{2})}_{p}\right\nonscript\;\middle|\nonscript\;\right]={D}\cdot\mathbb{E}\left[\left(1-\mathbb{S}^{(s_{1})}(0)\right)\cdot\mathbb{S}^{(s_{2})}(0)\prod_{\sigma\not\in\\{s_{1},s_{2}\\}}\left(1-\mathbb{S}^{(\sigma)}(0)\right)\right]\,.$
(2.13)

As claimed in section 2.1.1, these expressions for the expected picky data
throughput involve the higher order moments of $\vec{\mathbf{\mathbb{S}}}(0)$.

#####  2.2 Pearson Categorical Divergence

The Pearson categorical divergence (more often referred to as the “Pearson
$\chi^{2}$ divergence” or “Pearson $\chi^{2}$ statistic”) of
$\vec{\mathbf{v}}=(v^{(1)},\dots,v^{(I)})\in\mathbb{R}_{\geq 0}$ relative to
$\vec{\mathbf{w}}=(w^{(1)},\dots,w^{(I)})\in\mathbb{R}_{\geq 0}$ is
defined[CZPA09, p. 57] as

$\sum_{i=1}^{I}\frac{(v^{(i)}-w^{(i)})^{2}}{w^{(i)}}\,,$ (2.14)

where both $\vec{\mathbf{v}}$ and $\vec{\mathbf{w}}$ belong to the unit
simplex (entries sum to $1$ and are all $\geq 0$), and $\vec{\mathbf{w}}$
belongs to the interior of the unit simplex (entries are all $>0$). This
definition is not symmetric. In general the divergence of $\vec{\mathbf{w}}$
relative to $\vec{\mathbf{v}}$ differs from that of $\vec{\mathbf{v}}$
relative to $\vec{\mathbf{w}}$.

Herein I call this the “Pearson categorical divergence” to avoid confusion
with either

  1. (a)

the “Pearson divergence”, ${1-\rho_{\vec{x}\vec{y}}}$, where
$\rho_{\vec{x}\vec{y}}$ is the Pearson correlation of $\vec{x}$ and $\vec{y}$,
or

  2. (b)

the “Pearson $\chi^{2}$ distribution”202020Technically a parameterized family
of distributions and not a single distribution., the asymptotic sampling
distribution of this statistic under the null distribution of “Pearson-style”
goodness of fit tests.

In particular, it would sound confusing to say (in section 3.2) that the
sampling distribution of the “$\chi^{2}$ divergences” is not a “$\chi^{2}$
distribution”.

#####  2.3 Global Picky Goodness of Fit Test

The strains present in a droplet determine the microbial interactions it can
help characterize. Therefore we want to group the droplets according to the
strains that are present in each droplet. Different models for the initial
formation of droplets make different predictions about which such groups of
droplets will be most common. A hypothesis test whose null distribution is a
sampling distribution derived from the hPoMu working model compares the
predictions made by other (working) models with those made by hPoMu. Rejecting
the null hypothesis corresponds to asserting that hPoMu is unable to describe
the observed data, so the true model that generated the data has more
“descriptive power”, cf. section 1.

The definition of the hypothesis test is split into several parts. Section
2.3.1 defines how the droplets are grouped according to strains for the
hypothesis test. Section 2.3.2 gives probabilities for each group under the
null hypothesis (hPoMu), while section 2.3.3 clarifies the expected counts.
Section 2.3.4 defines the observed counts used in defining the goodness of fit
test. Section 2.3.5 explains why some groups are redundant. Finally, section
2.3.6 uses the conclusions of the previous sections to define the test
statistic.

######  2.3.1 Categories

There are technically exponentially many, $2^{S}$, potential combinations of
strains that could occur in a given droplet. All of these could be used in
defining mutually exclusive categories for a goodness of fit test. However,
for both reasons of computational feasibility and data availability212121 Most
combinations of three or more strains are highly unlikely and will be
represented by few or no data points. Cf. the related discussion in sections
2, 8.5, and 7.3. , herein only combinations with two or fewer strains are
given distinct categories, with all combinations of three or more strains
lumped into a category.

Thus the number of categories for the “global” goodness of fit test is

$1+S+\binom{S}{2}+1\,.$ (2.15)

This is one category for all droplets with zero strains (empty droplets), $S$
categories for each of the picky control groups (droplets with exactly one
strain), $\binom{S}{2}$ categories for each of the picky treatment groups
(droplets with exactly two strains), and one category for all droplets with
three or more (“multiple”) strains. The groups being “picky” and not
“gluttonous” ensures that categories do not overlap. This makes it feasible to
define a single “global” test.

Cf. section 2.1 for more details about “picky” and “gluttonous” groups.

######  2.3.2 Probabilities

The probabilities associated with each category are

$\begin{split}p^{(\emptyset)}_{\mathcal{G}}&:=\mathbb{P}(N(0)=0)=e^{-\lambda}\,,\\\
p^{(s)}_{\mathcal{G}}&:=\mathbb{P}(N^{(s)}(0)\geq
1,\forall\sigma\not=s:N^{(\sigma)}(0)=0)\\\ &\,=\mathbb{P}(N^{(s)}(0)\geq
1)\cdot\prod_{\sigma\not=s}\mathbb{P}(N^{(\sigma)}(0)=0)\\\
&\,=(1-e^{-f^{(s)}\lambda})\cdot e^{-(1-f^{(s)})\lambda}\,,\\\
p^{(s_{1},s_{2})}_{\mathcal{G}}&:=\mathbb{P}(N^{(s_{1})}(0)\geq
1,N^{(s_{2})}(0)\geq
1,\forall\sigma\not\in\\{s_{1},s_{2}\\}:N^{(\sigma)}(0)=0)\\\
&\,=\mathbb{P}(N^{(s_{1})}(0)\geq 1)\cdot\mathbb{P}(N^{(s_{2})}(0)\geq
1)\cdot\prod_{\sigma\not\in\\{s_{1},s_{2}\\}}\mathbb{P}(N^{(\sigma)}(0)=0)\\\
&\,=(1-e^{-f^{(s_{1})}\lambda})(1-e^{-f^{(s_{2})}\lambda})\cdot
e^{-(1-(f^{(s_{1})}+f^{(s_{2})}))\lambda}\,,\\\
p^{(3+)}_{\mathcal{G}}&:=1-p^{(\emptyset)}_{\mathcal{G}}-\sum_{{\sigma}=1}^{S}p^{(\sigma)}_{\mathcal{G}}-\sum_{{\sigma_{1}}=1}^{S}\sum_{{\sigma_{2}}>{\sigma_{1}}}p^{({\sigma_{1}},{\sigma_{2}})}_{\mathcal{G}}\,.\end{split}$
(2.16)

######  2.3.3 Expected Counts

The expected counts equal the total number of droplets ${D}$ multiplied by the
corresponding probabilities.

$\begin{split}M^{(\emptyset)}_{\mathcal{G}}&:={D}\cdot
p^{(\emptyset)}_{\mathcal{G}}\,,\\\ M^{(s)}_{\mathcal{G}}&:={D}\cdot
p^{(s)}_{\mathcal{G}}\,,\\\ M^{(s_{1},s_{2})}_{\mathcal{G}}&:={D}\cdot
p^{(s_{1},s_{2})}_{\mathcal{G}}\,,\\\ M^{(3+)}_{\mathcal{G}}&:={D}\cdot
p^{(3+)}_{\mathcal{G}}\,.\end{split}$ (2.17)

######  2.3.4 Observed Counts

The observed counts for the global goodness of fit test are

$\begin{split}\hat{M}^{(\emptyset)}_{\mathcal{G}}&:=\sum_{d\in[{D}]}\mathbb{I}\\!\left\\{N_{d}(0)=0\right\\}\,,\\\
\hat{M}^{(s)}_{\mathcal{G}}&:=\sum_{d\in[{D}]}\mathbb{I}\\!\left\\{N^{(s)}_{d}(0)\geq
1\right\\}\cdot\prod_{{\sigma}\not=s}\mathbb{I}\\!\left\\{N^{(\sigma)}_{d}(0)=0\right\\}\,,\\\
\hat{M}^{(s_{1},s_{2})}_{\mathcal{G}}&:=\sum_{d\in[{D}]}\left[\vphantom{\prod_{{\sigma}\not\in\\{s_{1},s_{2}\\}}}\mathbb{I}\\!\left\\{N^{(s_{1})}_{d}(0)\geq
1\right\\}\cdot\mathbb{I}\\!\left\\{N^{(s_{2})}_{d}(0)\geq
1\right\\}\right.\\\
&\phantom{:=\sum_{d\in[{D}]}}\left.\cdot\prod_{{\sigma}\not\in\\{s_{1},s_{2}\\}}\mathbb{I}\\!\left\\{N^{(\sigma)}_{d}(0)=0\right\\}\right]\,,\\\
\hat{M}^{(3+)}_{\mathcal{G}}&:={D}-\hat{M}^{(\emptyset)}_{\mathcal{G}}-\sum_{{\sigma}=1}^{S}\hat{M}^{(\sigma)}_{\mathcal{G}}-\sum_{{\sigma_{1}}=1}^{S}\sum_{{\sigma_{2}>\sigma_{1}}}\hat{M}^{(\sigma_{1},\sigma_{2})}_{\mathcal{G}}\,.\end{split}$
(2.18)

######  2.3.5 Symmetry of Definitions

Observe how, for any given pair of strains $s_{1},s_{2}$:

$\begin{array}[]{rcl}p^{(s_{1},s_{2})}_{\mathcal{G}}&=&p^{(s_{2},s_{1})}_{\mathcal{G}}\,,\\\
M^{(s_{1},s_{2})}_{\mathcal{G}}&=&M^{(s_{2},s_{1})}_{\mathcal{G}}\,,\\\
\hat{M}^{(s_{1},s_{2})}_{\mathcal{G}}&=&\hat{M}^{(s_{2},s_{1})}_{\mathcal{G}}\,.\end{array}$
(2.19)

Thus for any given pair of strains $s_{1},s_{2}$, only the left hand sides of
(2.19) are considered, avoiding redundancy and overlapping categories.

######  2.3.6 Test Statistic

The test statistic is the Pearson categorical divergence of

$\begin{split}&{D}^{-1}\\!\cdot\\!\left(\hat{M}^{(\emptyset)}_{\mathcal{G}},\hat{M}^{(1)}_{\mathcal{G}},\dots,\hat{M}^{(S)}_{\mathcal{G}},\hat{M}^{(1,2)}_{\mathcal{G}},\dots,\hat{M}^{(S\\!-\\!1,S)}_{\mathcal{G}},\hat{M}^{(3+)}_{\mathcal{G}}\right)\\\
&\text{relative to}\\\
&{D}^{-1}\\!\cdot\\!\left(M^{(\emptyset)}_{\mathcal{G}},M^{(1)}_{\mathcal{G}},\dots,M^{(S)}_{\mathcal{G}},M^{(1,2)}_{\mathcal{G}},\dots,M^{(S\\!-\\!1,S)}_{\mathcal{G}},M^{(3+)}_{\mathcal{G}}\right)\,.\end{split}$
(2.20)

We can get an asymptotic $p$-value in the standard way, via the survival
function of the $\chi^{2}$ distribution with
$\left(1+S+\binom{S}{2}+1\right)-1$ degrees of freedom.

#####  2.4 Pairwise Gluttonous Goodness of Fit Tests

In section 2.3 droplets were grouped according to strain in a non-overlapping
way. This made it possible to compare other distributions to hPoMu using a
single hypothesis test, but for practical reasons it also required lumping
together all droplets with three or more strains. However, if we are willing
to perform multiple hypothesis tests, we can better utilize the information
contained in those droplets. By performing a separate hypothesis test for each
pair of strains, we can include more droplets in both our “treatment” and
“control” groups.

The definition of the hypothesis test corresponding to each pair of strains is
split into several parts. Section 2.4.1 defines how the droplets are grouped.
Section 2.4.2 gives probabilities for each group under the null hypothesis
(hPoMu), while section 2.4.3 clarifies the expected counts. Section 2.4.4
defines the observed counts used in defining the goodness of fit test. Section
2.4.5 explains why some of the groups are redundant. Finally, section 2.4.6
combines the information from the previous sections to define the test
statistic.

######  2.4.1 Categories

Note that the above definitions correspond to “gluttonous” groups, unlike the
“global” test defined before which corresponds to “picky” groups. The
“gluttony” means that the categories for distinct combinations of strains can
overlap. This is what necessitates separate tests for each combination of
strains, i.e. why they are “pairwise”.

######  2.4.2 Probabilities

The probabilities $p^{(s_{1},s_{2})}_{0,0}$,$p^{(s_{1},s_{2})}_{1,0}$,
$p^{(s_{1},s_{2})}_{1,1}$, $p^{(s_{1},s_{2})}_{0,1}$, for the pairwise tests
are

$\begin{split}p^{(s_{1},s_{2})}_{0,0}&:=\mathbb{P}_{\scalebox{0.75}{$\mathrm{hPoMu}(\lambda,\vec{f})$}}\left(\scalebox{0.85}{$N^{(s_{1})}(0)=0,N^{(s_{2})}(0)=0$}\right)\\\
&=e^{-f^{(s_{1})}\lambda}\cdot e^{-f^{(s_{2})}\lambda}\,,\\\
p^{(s_{1},s_{2})}_{1,0}&:=\mathbb{P}_{\scalebox{0.75}{$\mathrm{hPoMu}(\lambda,\vec{f})$}}\left(\scalebox{0.85}{$N^{(s_{1})}(0)\geq
1,N^{(s_{2})}(0)=0$}\right)\\\ &=(1-e^{-f^{(s_{1})}\lambda})\cdot\cdot
e^{-f^{(s_{2})}\lambda}\,,\\\
p^{(s_{1},s_{2})}_{1,1}&:=\mathbb{P}_{\scalebox{0.75}{$\mathrm{hPoMu}(\lambda,\vec{f})$}}\left(\scalebox{0.85}{$N^{(s_{1})}(0)\geq
1,N^{(s_{2})}(0)\geq 1$}\right)\\\
&=(1-e^{-f^{(s_{1})}\lambda})\cdot(1-e^{-f^{(s_{2})}\lambda})\,,\\\
p^{(s_{1},s_{2})}_{0,1}&:=\mathbb{P}_{\scalebox{0.75}{$\mathrm{hPoMu}(\lambda,\vec{f})$}}\left(\scalebox{0.85}{$N^{(s_{1})}(0)=0,N^{(s_{2})}(0)\geq
1$}\right)\\\
&=e^{-f^{(s_{1})}\lambda}\cdot(1-e^{-f^{(s_{2})}\lambda})\,.\end{split}$
(2.21)

Note the implicit dependence on an assumed $\lambda$ in the definitions of
$p^{(s_{1},s_{2})}_{0,0}$, $p^{(s_{1},s_{2})}_{1,0}$,
$p^{(s_{1},s_{2})}_{1,1}$, $p^{(s_{1},s_{2})}_{0,1}$. Thus the scientist needs
an estimate, or some a priori knowledge, of the value of $\lambda$ to use
these tests.

######  2.4.3 Expected Counts

The expected counts for the pairwise tests are

$\begin{split}M^{(s_{1},s_{2})}_{0,0}&:={D}\cdot p^{(s_{1},s_{2})}_{0,0}\,,\\\
M^{(s_{1},s_{2})}_{1,0}&:={D}\cdot p^{(s_{1},s_{2})}_{1,0}\,,\\\
M^{(s_{1},s_{2})}_{1,1}&:={D}\cdot p^{(s_{1},s_{2})}_{1,1}\,,\\\
M^{(s_{1},s_{2})}_{0,1}&:={D}\cdot p^{(s_{1},s_{2})}_{0,1}\,.\end{split}$
(2.22)

######  2.4.4 Observed Counts

The observed counts are defined as

$\begin{split}\hat{M}^{(s_{1},s_{2})}_{0,0}&:=\left|\left\\{d:N^{(s_{1})}_{d}(0)=0,N^{(s_{2})}_{d}(0)=0\right\\}\right|\,,\\\
\hat{M}^{(s_{1},s_{2})}_{1,0}&:=\left|\left\\{d:N^{(s_{1})}_{d}(0)\geq
1,N^{(s_{2})}_{d}(0)=0\right\\}\right|\,,\\\
\hat{M}^{(s_{1},s_{2})}_{1,1}&:=\left|\left\\{d:N^{(s_{1})}_{d}(0)\geq
1,N^{(s_{2})}_{d}(0)\geq 1\right\\}\right|\,,\\\
\hat{M}^{(s_{1},s_{2})}_{0,1}&:=\left|\left\\{d:N^{(s_{1})}_{d}(0)=0,N^{(s_{2})}_{d}(0)\geq
1\right\\}\right|\,.\end{split}$ (2.23)

######  2.4.5 Symmetry of Definitions

Of course,

$\begin{array}[]{rcl}\hat{M}^{(s_{1},s_{2})}_{1,0}=\hat{M}^{(s_{2},s_{1})}_{0,1}\,,&&p^{(s_{1},s_{2})}_{1,0}=p^{(s_{2},s_{1})}_{0,1}\,,\\\
\hat{M}^{(s_{1},s_{2})}_{0,0}=\hat{M}^{(s_{2},s_{1})}_{0,0}\,,&&p^{(s_{1},s_{2})}_{0,0}=p^{(s_{2},s_{1})}_{0,0}\,,\\\
\hat{M}^{(s_{1},s_{2})}_{1,1}=\hat{M}^{(s_{2},s_{1})}_{1,1}\,,&&p^{(s_{1},s_{2})}_{1,1}=p^{(s_{2},s_{1})}_{1,1}\,.\end{array}$
(2.24)

Therefore the results of the test for $(s_{2},s_{1})$ will always be the same
as those for $(s_{1},s_{2})$, and so at most one of the two tests needs to be
performed.

######  2.4.6 Test Statistic

The test statistic is the Pearson categorical divergence of

$\begin{split}&{D}^{-1}\\!\cdot\\!\left(\hat{M}^{(s_{1},s_{2})}_{0,0},\hat{M}^{(s_{1},s_{2})}_{1,0},\hat{M}^{(s_{1},s_{2})}_{1,1},\hat{M}^{(s_{1},s_{2})}_{0,1}\right)\\\
\text{relative
to}\quad&{D}^{-1}\\!\cdot\\!\left(M^{(s_{1},s_{2})}_{0,0},M^{(s_{1},s_{2})}_{1,0},M^{(s_{1},s_{2})}_{1,1},M^{(s_{1},s_{2})}_{0,1}\right)\,,\end{split}$
(2.25)

namely

$\sum_{i=0}^{1}\sum_{j=0}^{1}\frac{\left(\hat{M}^{(s_{1},s_{2})}_{i,j}-M^{(s_{1},s_{2})}_{i,j}\right)^{2}}{M^{(s_{1},s_{2})}_{i,j}}\,.$
(2.26)

One can then get an asymptotic $p$-value for this test in the standard way,
via the survival function of the $\chi^{2}$ distribution with $3$ degrees of
freedom.

####  3 Methods

Section 3.1 explains which distributions were simulated, why, and how. Section
3.2 explains how categorical divergences and corresponding approximate
$p$-values were computed. Section 3.3 explains how percentages and percent
changes for droplets with a given number of strains were computed. Section 3.4
explains how I computed the percent changes for gluttonous groups. Section 3.5
explains how I computed the percent changes for picky groups. I plotted
results using Matplotlib [Hun07] version 3.4.1 and Seaborn [Was21] version
0.11.1. The bar plot was made using the microbiome [LS19], phyloseq [MH13],
and ggplot2 [Wic16] packages for the R programming language. Complete
implementation details can be found in the code at the relevant GitLab
repository. See https://gitlab.com/krinsman/droplets.

#####  3.1 Simulated Distributions

To test how much the assumptions of hPoMu may be violated with hPoMu remaining
an adequate working model, I simulated seven distributions. Below, I describe
the two new distributions beyond the five distributions that were already
described in section 3.1. Details of the simulations common to all seven
distributions were given already in section 3.1.4.

######  3.1.1 Simulating High Compositional Heterogeneity Only

This distribution was a hierarchical Exponential hPoDM model (hExhPoDM), with
$\zeta_{C}$ distributed as an $\operatorname{Exponential}(1)$ random variable.

######  3.1.2 Simulating High Compositional and Density Heterogeneities

This distribution was a hierarchical Exponential hNBDM model (hExhNBDM), and
$\zeta$ was distributed as an $\operatorname{Exponential}(1)$ random variable.

#####  3.2 Computing $p$-values and $\chi^{2}$ Divergences

Using the categories and expected counts defined for the “global” goodness of
fit hypothesis test of the null hPoMu using “picky” groups defined in section
2.1.3, I computed Pearson categorical (“$\chi^{2}$”, cf. section 2.2)
divergence values and $\chi^{2}$-approximated $p$-values for all $500$
simulations of all $7$ distributions.

Despite that each simulation corresponded to an extremely large number of
droplets, the expected counts for some categories were very small ($<5$),
making the asymptotic approximation provided by the $\chi^{2}$ assumption
possibly inadequate. To account for this, I generated $10^{9}$ (one billion)
independent replicates of the Multinomial distribution corresponding to the
number of droplets per simulation and the probabilities for each category in
the goodness of fit test. From this, I generated $10^{9}$ replicates from the
sampling distribution under the null of the Pearson categorical divergence
statistic. I then used this Monte Carlo distribution of Pearson categorical
divergence statistics to compute “Monte Carlo $p$-values” for all $500$
simulations from all $7$ distributions. Given the very large number of
replicates, I expect these $p$-values to be more accurate than those from the
$\chi^{2}$ asymptotic approximation. Even so the results using the $\chi^{2}$
approximation were similar, cf. section 5.A.

#####  3.3 Numbers of Droplets with $n$ Strains

Using the relationships:

$\begin{split}\mathbb{E}\left[\left\nonscript\;\middle|\nonscript\;\left\\{d\in[{D}]:\mathbb{S}_{d}(0)=n\right\\}\right\nonscript\;\middle|\nonscript\;\right]=&{D}\cdot\mathbb{P}(\mathbb{S}(0)=n)\\\
=&{D}\cdot\left[\sum_{\mathcal{S}\subseteq[S]:|\mathcal{S}|=n}\mathbb{P}(\operatorname{supp}(\vec{\mathbf{\mathbb{S}}}(0))=\mathcal{S})\right]\,,\end{split}$
(3.1)

I computed via brute force these quantities for $n\in[5]$ for hPoMu:

$\begin{split}&{D}\cdot\left[\sum_{\mathcal{S}\subseteq[S]:|\mathcal{S}|=n}e^{-\sum_{s\not\in\mathcal{S}}f^{(s)}\lambda}\cdot\prod_{s\in\mathcal{S}}\left(1-e^{-f^{(s)}\lambda}\right)\right]\\\
=&{D}\cdot\left[\sum_{\mathcal{S}\subseteq[S]:|\mathcal{S}|=n}e^{-(1-\sum_{s\in\mathcal{S}}f^{(s)})\lambda}\cdot\prod_{s\in\mathcal{S}}\left(1-e^{-f^{(s)}\lambda}\right)\right]\,.\end{split}$
(3.2)

The value for $6+$ strains is just ${D}$ minus the sum of the values for all
$n\in[5]$.

For each of the seven distributions, I stratified and then counted the
droplets from each of the 500 simulations according to how many strains were
present. I then computed the average values for each distribution as the
arithmetic mean over the 500 simulations. (By the law of large numbers the
arithmetic mean is consistent for the true expectation and the only reasonable
choice of measure of “central tendency” in this context.) I then computed the
average percent change for each distribution by subtracting the value expected
under hPoMu from the mean observed value, and then dividing by the expected
value.

#####  3.4 Gluttonous Groups

Using the formulae (2.21) and (2.22), I computed the expected count for each
treatment gluttonous group222222 For the gluttonous “controls”, i.e.
corresponding to the diagonals of the heatmaps, I actually used the sets
${\\{d\in[{D}]:N^{(s)}(0)\geq 1\\}}$, which is a superset of all gluttonous
controls for experiments measuring effects on the growth of strain $s$. The
expected size of these sets is
${{D}\cdot\mathbb{E}[\mathbb{S}^{(s)}(0)]={D}\cdot\mathbb{P}(N^{(s)}(0)\geq
1)}$, which under hPoMu equals ${D}\cdot(1-e^{-f^{(s)}\lambda})$.  under
hPoMu. For each of the seven distributions, for each of the 500 simulations
using NumPy [HMvdW+20] version 1.20.2 I computed via brute force the number of
occurrences of each of the gluttonous groups, and then reported the arithmetic
means over all 500 simulations. I again computed the average percent change as
the observed value minus the value expected under hPoMu divided by the value
expected under hPoMu.

Initially I considered “filtering” differences depending on whether they were
“significant” according to the (asymptotic) $p$-values of the hypothesis tests
defined in section 2.4. However, I found that doing so made the results more
difficult to interpret by obscuring overall trends. The inferred overall
trends were qualitatively the same in either case. Moreover, because each
gluttonous pair corresponds to a distinct hypothesis test, and the gluttonous
groups overlap, this leads to a difficult multiple comparisons problem. One
has to not only decide how to account for multiple testing between simulations
but also “within” simulations. Each of the multiple testing correction schemes
I tried seemed objectionable in some way. Given all of the above, I decided to
report the “raw” average percent change without “filtering” for
“significance”.

#####  3.5 Picky Groups

Using the formulae (2.16) and (2.17), cf. also the equations (2.12) and
(2.13), I computed the expected count for each picky group under hPoMu. Note
that, unlike for gluttonous treatment groups, picky treatment groups do not
overlap. Moreover, unlike for gluttonous control groups, the picky control
group for investigating the effect of strain $s_{1}$ on the rate of growth of
strain $s_{2}$ is the same regardless of what $s_{1}$ is. So the values shown
on the diagonals of the heat maps really correspond to the actual picky
control groups. The procedure was then otherwise exactly the same as for the
gluttonous groups. I did not “filter” according to “significance”.

####  4 Results

Section 4.1 examines the results of the global $\chi^{2}$ goodness of fit
tests computed from the picky contingency tables. Section 4.2 directly
examines both the picky and gluttonous contingency tables.

#####  4.1 Global $\chi^{2}$ Tests

Section 4.1.1 discusses evidence indicating that sampling with replacement is
the most innocuous of the assumptions behind hPoMu. Section 4.1.2 discusses
evidence indicating that ghNBDM distributions with low (but nonzero)
heterogeneity may be adequately modeled by hPoMu. Section 4.1.3 discusses
evidence indicating that ghNBDM distributions with only moderate heterogeneity
can be easily distinguished from hPoMu. Finally, section 4.1.4 discusses
evidence indicating that even further heterogeneity can make differences from
hPoMu extremely obvious.

Figure 5.2:

######  4.1.1 Sampling without Replacement

All available evidence suggests that, even after $15,000,000$ droplets have
been formed, the hTPMH distribution for the chosen value of $K$ is extremely
similar to the corresponding hPoMu distribution.

###### Goodness of Fit $p$-Values Are Approximately Uniformly Distributed

For the hTPMH distribution, figure 5.3 shows how the distribution of (Monte
Carlo) $p$-values over $500$ simulations is approximately uniform, as it would
be if we had actually sampled from the true hPoMu null distribution instead.

Figure 5.3: Monte Carlo $p$-value distribution for sampling without
replacement.

###### Distribution of Pearson Categorical Divergences for Sampling without
Replacement

Moreover, as seen clearly in figure 5.4 (and less so in figure 5.2), the
observed distribution of Pearson categorical divergences from the hTPMH
distribution overlaps strongly with approximately what would have been
expected under the hPoMu null distribution.

Figure 5.4:

######  4.1.2 Very Similar Distributions

All available evidence suggests that the hPoDM with $\zeta_{C}=100$ and hNBDM
with $\zeta=100$ distributions are very similar to the corresponding hPoMu
distribution. Nevertheless, they also appear to differ slightly more than the
hTPMH distribution.

###### Goodness of Fit $p$-Values Are Roughly Uniformly Distributed

For the hPoDM distribution with $\zeta_{C}=100$ and the hNBDM distribution
with $\zeta=100$, figures 5(a) and 5(b) show respectively how their
distributions of (Monte Carlo) $p$-values over $500$ simulations is roughly
similar to the uniform distribution.

(a)

(b)

Figure 5.5: Monte Carlo goodness of fit $p$-value distributions for very
similar distributions.

The $p$-value distribution for hNBDM with $\zeta=100$ in figure 5(b) appears
to fit the uniform distribution worse than the corresponding distribution for
hPoDM with $\zeta_{C}=100$ in figure 5(a). This could possibly be a reflection
of how the hNBDM family incorporates both kind of heterogeneities whereas the
hPoDM family does not. However, it could also possibly be a reflection merely
of the relatively small ($n=500$) sample size.

###### Distributions of Pearson Categorical Divergences for Very Similar
Distributions

Furthermore, figure 5.4 indicates how the distribution of their Pearson
categorical divergence statistics also overlaps well with roughly what would
be expected under the null distribution.

The overlap for hNBDM $\zeta=100$ seen in figure 5.4 does not seem to be
substantially worse (nor better) than that observed for hPoDM $\zeta_{C}=100$,
in contrast to the extent of difference that might have been anticipated after
comparing figures 5(a) and 5(b).

######  4.1.3 Slightly Different Distributions

All available evidence suggests that the hPoDM with $\zeta_{C}=1$ and hNBDM
with $\zeta=1$ distributions are (at least) slightly different from the
corresponding hPoMu distribution.

###### Goodness of Fit $p$-Values Are All Effectively Zero

For the hPoDM distribution with $\zeta_{C}=1$ and the hNBDM distribution with
$\zeta=1$, all of the approximate $p$-values were indistinguishable from $0$
up to numerical precision. Therefore their distributions of approximate
$p$-values, being a “single point spike” or “Dirac delta”, were not plotted.

###### Distributions of Pearson Categorical Divergences for Slightly Different
Distributions

The reason behind these low $p$-values is fairly obvious from either figure
5.2 or 5.6: their distributions of observed Pearson categorical divergences
overlap almost not at all with what would be anticipated under the hPoMu null
distribution. Also, the distribution of Pearson categorical divergences for
hNBDM with $\zeta=1$ is slightly further to the right than that of hPoDM with
$\zeta_{C}=1$, reflecting how the former incorporates both density and
compositional heterogeneity while the latter incorporates only compositional
heterogeneity.

Figure 5.6:

######  4.1.4 Very Different Distributions

All available evidence suggests that the hExhPoDM with
$\mathbb{E}[\zeta_{C}]=1$ and hExhNBDM with $\mathbb{E}[\zeta]=1$
distributions are very different from the corresponding hPoMu distribution.

###### Goodness of Fit $p$-Values Are All Effectively Zero

For the hExhPoDM with $\mathbb{E}[\zeta_{C}]=1$ and hExhNBDM with
$\mathbb{E}[\zeta]=1$ distributions, all of the (approximate) $p$-values were
again indistinguishable from $0$ up to numerical precision, as was also the
case for hPoDM $\zeta_{C}=1$ and hNBDM $\zeta=1$. The approximate $p$-values
were again not plotted for the same reasons as for hPoDM $\zeta_{C}=1$ and
hNBDM $\zeta=1$.

###### Distributions of Pearson Categorical Divergences for Very Different
Distributions

The reason for these small $p$-values is also the same as for hPoDM
$\zeta_{C}=1$ and hNBDM $\zeta=1$: as can be seen clearly in either figure 5.2
or figure 5.7, the peaks of their distributions of Pearson categorical
divergences completely fail to overlap with the range of values that would be
anticipated under the hPoMu null distribution. Both peaks are also much
further to the right than anything observed for any of the other five
distributions, indicating that the hExhPoDM $\mathbb{E}[\zeta_{C}]=1$ and
hExhNBDM $\mathbb{E}[\zeta]=1$ distributions are (by far) the most
heterogeneous of the seven distributions simulated.

Figure 5.7:

Analogous to what occurs for hPoDM $\zeta_{C}=1$ and hNBDM $\zeta=1$, the peak
for hExhNBDM $\mathbb{E}[\zeta]=1$ is further to the right than the peak for
hExhPoDM $\mathbb{E}[\zeta_{C}]=1$, reflecting how the former incorporates
both kind of heterogeneities whereas the latter does not. (Thus hExhNBDM
$\mathbb{E}[\zeta]=1$ can be considered the single most heterogeneous
distribution of the seven that were simulated.) However, in contrast to the
situation observed for hPoDM $\zeta_{C}=1$ and hNBDM $\zeta=1$, the peaks for
hExhNBDM $\mathbb{E}[\zeta]=1$ and hExhPoDM $\mathbb{E}[\zeta_{C}]=1$ do not
overlap at all, with the peak for hExhNBDM $\mathbb{E}[\zeta]=1$ being greatly
further to the right than that for hExhPoDM $\mathbb{E}[\zeta_{C}]=1$.

#####  4.2 Direct Examination of Contingency Tables

Section 4.2.1 explains how the observed effects on data throughput support the
“taxonomy” of distributions previously proposed in section 4. Section 4.2.2
explains the evidence indicating that increased heterogeneity decreases the
average number of droplets with 3 or more strains. Section 4.2.3 explains the
effect of increased heterogeneity on data throughput for droplets with 2 or
fewer strains, with section 4.2.3 explaining what happens for compositional
heterogeneity without density heterogeneity, and section 4.2.3 explaining what
happens when both types of heterogeneity are present.

(a)

(b)

Figure 5.8: Droplets stratified by numbers of strains found inside droplets.

######  4.2.1 Distributions Again Fall Into Same Groups

From figures 8(a), 8(b), 9(a), 10(a), 10(b), 11(a), 12(a), and 12(b) we see
that hTPMH, hPoDM with $\zeta_{C}=100$, and hNBDM with $\zeta_{C}=100$ are
again (cf. section 4) either indistinguishable from or very similar to hPoMu.
Any differences in average composition are only visible for medium and high
heterogeneities, and easily visible usually only for high heterogeneity.

(a) Figure 5.9: Gluttonous groups for low heterogeneity distributions.

(a)

(b)

Figure 5.10: Gluttonous groups for low heterogeneity distributions. (a) Figure
5.11: Picky groups for low heterogeneity distributions.

(a)

(b)

Figure 5.12: Picky groups for low heterogeneity distributions.

######  4.2.2 Increased Heterogeneity Decreases Number of Droplets with 3 or
More Strains

From figures 8(a) and 8(b) we see that, both for the hPoDM and hNBDM families,
as the compositional heterogeneity increases, the number of droplets with 3 or
more strains decreases in both absolute and relative terms. Since droplets
with 3 or more strains are a component of the gluttonous groups, their
decrease with increased compositional heterogeneity is also reflected in
figures 13(a), 14(a), 13(b), and 14(b). We can also see from figures 8(b),
13(a), 14(a), 13(b), and 14(b) that this decrease in the number of droplets
with 3 or more strains is more pronounced in the absence of density
heterogeneity.

(a)

(b)

Figure 5.13: Gluttonous groups for high heterogeneity distributions.

(a)

(b)

Figure 5.14: Gluttonous groups for high heterogeneity distributions.

######  4.2.3 Effect of Increased Heterogeneity on Droplets with Fewer than 3
Strains

For droplets with fewer than 3 strains, the effects that occur as
compositional heterogeneity increases depend on whether density heterogeneity
is present.

###### Increased Compositional Heterogeneity Only

In the absence of density heterogeneity, we see an increase in the number of
droplets with one strain, with a smaller increase in the number of droplets
with two strains. This is best evidenced in figure 8(b). This slight increase
in the number of droplets with two strains causes the slight increase in the
number of picky treatments shown in figures 15(a) and 15(b). As is clear also
from figure 8(a), there is little change in the number of empty droplets with
zero strains regardless of the level of compositional heterogeneity.

(a)

(b)

Figure 5.15: Picky groups for high heterogeneity distributions.

(a)

(b)

Figure 5.16: Picky groups for high heterogeneity distributions.

###### Increased Compositional and Density Heterogeneity

In the presence of density heterogeneity, as seen in figure 8(b) the number of
droplets with a single strain does still increase slightly. However, in
contrast to what happens in the absence of density heterogeneity, in the
presence of density heterogeneity there is an observable decrease in the
number of droplets with two strains. While this trend is less obvious based
only on figure 8(b), it is clear based on the observed decrease in the number
of picky treatments visible in figures 16(a) and 16(b). Nevertheless it is
also clear from figure 8(b) and from a comparison of figures 16(b) and 14(b)
(and 16(a) and 14(a)) that the decrease in the number of droplets with two
strains is not as pronounced (at least in relative terms) as the decrease in
the number of droplets with three or more strains. Also in contrast to what
happens in the absence of density heterogeneity, in the presence of density
heterogeneity figures 8(a) and 8(b) show a marked increase in the number of
empty droplets with zero strains as the compositional heterogeneity increases.

####  5 Discussion

Section 5.1 discusses the taxonomy of four groups of distributions identified
in section 4. Section 5.2 discusses the distributions which are most likely to
exhibit differences from hPoMu that are important in practice. Section 5.3
discusses how the interpretation of the results from section 4 depends
fundamentally on our opinion of which droplets constitute “good data”. Section
6 discusses how the results from section 4 indicate what the next step is in
order to recover true interaction networks from MOREI data.

#####  5.1 Overview of Groups of Distributions

The seven simulated distributions can be grouped into roughly four levels of
assumption violation.

For the first two levels, the first consisting of hTPMH and the second
consisting of hPoDM $\zeta_{C}=100$ and hNBDM $\zeta=100$, hPoMu is most
likely an adequate working model. Cf. again section 5.2 from the previous
chapter.

The third level consists of hPoDM $\zeta_{C}=1$ and hNBDM $\zeta=1$, which
produce data that can be relatively easily distinguished from data produced by
hPoMu. Cf. again section 5.3.

In all cases, the results using the contingency tables for the target
parameters and resulting $\chi^{2}$-divergences do nothing to change the
qualitative conclusions that were made in the previous chapter. Similarity (or
lack thereof) of these working models, as measured by log-likelihood ratios,
corresponds to similarity (or lack thereof) of their contingency tables of
target parameters.

The fourth level consists of hExhPoDM $\mathbb{E}[\zeta_{C}]=1$ and hExhNBDM
$\mathbb{E}[\zeta]=1$, which produce data so substantially different from that
produced by hPoMu that it is almost impossible to not distinguish. It is very
difficult to believe that hPoMu could possibly serve as an adequate working
model for such data.

#####  5.2 Very Different Distributions

The simulated hExhPoDM and hExhNBDM distributions technically do not belong to
the ghNBDM family of distributions, since their concentration “parameters” do
not have fixed values, and instead are exponentially distributed. This
distinction will be particularly relevant in sections 4.2.3 and 4.3.3.

Moreover, although the mean of these exponential distributions is $1$, the
behavior of the simulated hExhPoDM with $\mathbb{E}[\zeta_{C}]=1$ and hExhNBDM
with $\mathbb{E}[\zeta]=1$ distributions is substantially different from that
of the hPoDM $\zeta_{C}=1$ and hNBDM $\zeta=1$ distributions. This is most
likely because concentration parameter values sampled from the
$\operatorname{Exponential}(1)$ are $1-\frac{1}{e}\approx 63.2\%$ of the time
smaller than $1$. (The median of the distribution is $\ln(2)\approx.693<1$.)

Thus the simulated hExhPoDM with $\mathbb{E}[\zeta_{C}]=1$ and hExhNBDM with
$\mathbb{E}[\zeta]=1$ can be assumed a priori to resemble more closely hPoDM
and hNBDM distributions with values of $\zeta_{C}$ and $\zeta$ (substantially)
smaller than $1$. This corresponds to what we see in practice, because
heterogeneities, as measured by over-dispersion, are inversely proportional to
the concentration parameters. So for example $\zeta_{C}=\ln(2)$ would
correspond to an hPoDM distribution with $\approx 44\%$ more heterogeneity
than the hPoDM $\zeta_{C}=1$ distribution.

Considering that the simulated hExhPoDM with $\mathbb{E}[\zeta_{C}]=1$ and
hExhNBDM with $\mathbb{E}[\zeta]=1$ correspond $\approx 63.2\%$ to
distributions with heterogeneities no less than those of hPoDM $\zeta_{C}=1$
and hNBDM $\zeta=1$, and $1-\frac{1}{\sqrt{e}}\approx 39.3\%$ of the time to
distributions with heterogeneities no less than twice those of hPoDM
$\zeta_{C}=1$ and hNBDM $\zeta=1$, it is unsurprising that the
aggregate/average heterogeneities observed for these distributions is much
larger than those observed for either the hPoDM $\zeta_{C}=1$ or hNBDM
$\zeta=1$ distributions.

Therefore fact that data generated by these two distributions is very unlikely
to be adequately modeled by hPoMu further suggests that data generated by
hPoDM or hNBDM distributions with $\zeta_{C}$ or $\zeta$ values not much
smaller than $1$ will also be inadequately modeled by hPoMu. It is a priori
unclear which, if any, values of the concentration parameters smaller than $1$
are realistic in practice.

#####  5.3 The Effect of Heterogeneity on Data Throughput is Equivocal

The general effect that heterogeneities will have on the data throughput of
MOREI is not straightforward. The answer depends on the particular combination
of density heterogeneity and compositional heterogeneity present. Given both
density and compositional heterogeneity, it is clear that the effect on data
throughput is negative. Both the numbers of picky treatments and gluttonous
treatments decrease.

In the presence of only compositional heterogeneity but no density
heterogeneity, whether the effect on data throughput is positive depends
primarily on our opinion of what is useful data. If our primary concern is
avoiding “data starvation” and being able to say anything about the
interactions of the rarest strains, then the strong decrease in the number of
gluttonous treatment groups clearly means that the effect of heterogeneity is
negative. On the other hand, if we consider only the data from the picky
treatments useful, then in the absence of density heterogeneity the effect of
compositional heterogeneity is slightly positive.

####  6 Conclusion

###### Findings and Contributions

I confirmed that more severe failures of the assumptions from section 3.3
cause more severe discrepancies with the predictions derived from hPoMu. The
nature of the effect depends on the chosen grouping (picky or gluttonous) of
droplets defining the targeted estimands. This demonstrates that failures of
the homogeneity assumptions of hPoMu potentially could be severe enough to
require attention in practice.

###### Practical Implications

The practical implications of this chapter remain limited by ignorance of
which values of density and compositional concentration parameters are
realistic in practice. If realistic values correspond to sufficiently small
violations of the homogeneity assumptions, then this chapter actually provides
evidence that such violations could be ignored in practice. Likewise, if
realistic values correspond to sufficiently large violations of the
homogeneity assumptions, then this chapter provides evidence of the opposite.

###### Next Steps and Open Questions

What is also clear is that, if we wish to be able to predict generally what
effects heterogeneities may have on the data throughput of MOREI in practice,
we must be able to estimate both density heterogeneity and compositional
heterogeneity from empirical data. This is the topic of chapter 6. Given such
estimates from real-world data, we can get some sense for the amount of either
type of heterogeneity one is likely to encounter in practice. That makes it
possible to create more accurate simulations. More accurate simulations in
turn better assess which inference methods are most capable of recovering true
interaction networks from MOREI data.

### Chapter 5. Appendices

###  Appendix 5.A $p$-Values from $\chi^{2}$ Approximation

Even for 15 million droplets, some categories corresponding to the rarest
strains had very small expected counts ($<5$). Thus the $\chi^{2}$
approximation for the sampling distribution of the categorical divergence
under the null distribution was unlikely to be accurate. In contrast, the
Monte Carlo approximation to the null distribution is from $10^{9}$ (one
billion) replicates. Therefore the $p$-values from the Monte Carlo
approximation are more likely to be representative of the true $p$-values. For
completeness, here I report the $p$-value distributions computed from the
$\chi^{2}$ approximation. While similar to the $p$-values from the Monte Carlo
approximation, they appear to over-estimate the differences from the null.

(a)

(b)

(c)

Figure 5.A.1:  Approximate goodness of fit $p$-value distributions using
$\chi^{2}$ approximation.

###  Appendix 5.B Additional Pearson Categorical Divergence Distributions

Here are two additional plots, one focusing on those distributions with
compositional heterogeneity only, and the other on those distributions with
both heterogeneities.

(a)

(b)

###  Appendix 5.C Multiple Representatives Problem

Another potential challenge arising from the randomness of the initial droplet
formation process is the “multiple representatives problem”, where within a
given droplet there is initially a strain that is represented by multiple
cells. Cf. figure 5.C.1. Even if we were more concerned about the effects
randomness of the initial droplet formation process causes via the “multiple
representatives problem” than the effects it has on data throughput, I argue
we should still first focus on the effects it has on data throughput, for two
reasons.