A graphical exploration of the relationship between parasite aggregation indices
R. McVinish R. J. G. Lester


]A graphical exploration of the relationship between parasite aggregation indices 

]R. McVinish
School of Mathematics and Physics, University of Queensland
r.mcvinish@uq.edu.au
]R.J.G. Lester
School of Biological Sciences, University of Queensland



§ ABSTRACT 
    The level of aggregation in parasite populations is frequently incorporated into ecological studies.  It is measured in various ways including variance-to-mean ratio, mean crowding, the k parameter of the negative binomial distribution and indices based on the Lorenz curve such as the Gini index (Poulin's D) and the Hoover index. Assuming the frequency distributions follow a negative binomial, we use contour plots to clarify the relationships between aggregation indices, mean abundance and prevalence. The contour plots highlight the nonlinear nature of the relationships between these measures and suggest that correlations are not a suitable summary of these relationships.



[
    [
    September 5, 2024
=====================






§ INTRODUCTION


Investigations into parasite population dynamics frequently require an indicator of the level of aggregation in the parasite population <cit.>. As the concept of aggregation in parasites is poorly defined <cit.>, aggregation has been measured in various ways. Commonly used indices include prevalence, the Variance-to-Mean Ratio (VMR), and the k parameter of the negative binomial distribution. Closely related to VMR and k are mean crowding and patchiness <cit.> which can be seen as more direct measures of the competitive experience of parasites within a host <cit.>. Two other indices are derived from the Lorenz curve <cit.>, the most widely accepted quantification of inequality. <cit.> proposed using the Gini index <cit.>, which has since become widely used in parasitology <cit.>.  The Hoover index (aka Pietra index) has more recently been proposed to measure parasite aggregation <cit.>.  

This paper clarifies and extends our previous work on aggregation. It was stimulated by a recent paper by <cit.> which correlated aggregation indices with mean abundance and prevalence using simulated data.  We present a more accurate representation using ‘contour plots’, calculated directly from the parameters of the negative binomial distributions.  The plots provide a simple and more insightful way to comprehend the relationships.





§ CONTOUR PLOTS


The contours show combinations of two indices, specified on the vertical and horizontal axes, that give rise to similar values of the third index. Contour plots, developed in the 16th century <cit.>, are widely used in other disciplines but rarely in parasitology (e.g. <cit.>).

Our analysis assumed that parasite burden is adequately modelled by a negative binomial distribution <cit.>. Following the typical practice in parasitology, we parameterised the negative binomial distribution in terms of mean abundance, m, and the parameter k which controls the shape of the distribution. We did not make any assumption on the distribution of m and k. We used the range of values for m and k suggested by the extensive data of <cit.>. Their values for m, k and prevalence are superimposed on several of the contour plots as dot points.  

To construct a contour plot of an aggregation index against m and k, we expressed the aggregation index as a function of m and k. The population values of several indices can be expressed simply in terms of m and k: 
    prevalence = 1 - (k/(k+m))^k,
 VMR = 1 + m/k , mean crowding = m +m/k, and patchiness = 1 + 1/k. The Gini and Hoover indices lack simple expressions in terms of m and k however, they can still be evaluated numerically. The Hoover index can be expressed in terms of m and k by applying <cit.>,

    H=F(m;k,m)-F(m-1;k+1,m+m/k),

where F(x; m, k) is the cumulative distribution function of the negative binomial distribution with k and mean m evaluated at x.  Further details are given in the Appendix. The cumulative distribution function of the negative binomial distribution, F, is available in statistical packages such as R <cit.>. The Gini index can be expressed as

    G= (1+m/k)  _2F_1 (k+1,1/2,2;-4 m/k(1+m/k) ),

where _2F_1 is the Gaussian hypergeometric function <cit.>. This can be evaluated in R using the hypergeo package <cit.>. Calculating indices directly from the parameters of the negative binomial distribution rather than using simulated data obviates the need to consider the uncertainty of estimates and the effects of different sample sizes. 

We also employed contour plots to examine the relationship between aggregation indices, k, and prevalence. This required first solving the equation 
    prevalence = 1 – (k/(k+m))^k
 in terms of k for each pair of m and prevalence in the contour plot. This equation has a unique solution if m + ln(1 - prevalence) > 0. On the other hand, if m + ln(1 - prevalence) < 0, there is no solution to the equation. The solution was found numerically using the uniroot function in R. The expressions for the aggregation indices in terms of m and k is then used to construct the contour plot. Regions of m and prevalence that are inconsistent with a negative binomial distribution are represented as white in the contour plot.  

All contour plots were produced in R using the ggplot2 package <cit.>. The values for m and k reported in <cit.> were heavily skewed and spanned several orders of magnitude with m ranging between 0.1 and 5200 and k ranging between 0.001 and 16.5. To make the plots clearer, log scaling has been applied to these variables. 



§ RELATIONSHIP BETWEEN MEAN ABUNDANCE, K, AND PREVALENCE


The relationship between m, k, and prevalence in wild parasite populations has been examined by several authors with conflicting results <cit.>. While the expression of prevalence in terms of m and k is sufficiently simple to analyse, it is still instructive to construct the contour plot (Fig. <ref> left).  In it, each colour represents a region of values of m and k that give rise to similar values of prevalence. 

We see that prevalence is increasing in both m and k leading to contours that are roughly L-shaped on the range of m and k plotted so prevalence is small when either m or k are small, and prevalence is large when both m and k are large. The contours also show that there is a non-linear relationship between m and k when prevalence is considered fixed. The contours become almost parallel to the horizontal axis as k increases, a consequence of lim_k→∞prevalence = 1 - e^-m. On the other hand, the contours continue to move left as m increases, a consequence of lim_m→∞prevalence = 1.  The contour plot shows that the rate at which prevalence approaches one as m increases is slow when k is small. 


If we restrict our attention to a single-coloured band, i.e. those values of m and k giving rise to similar values of prevalence, we see that, after controlling for prevalence, there is a negative relationship between m and k.  This relationship is forced by the negative binomial distribution, so it will hold true in natural systems to the extent that those systems are well modelled by the negative binomial distribution. The different widths of the contour lines show the non-linearity of the relationship between m, k and prevalence.

The dot points represent estimates of m and k from the 269 parasite-host systems reported in <cit.>. Although several parasite-host systems lie in a region of very high prevalence (both m and k large) or very small prevalence (either m or k small), many others occupy a region of the parameter space where a moderate change in the parameter values would result in a significant change in prevalence assuming a negative binomial distribution.

As Shaw & Dobson reported prevalences in their review, it is possible to compare these with the prevalence values implied by the negative binomial distribution (Fig. <ref> right). In general, there is good agreement; most points within a given contour having the same colour.  This demonstrates the accuracy of the contour plots to interpret relationships in real life situations. The few points where the observed prevalences don’t agree with that determined by the negative binomial could be because these distributions did not conform to a negative binomial.  







§ RELATIONSHIP OF HOOVER & GINI INDICES WITH MEAN ABUNDANCE, K, AND PREVALENCE


Contour plots of the Hoover index and Gini index as functions of m and k are shown in Fig <ref> left and right. The contour plots are qualitatively very similar and share some similarities with the contour plot of prevalence (Fig. <ref>). Both Hoover and Gini indices decrease in both m and k, taking values close to one when either m or k were small, and taking values close to zero when both m and k were large. The contours are L-shaped becoming almost parallel to the horizontal axis as m increases and almost parallel to the vertical axis as k increases. 




The plots show both indices display some stability over a wide range of m and k. Restricting our attention to the Hoover index (Figure <ref> left), we see that for m>5 the value of the index is largely determined by the size of k. For m< 5 the value is less affected by k but more affected by m, as indicated by the number of contours crossed as m decreases. For example, starting from k=1 and m=6, as m decreases the value of the index increases quickly crossing several contours from 0.4 to 1. On the other hand, when m increases from the same point (1,6) the index stays in the same colour band and there is little change in the Hoover value (0.4 to 0.5). For many of the parasite-host systems reported in <cit.>, shown on the figure as dot points, an increase in m, that is moving the points vertically on the contour plot, does not appear to impact the Hoover index since the point would remain in the same-coloured region. On the other hand, in many of the samples, a moderate change in k, that is moving the point horizontally, has a large impact on the Hoover index. Similar behaviour is observed in the contour plot of the Gini index (Fig. <ref> right), with the Gini index appearing to be even less affected by changes in m.

There are two noticeable differences between the contour plots for the Hoover and Gini indices (Fig. <ref>). Firstly, the Gini index is always larger than the Hoover index <cit.> <cit.>. This causes the Gini index to have a smaller range over the region of values for m and k observed in wild populations. Specifically, for the values of m and k reported in <cit.>, the Gini index exceeds 0.9 in 42% (113/269) of cases compared to 20% (54/269) of cases exceeding 0.9 Hoover index. Second, the contours of the Hoover index are not smooth, unlike those of the Gini index. The bumps that occur on the contours of the Hoover index occur at integer values of the mean, the most prominent occurring when the mean is 1. These bumps quickly become much less noticeable as the mean increases.

The contour plots of the Gini and Hoover indices exhibit greater differences when considered as functions of m and prevalence (Fig. <ref>). First, unlike the Gini index, the contour lines of the Hoover index are parallel to the vertical axis when m is less than one. As noted by <cit.>, when all infected hosts harbour infrapopulations larger than or equal to the overall mean, the Hoover index is equal to one minus prevalence. For the negative binomial distribution, this implies the Hoover index is equal to one minus prevalence when the mean is less than or equal to one. Second, there is less variability in the widths of the contours for the Hoover index compared to the Gini index. This suggests the dependence of the Hoover index on prevalence is more regular. At a given m, a change of 0.1 in the prevalence will have roughly the same effect on the value of the Hoover index, regardless of the initial value of prevalence. In contrast, much of the contour plot of the Gini index is coloured yellow, corresponding to values greater than 0.9. Values of the Gini index less than 0.6 are restricted to small region of the plot, indicating that small changes in prevalence in that region will result in a large change in the Gini index.




The parasite data from Shaw and Dobson were taken from five taxonomic groups.  The data, divided into taxa, were superimposed on the plots of m vs k with contour lines of prevalence and Hoover index.  They did not show any obvious grouping.



§ LORENZ ORDER AND THE NEGATIVE BINOMIAL DISTRIBUTION


Both the Hoover and Gini indices are seen in Figure <ref> to be decreasing functions of m and k, as is 1 - prevalence (Figure <ref>). This behaviour is due to how these indices relate to the Lorenz curve and how the parameters m and k affect the Lorenz curve of the negative binomial distribution.

The Lorenz curve of a distribution with cumulative distribution function F is given by 

    L(u) =  ∫^u_0 F^-1 (y)  dy/m,    u ∈ [0,1],

where m is the mean of the distribution and F^-1(x) = sup_x{x:F(x)≤ y} for y∈ (0,1) (Gastwirth, 1971). In our context, the Lorenz curve describes the proportion u of the host population that is infected with a proportion L(u) of the parasite population. When all hosts have the same parasite burden, the Lorenz curve is given by L(u) = u for all u in [0,1]. This is called the egalitarian line. Several indices can be defined in terms of the Lorenz curve. Specifically, the Gini index is twice the area between the Lorenz curve and the egalitarian line, and the Hoover index is the greatest vertical distance between the Lorenz curve and the egalitarian line. Even 1 -prevalence can be viewed as the largest value of u such that L(u) = 0.
The Lorenz curve induces a partial ordering of distributions. Assume F_A and F_B are two distribution functions with finite means. If the Lorenz curve of F_A is greater than the Lorenz curve of F_B for all u, then we say that F_A is smaller than F_B in the Lorenz order and write F_A≤ F_B. This ordering corresponds to the notion of aggregation put forward by <cit.>. From their connections with the Lorenz curve, we see that if F_A≤ F_B, then the Gini and Hoover indices as well as 1-prevalence will be smaller for F_A than for F_B.

The following result shows that the negative binomial distribution decreases in the Lorenz order as m increases and as k increases. 

 
 Let 𝖭𝖡(m,k) denote the negative binomial distribution with parameters m and k. If m_1 < m_2, then 

    𝖭𝖡(k,m_2) ≤_L𝖭𝖡(k, m_1).
 
If k_1 < k_2, then 

    𝖭𝖡(k_2,m) ≤_L𝖭𝖡(k_1,m).


The proof is provided in the Appendix.

The above result explains why Gini and Hoover indices and 1 – prevalence are all decreasing functions of m and k. Figure <ref> also shows that that the contours of both the Gini and Hoover indices become parallel with the axes. This is due to the limiting behaviour of the negative binomial distribution. Depending on how the parameters are allowed to vary, it is known that the negative binomial distribution will converge to either a Poisson distribution or a gamma distribution <cit.>. Fixing m and letting k increase, the negative binomial distribution converges to a Poisson distribution with mean m. This causes the contour lines to become parallel with the horizontal axis as k increases. Similarly, fixing k and letting m increase, an appropriately scaled negative binomial distribution converges to a Gamma distribution with shape and rate parameters both equal to k. Since the Gini and Hoover indices are scale invariant <cit.>, these indices approach their respective values for a Gamma (k, k) distribution as m increases. This causes the contour lines to become parallel with the vertical axis as m increases.




§ DISCUSSION


In choosing the index to use to measure aggregation, those based on Lorenz curves seem to be the favoured, such as the Hoover and Gini.  The Gini returns closer values over a wider range of means, k and prevalence compared to the Hoover, making differences less discernible.  The Hoover has a biological interpretation and may be easier to calculate.  When mean abundances are below one, the Hoover index has restricted values whereas the Gini has no such restriction, suggesting that Gini may be preferred in such a situation. Nevertheless, both indices provide a figure that seems to measure the same phenomenon, a phenomenon that is still undefined.

The contour graphs provide an easily interpreted demonstration of the effects of the various parameters on the Hoover and Gini indices.  These could be deduced by an analysis of the formulae used to calculate the indices but this is not straightforward; indices do not correlate with a particular parameter.  When applying an index to compare aggregation between samples or species, it is useful to know which parameter is having the greatest effect on the index.  The contour graphs provide the answer.

In producing the graphs we calculated indices directly from the parameters of the negative binomial distribution rather than using simulated data as done by <cit.>.  This obviated the need to consider the uncertainty of estimates and the effects of different sample sizes. Our results demonstrated the deterministic functional relationships between the aggregation indices, and the parameters, mean abundance and prevalence. The relationships were not linear indicating that correlation and principal components analysis may not be the best methods to analyse the relationships <cit.>.  


Listing the advantages and disadvantages of Hoover and Gini indices, <cit.> describe them as having the disadvantages of being “strongly negatively correlated with prevalence” and “weakly negatively correlated with mean abundance.” In contrast, the k parameter of the negative binomial distribution and patchiness are described as having the advantages of being “not necessarily correlated with mean abundance” and “only weakly correlated with prevalence.”  These comments ignore the fact that the negative binomial distribution, and hence any index computed on that distribution, is completely specified by the mean and prevalence. In other words, the dependence of any index on mean and prevalence is perfectly deterministic. In fact, the dependence on any pair of quantities that can be used to parameterise the negative binomial distribution, like m and k is perfectly deterministic. 


<cit.> argue that the Gini index is to be preferred over the Hoover index on the basis that Hoover index equals one minus prevalence when the mean is less than or equal to one whereas the Gini index has no such restriction. To decide between the Hoover and Gini indices, if one must choose, then the relationship between these indices and m, k and prevalence need to be considered more closely. Our contour plots (Fig. <ref> & <ref>) have shown other differences in the behaviour of the Gini and Hoover indices. Compared to the Gini index, the Hoover index has a greater range over the region of values for m and k (or prevalence) observed in wild populations and has more regular dependence on prevalence. Given these properties and the Hoover index’s clear biological interpretation, we argue that the Hoover index should be preferred over the Gini index, at least when m is greater than one.


Our analysis has used contour plots to examine how the Gini and Hoover indices are affected by changes in m, k, and prevalence. This approach could, in principle, be applied to construct contour plots from any three indices, provided two of these can be used to parameterise the negative binomial distribution. For example, one could construct a contour plot of the Gini index as a function of VMR and mean crowding as both m and k can be expressed in terms of VMR and mean crowding: 

    m = mean crowding - VMR + 1

and

    k = mean crowding/VMR -1 -1.

The contour plot could then be constructed using the expression for the Gini index in terms of m and k given in Section 2. Further application of contour plots may unravel other complex relationships in ecological parasitology.








§ HOOVER INDEX OF THE NEGATIVE BINOMIAL DISTRIBUTION


In parasitology the negative binomial distribution is usually parameterised in terms of the the mean m and k. The probability mass function is then

    f(x;k,m) =  k+x-1  k-1(k/k+m)^k(m/k+m)^x,    x ∈ℕ_0,

and we write 𝖭𝖡(k,m). Let F(·;k,m) denote the cumulative distribution function of the 𝖭𝖡(k,m) distribution. The first moment distribution of the 𝖭𝖡(k,m) distribution, F^(1)(·;k,m), is 

    F^(1)(x;k,m) = ∑_y≤ x y  f(y;k,m)/m.

For any non-negative integer x

    x  f(x)/m    = x/mk+x-1  k-1(k/k+m)^k(m/k+m)^x
        = x/m(k+x-1)!/(k-1)! x!(k/k+m)^k(m/k+m)^x
    
        
            = (k+x-1)!/k! (x-1)!(k/k+m)^k+1(m/k+m)^x-1
        = (k+x-1)!/k! (x-1)!(k(1+1/k)/(k+m)(1+1/k))^k+1(m(1+1/k)/(k+m)(1+1/k))^x-1
    
            = (k+1) + (x-1) +1  (k+1) -1(k +1 /k+1 + (m+m/k))^k+1(m+m/k/k+1 +(m+m/k))^x-1,

which is the probability mass function of the 𝖭𝖡(k+1,m + m/k) distribution evaluated at x-1. Hence,

    F^(1)(x;k,m)  = F(x-1;k+1,m+m/k).

<cit.> states that the Hoover index can be expressed as

    H = F(m;k,m) - F^(1)(m;k,m).

Hence,

    H=  F(m;k,m) - F(m-1;k+1, m+m/k).





§ PROOF OF THEOREM <REF>


We first recall the definition of convex order, which is closely related to the Lorenz order <cit.>. 


Definition: For random variables X and Y such that 𝔼 ϕ(X) ≤𝔼 ϕ(Y) for all convex functions ϕ : ℝ→ℝ for which the expectations exist. Then we say that X is smaller than Y in the convex order, denoted X ≤_cx Y 



The convex order relates to the Lorenz order in the sense that 

    X/𝔼 X≤_cxY/𝔼 Y

if and only if X ≤_L Y, provided the expectations exist <cit.>  or <cit.>.



For part (a), let X_2∼𝖭𝖡(k,m_2). Conditional on X_2, let X_1∼𝖡𝗂𝗇𝗈𝗆𝗂𝖺𝗅(X_2,m_1/m_2). Then X_1∼𝖭𝖡(k, m_1). As 𝔼 (X_1|X_2) = (m_1/m_2) X_2 and 𝔼X_1  = (m_1/m_2) 𝔼X_2, <cit.> implies (m_1/m_2)  X_2≤_L X_1. Since the Lorenz order is invariant under a change of scale, 𝖭𝖡(k,m_2) ≤_L𝖭𝖡(k,m_1).

For part (b), standard conditioning arguments show that if (N, t≥ 0) is a standard Poisson process and Λ∼𝖦𝖺𝗆𝗆𝖺(α,β) (Gamma distribution with shape parameter α and rate parameter β), then N_Λ∼𝖭𝖡(α, α/β). Let T_i∼𝖦𝖺𝗆𝗆𝖺(k_i, k_i/m). 

It is known that for every convex function ϕ, 𝔼ϕ(N_t) is a convex it t <cit.>. If we can show that T_1≤_cx T_2, then the result will follow from <cit.>. 


By construction 𝔼 T_1 = 𝔼 T_2. Let g_i be the probability density function of 𝖦𝖺𝗆𝗆𝖺(k_i, β_i). Then T_1≤_cx T_2 if g_2 - g_1 exhibits exactly two sign changes in the sequence  +,  -, + <cit.>. As the log function is increasing, log g_2  - log g_1 has the same sequence of sign changes as g_2 - g_1. Then

    log g_2(x) - log g_1(x)   
        = (k_2-1)log x - k_2/m x - ((k_1 -1) log x - k_1/m x ) + C 
        = (k_2 - k_1) log x - (k_1 - k_2)/m x + C,

where C is depends on k_1, k_2 and m but not x. There must be at least one sign change since both g_1 and g_2 integrate to one. For k_2 > k_1 this function is concave so there must be two sign changes. As this function is positive for x → 0 and x →∞ when k_2 > k_1 we have shown T_1≤_cx T_2.  This completes the proof.





4


urlstyle

[Adell and De La Cal(1994)]AdlC:94
J.A. Adell and J. De La Cal.
Approximating gamma distributions by normalized negative binomial distributions.
Journal of Applied Probability, 31:0 391-400, 1994.

[Arnold and Sarabia(2018)]Arnold:87
B.C. Arnold and J.M. Sarabia.
Majorization and the Lorenz order with applications in applied mathematics and economics.
Springer, New York, 2018.

[Bezerra and Bocchiglieri(2023)]BB:2023
R.H.S. Bezerra and A. Bocchiglieri.
Ectoparasitic flies of bats (Mammalia: Chiroptera) in urban green areas of northeastern Brazil.
Parasitol. Res., 122:0 117–126, 2023.

[Crofton(1971)]Crofton:71
H.D. Crofton. 
A quantitative approach to parasitism. 
Parasitology, 62:0 179–193, 1971. 

[Gastwirth(1971)]Gastwirth:1971
J.L. Gastwirth.
A general definition of the Lorenz curve.
Econometrica, 39:0 1037–1039, 1971.

[Gini(1914)]Gini:1914
C. Gini
Sulla misura della concentrazione e della variabilita dei caratteri. 
Atti del R. Instituto Veneto di Scienze, Lettere ed Arti, 73:0 1203–1248, 1914.

[Hankin(2015)]Hankin:2015
R.K.S. Hankin. 
Numerical Evaluation of the Gauss Hypergeometric Function with the hypergeo Package.
The R Journal, 7:0 81-88, 2015.

[Kura et al(2022)]KTCPGA:2022
K. Kura, J.E. Truscott, B.S. Collyera, A. Phillips, A. Garbae and R.M Anderson.
The observed relationship between the degree of parasite aggregation and the prevalence of infection within human host populations for soil-transmitted helminth and schistosome infections. 
Trans. R. Soc. Trop. Med. Hyg., 116:0 1226–1229, 2022. 

[Lester and Blomberg(2021)]LB:2021
R.J.G. Lester and S.P. Blomberg.
Three methods to measure parasite aggregation using examples from Australian fish parasites. 
Methods Ecol. Evol., 12:0 1999–2007, 2021.

[Lloyd(1967)]Lloyd:67
M. Lloyd, 
Mean crowding.
J. Anim. Ecol., 36:0 1–30, (1967).

[Lorenz(1905)]Lorenz:1905
M.O. Lorenz,
Methods of measuring the concentration of wealth. 
Publication of the American Statistical Association, 9:0 209–219, 1905.

[Matos et al(2023)]Matos:2023
I. Matos, D. Silva, J. Oliveira, C. Gonçalves, R. Alves, N. Pereira, F. Catarino, O.M.C.C. Ameixa, J.A. Sousa, L.F. Rangel, M.J. Santos and C. Ayra-Pardo.
Body size-dependent effects on the distribution patterns of phoretic mite species assemblages on Rhynchophorus ferrugineus (Olivier, 1790).
Ecology and Evolution, 13:0 e10338, 2023.

[McVinish and Lester(2020)]ML:2020
R. McVinish and R.J.G Lester.
Measuring aggregation in parasite populations. 
J. R. Soc. Interface, 17:0 20190886, 2020.

[Morato-Moreno(2017)]MM:2017
Morato-Moreno, Manuel.
Origins of the two-dimensional relief representation on some spanish american maps in the sixteenth century.
Boletín de la Asociación de Geógrafos Españoles, 73:0 493-499, 2017.

[Morrill et al(2023)]MPF:2023
A. Morrill, R. Poulin and M.R. Forbes. 
Interrelationships and properties of parasite aggregation measures: a user’s guide.
Int. J. Parasitol., 53:0 763-776, 2023.

[Pennycuick(1971)]Pennycuick:71
P. Pennycuick. 
Frequency distributions of parasites in a population of three-spined sticklebacks, Gasterosteus aculeatus L., with particular reference to the negative binomial distribution. 
Parasitology, 63:0 389-406, 1971.

[Pielou(1977)]Pielou:77
E.C. Pielou.
The measurement of aggregation. In: Mathematical Ecology. 
Wiley Interscience, New York, 1977. 

[Poulin(1993)]Poulin:93
R. Poulin. 
The disparity between observed and uniform distributions: A new look at parasite aggregation. 
Int. J. Parasitol., 23:0 937–944, 1993.

[Poulin(2011)]Poulin:2011
R. Poulin.
Evolutionary Ecology of Parasites: (Second Edition). 
Princeton University Press, 2011. 

[R Core Team(2023)]RCT:2023
R Core Team
R: A language and environment for statistical computing. 
R Foundation for Statistical Computing, Vienna, Austria, 2023.

[Ramasubban(1958)]Ramasubban:58
T.A. Ramasubban.
The mean difference and the mean deviation of some discontinuous distributions.
Biometrika, 45:0 549-556, 1958.

[Rodríguez-Hernández et al(2021)]Rod:2021
K. Rodríguez-Hernández, P. Álvarez-Mendizábal, P. Chapa-Vargas, F. Escobar, F. González-García and D. Santiago-Alarcon.
Haemosporidian prevalence, parasitaemia and aggregation in relation to avian assemblage life history traits at different elevations. 
Int. J. Parasitol., 51:0 365-378, 2021.


[Schweder(1982)]Schweder:1982
T. Schweder.
On the dispersion of mixtures.
Scandinavian Journal of Statistics, 9:0 165–169, 1982.

[Scott(1987)]Scott:87
M.E. Scott. 
Temporal changes in aggregation: a laboratory study.
Parasitology, 94:0 583-595, 1987.

[Shaked and Shanthikumar(2007)]SS:07
M. Shaked and J.G. Shanthikumar.
Stochastic orders.
Springer, New York, 2007.

[Shaw and Dobson(1995)]SD:95
D.J Shaw and A.P. Dobson.
Patterns of macroparasite abundance and aggregation in wildlife populations: A quantitative review. Parasitology, 111:0 S111–S133, 1995.

[Shaw et al(1998)]SGD:98
D.J. Shaw, B.T. Grenfell and A.P. Dobson.
Patterns of macroparasite aggregation in wildlife host populations. 
Parasitology, 117:0 597-610, 1998.

[Taguchi(1968)]Taguchi:68
T. Taguchi. 
Concentration-curve methods and structures of skew populations. 
Annals of Institute of Statistical Mathematics, 20:0 107–141, 1968.


[Tinsley et al(2020)]Tinsley:2020
R.C. Tinsley, H.R. Vineer, R. Grainger-Wood and E.R. Morgan, ER.
Heterogeneity in helminth infections: factors influencing aggregation in a simple host-parasite system. Parasitology, 147:0 65-77, 2020.

[Wade et al(2018)]WFL:2018
M.J. Wade, C.L. Fitzpatrick and C.M. Lively.
50-year anniversary of Lloyd's “mean crowding”: Ideas on patchy distributions. 
J. Anim. Ecol., 87:0 1221–1226, 2018.

[Wickham(2016)]Wickham:2016
H. Wickham. 
ggplot2: Elegant Graphics for Data Analysis. 
Springer-Verlag, New York, 2016.





