Re-examining aggregation in the Tallis-Leyton model of parasite acquisition

Tallis and Leyton (1969) proposed the following model for the parasite burden $M(a)$ of a definitive host at age $a$, conditional on survival of the host to age $a$. The host is assumed to be parasite free at birth so $M(0)=0$. The host makes infective contacts during its lifetime following a Poisson process with rate parameter $\lambda$. At each infective contact, the number of parasites that enter the host is a random variable $N$ and once a parasite enters the host it survives for a random period $T$. The parasites have no effect on the host mortality. The lifetimes of parasites, numbers of parasites entering the host at infective contacts, and the Poisson process of infective contacts are all assumed to be independent. Although we won’t make use of this fact, we note that this process also describes an infinite server queue with bulk arrivals and general independent service times (Holman et al., 1983).

We let $G_{X}$, $F_{X}$ and $\bar{F}_{X}=1-F_{X}$ denote the probability generating function (PGF), distribution function, and survival function of a random variable $X$. We write $G_{M}(\cdot;a)$ for the PGF of $M(a)$, conditional on survival of the host to age $a$. Tallis and Leyton (1969) showed

$$G_{M}(z;a)=\exp\left(\lambda\int_{0}^{a}\left[G_{N}(1+\bar{F}_{T}(a-s)(z-1))-1\right]ds\right).$$

Standard arguments show that the mean and variance of $M(a)$ are

	$\displaystyle\mu(a)$	$\displaystyle=\lambda\mathbb{E}N\int_{0}^{a}\bar{F}_{T}(s)\,ds$
	$\displaystyle\sigma^{2}(a)$	$\displaystyle=\lambda\mathbb{E}N(N-1)\int_{0}^{a}\bar{F}_{T}^{2}(s)\,ds+\lambda\mathbb{E}N\int_{0}^{a}\bar{F}_{T}(s)\,ds.$

Hence, the variance-to-mean ratio is

$$\text{VMR}(M(a))=1+\frac{\mathbb{E}N(N-1)}{\mathbb{E}N}\frac{\int^{a}_{0}\bar{F}^{2}_{T}(s)ds}{\int^{a}_{0}\bar{F}_{T}(s)ds}.$$

(1)

Assuming $\mathbb{E}N<\infty$ and $\mathbb{E}T<\infty$, the limiting distribution of parasite burden as $a\to\infty$ exists and has PGF

$$G_{M}(z,\infty)=\exp\left(\lambda\int_{0}^{\infty}\left[G_{N}(1+\bar{F}_{T}(s)(z-1))-1\right]ds\right).$$

Since scaling of the host age, the lifetime distribution of parasites, and the inverse of the rate of infective contacts by a common factor results in the same distribution for the host’s parasite burden, we assume that the expected lifetime of a parasite is 1.

Lorenz (1905) proposed the Lorenz curve as a graphical measure of inequality. The following general definition of the Lorenz curve was given by Gastwirth (1971).

Adapting the description in Arnold and Sarabia (2018, Section 3.1) to a parasitology context, the Lorenz curve $L(u)$ represents the proportion of the parasite population infecting the least infected $u$ proportion of the host population. When all hosts are infected with the same number of parasite, the Lorenz curve is given by $L(u)=u$ and is called the egalitarian line. Note that $L_{X}(u)\geq u$ for all $u\in[0,1]$.

The Lorenz curve defines a partial order on the class of all distributions on $[0,\infty)$ with finite mean (Arnold and Sarabia, 2018, Definition 3.2.1).

The Lorenz curves of some standard distributions are given in Arnold and Sarabia (2018, Section 6.1). The negative binomial distribution is extensively used in parasitology. When negative binomial distributions are parameterised in terms of the mean $m$ and $k$, they can be compared in the Lorenz order (McVinish and Lester, 2024). Specifically, let $\mathsf{NB}(m,k)$ denote the negative binomial distribution with PGF

$$G(z)=\left(\frac{k}{k+m-mz}\right)^{k}.$$

Then

• (i)

  for any $k>0$ and $0<m_{1}<m_{2}$, $\mathsf{NB}(k,m_{2})\leq_{\rm Lorenz}\mathsf{NB}(k,m_{1})$, and

• (ii)

  for any $m>0$ and $0<k_{1}<k_{2}$, $\mathsf{NB}(k_{2},m)\leq_{\rm Lorenz}\mathsf{NB}(k_{1},m)$.

Closely related to the Lorenz order is the convex order of distributions.

These two orderings are related since $X\leq_{\rm Lorenz}Y$ if and only if

$$\mathbb{E}\,\phi\left(\frac{X}{\mathbb{E}X}\right)\leq\mathbb{E}\,\phi\left(\frac{Y}{\mathbb{E}Y}\right)$$

for every continuous convex function $\phi$ (Arnold and Sarabia, 2018, Corollary 3.2.1). In other words,

$\displaystyle X$

$\displaystyle\leq_{\rm Lorenz}Y$

is equivalent to

$\displaystyle\frac{X}{\mathbb{E}X}$

$\displaystyle\leq_{\rm cx}\frac{Y}{\mathbb{E}Y}.$

Shaked and Shanthikumar (2007, Section 3.A) provide an extensive review of results on the convex order. We briefly mention some of the important results that are used in our analysis.

• •

  The convex order is closed under weak limits provided the expectations also converge (Shaked and Shanthikumar, 2007, Theorem 3.A.12 (c)).

• •

  The convex order is closed under mixtures (Shaked and Shanthikumar, 2007, Theorem 3.A.12 (b)). Let $X$, $Y$, and $\Theta$ be random variables and write $[X\mid\Theta=\theta]$ and $[Y\mid\Theta=\theta]$ for the conditional distributions of $X$ and $Y$ given $\Theta=\theta$. If $[X\mid\Theta=\theta]\leq_{\rm cx}[Y\mid\Theta=\theta]$ for all $\theta$ in the support of $\Theta$, then $X\leq_{\rm cx}Y$. As an application of this property we can say that if $X\leq_{\rm cx}Y$ and $Z$ is an independent non-negative random variable, then $ZX\leq_{\rm cx}ZY$.

• •

  The convex order is closed under convolutions (Shaked and Shanthikumar, 2007, Theorem 3.A.12 (d)). Let $X_{1},X_{2},\ldots,X_{k}$ and $Y_{1},Y_{2},\ldots,Y_{k}$ be two sets of independent random variables . If $X_{j}\leq_{\rm cx}Y_{j}$ for $j=1,2,\ldots,k$, then

  $$\sum_{j=1}^{k}X_{j}\leq_{\rm cx}\sum_{j=1}^{k}Y_{j}.$$

• •

  Combining the properties of closure under mixtures and closure under convolutions, we see the convex order is closed under random sums so

  $$\sum_{j=1}^{K}X_{j}\leq_{\rm cx}\sum_{j=1}^{K}Y_{j},$$

  for any non-negative integer random variable $K$. As an application of the closure under random sums property of the convex order, consider two random variables $K$ and $\tilde{K}$ that related by binomial thinning. That is, $G_{\tilde{K}}(z)=G_{K}(1-p+pz)$ for some $p\in(0,1)$. Then $K\leq_{\rm Lorenz}\tilde{K}$ (McVinish and Lester, 2020, Section 3)

• •

  The closure under random sums property can be adapted to the case where the $X_{1},X_{2},\ldots$ and $Y_{1},Y_{2},\ldots$ are two iid sequences with $X\leq_{\rm cx}Y$, and $K_{1}$ and $K_{2}$ are non-negative integer random variables such that $K_{1}\leq_{\rm cx}K_{2}$. In this case, Shaked and Shanthikumar (2007, Theorem 3.A.13) implies

  $$\sum_{j=1}^{K_{1}}X_{j}\leq_{\rm cx}\sum_{j=1}^{K_{2}}Y_{j}.$$

• •

  The survival function can be used to establish if two random variables can be compared in the convex order. If $X$ and $Y$ are two random variables with the same mean, then $X\leq_{\rm cx}Y$ if $\bar{F}_{X}-\bar{F}_{Y}$ has a single sign change and the sign sequence is $+,-$ (Shaked and Shanthikumar, 2007, Theorem 3.A.44(b)). This property can also be used to characterise the convex order (Shaked and Shanthikumar, 2007, Theorem 3.A.45).

In practice, levels of aggregation are compared with numerical summaries rather than using the entire Lorenz curve. A useful measure of aggregation respects the Lorenz ordering. That is, if $I(\cdot)$ is a measure of aggregation that respects the Lorenz ordering and $X\leq_{\rm Lorenz}Y$, then $I(X)\leq I(Y)$. Arnold and Sarabia (2018, Chapter 5) review several inequality measures and these can be applied as measures of aggregation. We restrict our attention in this paper to the Gini index, the Hoover index (also known as the Pietra index, or the Robin-Hood index), $1-\text{prevalence}$, and the coefficient of variation.

The Gini index (Gini, 1924) is given by twice the area between the egalitarian line and the Lorenz curve. For a random variable $X$, the Gini index can be expressed as

$$G(X)=\frac{\mathbb{E}\,|X-\tilde{X}|}{2\,\mathbb{E}X},$$

where $\tilde{X}$ is an independent random variable with $\tilde{X}\stackrel{{\scriptstyle d}}{{=}}X$. The Hoover index is given by the maximum vertical distance between the egalitarian line and the Lorenz curve. McVinish and Lester (2020) argue that this index could be useful due to its simple interpretation as the proportion of the parasite population that would need to be redistributed among the hosts in order for all hosts to have the same parasite burden. The Hoover index can be expressed as

$$H(X)=\frac{\mathbb{E}\left|X-\mathbb{E}X\right|}{2\,\mathbb{E}X}$$

Prevalence, the probability that a host is infected by at least one parasite, is an important quantity in parasitology. Although prevalence is not usually thought of as a measure of aggregation, we may express $1-\text{prevalence}$ in terms of the Lorenz curve as

$$1-\text{prevalence}=\max\{u:L(u)=0\}.$$

For the Tallis-Leyton model,

$$1-\text{prevalence}=G_{M}(0;a)=\exp\left(\lambda\int_{0}^{a}\left[G_{N}(1-\bar{F}_{T}(a-s))-1\right]ds\right).$$

(2)

Finally, the coefficient of variation is given by

$$CV(X)=\frac{\sqrt{\text{Var}(X)}}{\mathbb{E}X}.$$

This measure is rarely used in parasitology, though it is mentioned in some reviews on parasite aggregation such as Wilson et al. (2001) and McVinish and Lester (2020). As means and variances are commonly reported in empirical studies and are often easily calculated for theoretical models, it may be useful in some contexts. For example, the squared coefficient of variation for the Tallis-Leyton model is

$$\text{CV}^{2}(M(a))=\frac{1}{\lambda}\frac{\mathbb{E}N(N-1)}{(\mathbb{E}N)^{2}}\frac{\int^{a}_{0}\bar{F}^{2}_{T}(s)ds}{\left(\int^{a}_{0}\bar{F}_{T}(s)ds\right)^{2}}+\frac{1}{\lambda\mathbb{E}N\int^{a}_{0}\bar{F}_{T}(s)ds}.$$

(3)

These indices are related by the following inequality

$$1-\text{prevalence}\leq H(X)\leq G(X)\leq H(X)(2-H(X))\leq CV(X)$$

(Taguchi, 1968, McVinish and Lester, 2020). In particular, if $\mathbb{E}X\leq 1$, then we have the equality $1-\text{prevalence}=H(X)$. The Gini index and Hoover index can be further related to the coefficient of variation when the distribution of parasites is approximately normal. Suppose $X_{1},X_{2},\ldots$ is a sequence of random variables such that

$$\frac{X_{n}-\mathbb{E}X_{n}}{\sqrt{\text{Var}(X_{n})}}=:Z_{n}\stackrel{{\scriptstyle d}}{{\to}}Z,$$

where $Z\sim\mathsf{N}(0,1)$. As $X_{n}\geq 0$ with probability one, the above limit is only possible if $CV(X_{n})\to 0$. Nevertheless, the ratio of the Hoover index to the coefficient of variation still has a well defined limit. The Hoover index of $X_{n}$ can be expressed as

$$H(X_{n})=\frac{\mathbb{E}\left|X_{n}-\mathbb{E}X_{n}\right|}{2\mathbb{E}X_{n}}=\frac{\sqrt{\text{Var}(X_{n})}}{2\mathbb{E}X_{n}}\mathbb{E}|Z_{n}|.$$

Since $\mathbb{E}Z_{n}^{2}=1$, this collection of random variables is uniformly integrable and $\mathbb{E}|Z_{n}|\to\mathbb{E}|Z|=\sqrt{2/\pi}$. Hence,

$$\frac{H(X_{n})}{CV(X_{n})}\to\frac{1}{\sqrt{2\pi}}.$$

(4)

Similarly, the Gini index of $X_{n}$ can be expressed as

$$G(X_{n})=\frac{\mathbb{E}|X_{n}-\tilde{X}_{n}|}{2\mathbb{E}X_{n}}=\frac{\sqrt{\text{Var}(X_{n})}}{2\mathbb{E}X_{n}}\mathbb{E}|Z_{n}-\tilde{Z}_{n}|,$$

where $\tilde{Z}_{n}$ is an independent random variable with $\tilde{Z}_{n}\stackrel{{\scriptstyle d}}{{=}}Z_{n}$. Applying the asymptotic normality and uniform integrability of the $Z_{n}$,

$$\frac{G(X_{n})}{CV(X_{n})}\to\frac{1}{\sqrt{\pi}}.$$

(5)

From equation (3), the coefficient of variation can be relatively easily evaluated for the Tallis-Leyton model. Numerical integration of $\bar{F}_{T}(s)$ and $\bar{F}_{T}^{2}(s)$ may be required, but the dependence on age and $\lambda$ is explicit. Similarly, $1-\text{prevalence}$ could be evaluated with a single numerical integration using (2). On the other hand, evaluation of the Hoover and Gini indices require evaluation of the probability mass function. In the examples of the next section, we numerically evaluate the probability mass function of $M(a)$ by inverting $G_{M}(z;a)$ using the Abate-Whitte algorithm (Abate and Whitt, 1992). The algorithm was implemented in MATLAB (The MathWorks Inc., 2022a) using the vpa function in the Symbolic Math Toolbox (The MathWorks Inc., 2022b) for high precision arithmetic.

Let $U_{1},U_{2},\ldots$ be a sequence of independent standard uniform random variables and define $X:\mathbb{Z}_{\geq 0}\times[0,1]\to\mathbb{Z}_{\geq 0}$ by

$$X(n,v):=\sum_{k=1}^{n}\mathbb{I}(U_{k}\leq v),$$

with $X(n,v)=0$ when $n=0$. For given $n$ and $v$ the distribution of $X(n,v)$ is $\mathsf{Bin}(n,v)$.

Our first comparison result concerns the effect of the rate of infective contacts on the distribution of the host’s parasite burden. The rate of infective contacts has no effect on the variance-to-mean ratio (1), whereas the coefficient of variation is strictly decreasing as the rate of infective contacts increases (3). The following result shows that any index respecting the Lorenz order is decreasing as a function of the rate of infective contacts.

Figure 1 shows the four indices (Gini, Hoover, $1-\text{prevalence}$, and coefficient of variation) for a host aged 3 where $\lambda\in[0.25,128]$, $N\sim\mathsf{NB}(1,1)$, and $T\sim\mathsf{Exp}(1)$. Since the coefficient of variation (3) is proportional to $\lambda^{-1/2}$, it is not displayed for small values of $\lambda$. We see that all four indices are strictly decreasing as a function of $\lambda$. For $\lambda\leq 1$, $\mu(3)\leq 1$ so the Hoover index and $1-\text{prevalence}$. The Hoover index appears to display some discontinuity in the first derivative at points where the expectation is integer valued. This behaviour is less apparent at larger values of $\lambda$.

We now consider the role of the distribution of the number of parasites that enter the host during an infective contact. As a concrete example, suppose $N\sim\mathsf{NB}(m,k)$. Then the variance-to-mean ratio is

$$\text{VMR}(M(a))=1+cm(1+1/k),$$

where $c>0$ is a constant depending on $F_{T}$. From this expression we see that the variance-to-mean ratio is increasing in $m$ but decreasing in $k$. In contrast, the coefficient of variation of $M(a)$ is decreasing in both $m$ and $k$. The next two results show that the distribution of the host’s parasite burden is decreasing in the Lorenz order as functions of both $m$ and $k$. The first of these results requires the distributions being compared to have the same expectation.

For distributions with different means, we consider only the case where $N$ and $\tilde{N}$ are related by binomial thinning. Recall that if $G_{\tilde{N}}(z)=G_{N}(1-p+pz)$ for some $p\in(0,1)$, then $\tilde{N}\leq_{\rm Lorenz}N$.

As noted previously, $\mathsf{NB}(k,m)\leq_{\rm Lorenz}\mathsf{NB}(\tilde{k},\tilde{m})$ if $\tilde{k}\leq k$ and $\tilde{m}\leq m$. If $N\sim\mathsf{NB}(k,m)$, $\tilde{N}\sim\mathsf{NB}(\tilde{k},\tilde{m})$ and all other model parameters are equal, then Theorems 3 and 4 together imply that $M(a)\leq_{\rm Lorenz}\tilde{M}(a)$. Figure 2 shows the Hoover and Gini indices as functions of the negative binomial parameters $m$ and $k$ for a parasite host system with host aged 10, rate of infective contacts 5, and parasite lifetimes following an exponential distribution with mean 1. Both indices are decreasing in both $m$ and $k$ as we expect from the above results. The contours of the Hoover index display some discontinuity in the first derivative for $m=1/5$ ($\log_{2}(m)\approx-2.3$), which corresponds to a mean of 1 for the host. The contours for both the Hoover and Gini indices tend to become parallel to the respective axes as $m\to\infty$ and $k\to\infty$. This is a consequence of the limiting behaviour of the negative binomial distribution (Adell and Cal, 1994).

It is natural to consider which distribution for $N$ results in the least aggregated distribution for the host’s parasite burden. This requires determining the smallest distribution in the convex ordering. The convex order requires that the distributions compared have the same expected value so let $n=\mathbb{E}N$. Define the random variable $N^{\prime}$ such that

$\displaystyle\mathbb{P}(N^{\prime}=\lfloor n\rfloor)$

$\displaystyle=\lfloor n\rfloor+1-n$

$\displaystyle\mathbb{P}(N^{\prime}=\lfloor n\rfloor+1)$

$\displaystyle=n-\lfloor n\rfloor.$

In the supplementary material of McVinish and Lester (2020) it was shown that $N^{\prime}\leq_{\rm cx}N$ so $N^{\prime}$ is smallest distribution in convex order with expectation $n$. When $n\leq 1$, the smallest distribution in convex order for $N$ leads to $M(a)$ having a Poisson distribution. There is no largest distribution in the convex order.

Differentiating (1) with respect to age shows the variance-to-mean ratio is a decreasing function of age. Since the expected parasite burden is increasing in age, the coefficient of variation is also decreasing in age. The following result shows the host’s parasite burden is decreasing in the Lorenz order as a function of age.

The proof is built from the following lemmas.

Figure 3 shows the four indices (Gini, Hoover, $1-\text{prevalence}$, and coefficient of variation) for a host in the Tallis-Leyton model with $\lambda=5$, $N\sim\mathsf{NB}(1,1)$, and $T\sim\mathsf{Exp}(1)$. All four indices are strictly decreasing in host age. As in Figure 1, the Hoover index coincides with $1-\text{prevalence}$ for $\mu(a)\leq 1$ and displays some discontinuity in the first derivative at points where the expectation is integer valued.

Our final comparison result concerns the parasite lifetime distribution. The result shows that variability in the parasite lifetimes decreases the variability in the host’s parasite burden. In particular, the result implies that the host’s parasite burden is most aggregated when parasites have constant lifetimes.

As noted previously, for the host’s parasite burden to converge to a normal distribution, the coefficient of variation must tend to 0. In the Tallis-Leyton model, this is only possible when the rate of infective contacts tends to infinity.