Received: 12 October 2022 Revised: 16 April 2023 Accepted: 7 July 2023 DOI: 10.1111/sjos.12677 ORIGINAL ARTICLE Some mechanisms leading to underdispersion: Old and new proposals Pedro Puig1 Jordi Valero2 Amanda Fernández-Fontelo3 1 Departament de Matemàtiques, Universitat Autònoma de Barcelona, Abstract Centre de Recerca Matemàtica (CRM), In statistical modeling, it is important to know the Barcelona, Spain mechanisms that cause underdispersion. Several mech- 2 Escola Superior d’Agricultura de anisms that lead to underdispersed count distributions Barcelona, Universitat Politècnica de Catalunya, Barcelona, Spain are revisited from new perspectives, and new ones 3 School of Business and Economics, are introduced. These include procedures based on Humboldt University of Berlin, Berlin, the number of arrivals in arrival processes, such as Germany renewal and pure birth processes and steady-state dis- Correspondence tributions of birth-death processes, like queues with Pedro Puig, Departament de state-dependent service rates. Weighted Poisson and Matemàtiques, Universitat Autònoma de Barcelona, Centre de Recerca Matemàtica other well-known underdispersed distributions are also (CRM), Barcelona, Spain. related to birth-death processes. Classical and vari- Email: ppuig@mat.uab.cat able binomial thinning mechanisms are also viewed Funding information as important procedures for generating underdispersed Agencia Estatal de Investigación, distributions, which can also generate bivariate count Grant/Award Number: CEX2020-001084-M; Ministerio de distributions with negative correlation. Some example Ciencia e Innovación, Grant/Award applications are shown, one of which is related to Bio- Number: PID2022-137414OB-I00 dosimetry. KEYWORDS arrival process, binomial thinning, birth-death process, COM-Poisson distribution, count distributions, state-dependent service rates, weighted Poisson distribution This is an open access article under the terms of the Creative Commons Attribution-NonCommercial-NoDerivs License, which permits use and distribution in any medium, provided the original work is properly cited, the use is non-commercial and no modifications or adaptations are made. © 2023 The Authors. Scandinavian Journal of Statistics published by John Wiley & Sons Ltd on behalf of The Board of the Foundation of the Scandinavian Journal of Statistics. Scand J Statist. 2023;1–23. wileyonlinelibrary.com/journal/sjos 1 14679469, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sjos.12677 by Readcube (Labtiva Inc.), Wiley Online Library on [24/10/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 2 PUIG et al. 1 I NT RO DU CT ION The Poisson distribution is by far the most widely recognized and commonly used type of distribution for count data analysis. One of its key properties is that the population variance equals the population mean. However, in practice, many count data are overdispersed, that is, the variance is greater than the mean. Conversely, the case of underdispersion (variance lower than the mean) is also common but less frequent. In practice, dispersion is measured with the Fisher dispersion index, defined as the ratio of variance to mean. There have also been some attempts to generalize this index to multivariate count distributions (see Kokonendji & Puig, 2018). Why does over(under)dispersion occur? Of course, there are many count distribu- tions in the literature that are over or underdispersed, but we are not interested in these distributions themselves but in the mechanisms leading to them. There are many mech- anisms explaining overdispersion. For instance, the mechanism of compounding (or clustering) leads to the large family of Compound-Poisson distributions. Similarly, the mechanism of mixing (or clustering) leads to the other large family of Mixed-Poisson distributions. Both mechanisms are understandable and meaningful in many practical situations. In contrast, the mechanisms leading to underdispersion are few, not very well-known, and not always able to explain the cause of the underdispersion observed in a particular dataset. Consider for instance the biodosimetry data described in Pujol et al. (2014). The authors stud- ied the number of a certain type of chromosome aberrations, called dicentrics, in samples of blood lymphocytes when exposed to several doses of X-rays. Table 1 shows the frequency of dicentrics for three different doses. Note that the variances are lower than their respective means. Underdispersion is an extraordinary occurrence in Biodosimetry, as explained by the manual of the International Atomic Energy Agency (IAEA) (IAEA, 2011) on p. 48: Biologically, under- dispersion is very unlikely to occur…(underdispersion) may be indicative of a problem in data sampling. However, Pujol et al. (2014) claimed that the underdispersion detected in their datasets was caused by a real and unexplained phenomenon, and was not due to a sampling problem. They used a weighted-Poisson distribution to fit their data as an empirical solution, but they were unable to explain why underdispersion occurred. In this paper, we will revisit some previously considered mechanisms leading to under- dispersion in the literature, as well as introduce some new ones. Section 2 reviews such mechanisms connected with arrival processes, in particular those based on renewal and pure-birth processes. These are the best known mechanisms although, as we shall see, some of their results have not been correctly understood. Procedures based on the steady state distributions of birth-death processes are analyzed in Section 3. We show that an almost unknown and forgotten proposition by Wise (1962) can explain many of the underdispersed count distributions studied in the past by several authors, including the COM-Poisson and weighted-Poisson distributions (see Sellers & Morris, 2017). Thinning mechanisms and their relation with underdispersion are analyzed in Section 4. In particular, Proposition 6 presents a very meaningful variable thinning mechanism leading to underdispersion. In addition, Section 4.2 explores the relationship between variable thinning, underdispersion and nega- tive correlation. Finally, some conclusions, comments and further research are presented in Section 5. 14679469, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sjos.12677 by Readcube (Labtiva Inc.), Wiley Online Library on [24/10/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License PUIG et al. 3 T A B L E 1 Frequencies of dicentrics among cells for several doses of rays. Number of dicentrics Dose (grays) 0 1 2 3 4 5 6 7 8 9 x s2 3 213 192 85 9 1 0.786 0.641 5 3 23 58 38 15 10 2 1 2.547 1.578 7 4 23 35 35 29 10 9 4 1 4.027 2.697 2 MECHANISMS BA SED O N ARRIVAL PROCESSES An arrival process is a sequence {Tn } of nonnegative and increasing random variables 0 < T1 < T2 < …, representing the times when some phenomenon occurs, called arrival times. Interarrival times are defined as X1 = T1 , and Xi = Ti − Ti−1 for i > 1. The process starts at time 0 and, in principle, multiple arrivals cannot occur at the same time. Let N(t) be the number of arrivals in the interval (0, t]. This count random variable will be our source of interest. The literature on arrival processes and the distribution of N(t) is extensive; see, for example, Cox and Isham (1980) and Ross (1995). We are going to focus on two specific arrival processes that will be described in the following subsections. 2.1 Renewal processes A renewal process is an arrival process in which the interarrival times Xi are iid random variables. Let F be the distribution of X0 , and let F (k) be the distribution of the kth arrival time Tk = X0 + X1 + · · · + Xk−1 . It is clear that F (k) is the k-fold convolution of F with itself. Consider N(t), the number of arrivals in the interval (0, t]. Note that N(t) < k ⇔ Tk ≥ t. Then, P(N(t) < n) = 1 − F (n) (t), and the probability function of the number of arrivals or events is, P(N(t) = n) = P(N(t) < n + 1) − P(N(t) < n) = F (n) (t) − F (n+1) (t). (1) To be coherent with this notation we define F (0) (t) = 1 for t ≥ 0 and F (0) (t) = 0 otherwise. Accord- ing to (1), the distribution of the interarrival times determines the probability function of the number of arrivals. On the other hand, knowledge of the probability function also determines the distribution of the interarrival times because P(N(t) = 0) = 1 − F(t). When the interarrival times are exponentially distributed, the number of arrivals follows a Poisson distribution and, in fact, this is an important and meaningful characterization of the Poisson process. Other distributions have been used for the interarrival times, like Gamma (Winkelmann, 1995; Zeviani et al., 2014), Weibull (McShane et al., 2008; Moriña et al., 2019) or Inverse Gaussian (Levine, 1991; and Seshadri, 1998). In general, if F belongs to a family of distri- butions closed under convolutions (like Gamma and Inverse Gaussian), then the expression for the probability function (1) is tractable. The hazard function (or failure rate function) of an interarrival time distribution F is defined as, f (t) r(t) = , 1 − F(t) 14679469, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sjos.12677 by Readcube (Labtiva Inc.), Wiley Online Library on [24/10/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 4 PUIG et al. where f (t) is the density of F. It is worth mentioning that the hazard function r(t), as well as F(t) and f (t), characterizes the distribution. This function, widely used in Reliability and Survival Analysis, expresses the force of occurrence of events at time t. Distributions where r(t) is increasing(decreasing) for all t, are called Increasing(Decreasing) Failure (or Hazard) Rate dis- tributions, denoted as IFR(DFR). Exponential distribution has a constant hazard function, and hence belongs to both families IFR and DFR. There is an interesting relationship between the behavior of the hazard function of F and the overdispersion or underdispersion of the number of arrivals, stated in the following proposition (Barlow & Proschan, 1965, p. 54): Proposition 1. If F is IFR (DFR), then Var(N(t)) ≤ (≥)E(N(t)). A similar result was rediscovered by Winkelmann (1995), which is only valid asymptotically, and uses different proof. In fact, this proof is very interesting because it permits a different and useful proposition based on the coefficient of variation of F: Proposition 2. Let 𝜇 and 𝜎 be the expectation and standard deviation of F, and c = 𝜎∕𝜇 its coefficient of variation. Then, if c < (>)1, there exists t0 ≥ 0 such that Var(N(t)) ≤ (≥)E(N(t)), for t > t0 . Proof. The proof, outlined in Winkelmann (1995), is a direct consequence of the fact that, as t tends to infinity, E(N(t)) ≈ t∕𝜇 and Var(N(t)) ≈ t𝜎 2 ∕𝜇3 . ▪ Proposition 1 can be directly applied to explore the number of arrivals when F is Gamma or Weibull distributed because their hazard functions are monotone, IFR or DFR being dependent on the value of their shape parameter. However if F follows an Inverse Gaussian distribu- tion, its hazard function is not monotone, because its profile first increases and then decreases (Seshadri, 1998), so Proposition 1 cannot be applied. However, its coefficient of variation has a simple expression so Proposition 2 is useful in this case. Although Proposition 1 only shows a sufficient condition for underdispersion of the num- ber of arrivals, this condition is physically very meaningful and highly explanatory: When the events are more likely to happen as the time between occurrences passes, an underdispersed number of arrivals (occurrences) is produced. Unfortunately, in general, the probability func- tion of the number of arrivals (1) associated with this mechanism does not have a closed form. On the other hand, it is not always possible to express any count distribution as the number of arrivals (for instance at t = 1) of a certain renewal process, and, even when possible, it is very difficult. 2.2 Pure birth processes A pure birth process is an arrival process where the interarrival times Xi are independent and exponentially distributed random variables with parameters 𝜆i (the mean is 1∕𝜆i ). Note that the case where 𝜆0 = 𝜆1 = · · · = 𝜆 corresponds to the Poisson process. A remarkable result of Ball (1995), previously conjectured by Faddy (1994), establishes that the distribution of the number of arrivals is under(over)-dispersed depending on the behavior of the rates 𝜆i : Proposition 3. If the sequence 𝜆0 , 𝜆1 ,…is decreasing (increasing) then, Var(N(t)) < (>)E(N(t)). 14679469, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sjos.12677 by Readcube (Labtiva Inc.), Wiley Online Library on [24/10/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License PUIG et al. 5 The probabilities of the number of arrivals can be calculated from the Kolmogorov forward differential equations, obtaining P(N(t) = 0) = exp(−𝜆0 t) and for n ≥ 1, t P(N(t) = n) = 𝜆n−1 exp(−𝜆n (t − x))P(N(x) = n − 1) dx. ∫0 Considering that time scale is arbitrary and fixing t = 1, we denote pn = pn (1) = P(N(1) = n). This distribution can be expressed in terms of an exponential matrix of the form: (p0 , p1 ,…, pn ) = (1, 0, 0,…, 0) exp(Q), (2) where, ⎛−𝜆0 𝜆0 0 ⋅ 0 ⎞ ⎜ ⎟ ⎜ 0 −𝜆1 𝜆1 ⋅ 0 ⎟ ⎜ ⎟ Q=⎜ 0 0 −𝜆2 ⋅ 0 ⎟ ⎜ ⎟ ⎜ ⋮ ⋮ ⋮ ⋮ ⎟ ⎜ 0 0 0 · · · −𝜆n ⎟⎠ ⎝ It is clear that any sequence of 𝜆n ’s determines a count distribution of the number of arrivals pn ’s. For instance, a linearly increasing sequence of the form 𝜆n = a(b + n) generates the Negative Binomial distribution. A linearly decreasing sequence of the form 𝜆n = a(M − n), n = 0,…, M, generates a Binomial distribution with size parameter M. Unfortunately, except for these linear sequences, the distribution of arrivals pn does not generally have a closed form for any arbitrary sequence of 𝜆n ’s. However, a simple two-parameter distribution can be obtained considering a simple pure-birth process with rates 𝜆0 = a and 𝜆1 = 𝜆2 = · · · = b (Ball (1995)). When a = b this is a Poisson process and, according to Proposition 3, when a > (<)b the distribution of the number of arrivals is under(over)dispersed. Straightforward calculations show that, a p0 = e−a , p1 = a−b (e−b − e−a ) −b abn−1 e−b pn = p a−b n−1 + (a−b)(n−1)! , n ≥ 2. Underdispersion occurs (a > b) when the first birth arises at a higher rate than those coming after. Faddy (1997) considered a family of count distributions given by sequences of the form, 𝜆n = a(b + n)c . When c = 0 this is a constant sequence leading to a Poisson distribution. When c < 0(c > 0), the sequence is decreasing(increasing) and pn is underdispersed(overdispersed). The case of c = 1 corresponds to the Negative Binomial distribution. Podlich et al. (2004) used this family to construct semi-parametric count models that can describe several examples like the number of surviving fetal implants, and the number of cases of Toxoplasmosis. Faddy (1997) showed that for any count distribution 𝜋n there is a sequence of 𝜆n ’s such that 𝜋n = P(N(1) = n). This allows us to represent any count distribution by the number of arrivals of a pure birth process. From the empirical version of Equation (2) we have the system of equations, (p̂ 0 , p̂ 1 ,…, p̂ n ) exp(−Q) = (1, 0, 0,…, 0), 14679469, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sjos.12677 by Readcube (Labtiva Inc.), Wiley Online Library on [24/10/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 6 PUIG et al. Faddy’s representation 6 5 4 λ’s 3 2 1 Dose 7 Dose 5 Dose 3 0 0 1 2 3 4 5 6 7 8 F I G U R E 1 Representation of the distributions of dicentrics for three doses. The error bars represent the first and third quartiles of the bootstrap values. which can be solved numerically for the estimated rates 𝜆̂ n ’s. Figure 1 shows the values of 𝜆̂ n ’s obtained for the distributions of dicentrics shown in Table 1. The uncertainty of the 𝜆̂ n ’s is esti- mated by resampling with 1000 replications. The vertical bars in Figure 1 extend from the first to the third quartiles of the bootstrap values. Note that the interquartile range is large for those 𝜆̂ n ’s corresponding to low frequency observations. In the example shown in Figure 1, Faddy’s representation provides no further information other than a slightly decreasing trend for each dose, which is consistent with the observed under- dispersion. Nonetheless, for other uses, Faddy’s representation can be a meaningful and valuable modeling tool. 3 STEADY STAT E D ISTRIBUTIONS O F BIRTH-DEATH P RO C E S S E S A birth-death process is a special case of a continuous-time Markov process where the states take values equal to 0, 1, 2, …. These processes are used to describe the size of a population, or the number of individuals in a queue, over time. When the process increases by 1, a birth event happens, and when the process decreases by 1, a death event occurs. When the process is at state k (the number of individuals in the queue is k) there are two independent exponential occurrence 14679469, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sjos.12677 by Readcube (Labtiva Inc.), Wiley Online Library on [24/10/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License PUIG et al. 7 times running simultaneously: one with rate 𝜆k that will take us to state k + 1 (birth) if it happens first, and the other with rate 𝜇k which will take us to state k − 1 (death) if it happens first. Let pn (t) the probability that the birth-death process is in state n at time t. These probabilities satisfy the infinite system of differential equations, p′0 (t) = 𝜇1 p1 (t) − 𝜆0 p0 (t) p′n (t) = 𝜆n−1 pn−1 (t) + 𝜇n+1 pn+1 (t) − (𝜆n + 𝜇n )pn (t), n = 1, 2,…, ∑∞ and the normalizing condition n=0 pn (t) = 1. Analytical solutions of these differential equations are difficult and often an impossible task. Nevertheless, it is possible to estimate the parameters of a general birth-death process by observing data at discrete or continuous times (see Crawford et al. (2014) and the references therein). However, consistently with our applications, we will focus the attention on their steady state distributions. A well known result of the theory of birth-death processes establishes that there is a stationary distribution iff, ∑ ∞ 𝜆0 𝜆1 …𝜆n−1 < ∞. (3) n=1 𝜇1 𝜇2 …𝜇n In this case, the probabilities of the stationary distribution are, 𝜆0 𝜆1 …𝜆n−1 pn = p0 , (4) 𝜇1 𝜇2 …𝜇n and, ( )−1 ∑ ∞ 𝜆0 𝜆1 …𝜆n−1 p0 = 1+ , n=1 𝜇1 𝜇2 …𝜇n where pn = limt→∞ pn (t). In general, this distribution is overparameterized as it depends on the ratios 𝜆n ∕𝜇n+1 . But when is this equilibrium distribution underdispersed? To answer this question let us define 𝜃r = (r + 1)pr+1 ∕pr = (r + 1)𝜆r ∕𝜇r+1 , r = 0, 1, 2,…. (5) The behavior of 𝜃r has been used to characterize some families of distributions. For instance, this is the case of the Katz family defined by a linear relationship of the form, 𝜃r = 𝛼 + r𝛽, which comprises Poisson (𝛼 > 0, 𝛽 = 0), Binomial (𝛼 > 0, 𝛽 < 0) and Negative Binomial (𝛼 > 0, 0 < 𝛽 < 1) distributions. The sequence 𝜃r is also related to the Böhning (2016) ratio plot for the case of Poisson distribution. A little known result of Wise (1962) shows an important relationship between the behavior of 𝜃r and the dispersion of pn . Proposition 4. Let N be a random variable following the equilibrium distri- bution of the birth-death process, such that E(N) = 𝜇 and Var(N) = 𝜎 2 . Then, cov(N, 𝜃N ) = 𝜎 2 − 𝜇. 14679469, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sjos.12677 by Readcube (Labtiva Inc.), Wiley Online Library on [24/10/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 8 PUIG et al. A direct consequence of Proposition 4 is a sufficient condition on 𝜃r (5) to assure underdispersion: Corollary 1. Consider a decreasing sequence 𝜃0 ≥ 𝜃1 ≥ 𝜃2 ≥ …. Then, 𝜎 2 ≤ 𝜇. When the sequence is increasing the inequality gets reversed (𝜎 2 ≥ 𝜇). The Wise (1962) result can be used to model many of the underdispersed count distributions studied in the past by several authors, as we shall analyze in the following subsections. 3.1 Queues with state-dependent service rates: the COM-Poisson distribution and others The classic and simplest queuing model, denoted as M/M/1, consists of a queue where a single server attends to certain events (births, jobs, persons, etc) arriving according to a Poisson process at a constant rate 𝜆. It is usually assumed that service times follow exponential distributions that are independent of the interarrival times. Suppose that service times are state dependent, with rates given by 𝜇n = f (n), where n indicates the number of events waiting in the queue. Conway and Maxwell (1961) introduced a distribution that would later be called COM-Poisson (named after R. W. Conway, W. L. Maxwell, and S. D. Poisson), considering f (n) = anc , being studied later by Shmueli et al. (2005). According to (4), the probabilities take the form, 𝜌n pn = p0 , n ≥ 1, (n!)c and, ( )−1 ∑∞ 𝜌i p0 = , i=0 (i!)c where 𝜌 = 𝜆∕a. In this case, sequence (5) remains, 𝜃r ∝ (r + 1)1−c . Consequently, using Corol- lary 1, for c > 1(c < 1) the distribution is underdispersed(overdispersed). Note that for c = 1 this is just the Poisson distribution, which has a linear service-rate function (f (n) = an). Under(over)dispersion is therefore related with a more(less) stressed service-rate function than that of Poisson distribution. Note that COM-Poisson distribution has the property that log(𝜃n ) and log(pn+1 ∕pn ) are linear functions of log(n). In fact, Shmueli et al. (2005) proposed plotting the logarithms of the ratios of successive frequencies against log(n), as a method of validation and parameter estimation (using least squares) of the COM-Poisson distribution. Insufficient attention has been paid to this fruitful methodology for generating underdis- persed count distributions. For instance, consider now service rates of the form f (n) = aebn . Using (4) we obtain the probabilities, 𝜌n p0 = e𝛽0 n−𝛽1 n p0 , n ≥ 1, 2 pn = n(n+1) (6) b e 2 and, ( )−1 ∑ ∞ e𝛽0 i−𝛽1 i 2 p0 = , i=0 14679469, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sjos.12677 by Readcube (Labtiva Inc.), Wiley Online Library on [24/10/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License PUIG et al. 9 where 𝜌 = 𝜆∕a, 𝛽0 = log(𝜌) − b∕2 and 𝛽1 = b∕2. The domain of the parameters is 𝛽1 > 0 and 𝛽0 ∈ R, with the boundary case 𝛽1 = 0 and 𝛽0 < 0, corresponding to the geometric distri- bution. This distribution was first introduced by Gelfand and Dalal (1990), although they wrongly commented that this distribution was always overdispersed. Observe that the sequence 𝜃r ∝ (r + 1)e−2𝛽1 r , r = 0, 1,…, is decreasing for 𝛽1 > log(2)∕2, so, according to Corollary 1, the dis- tribution is underdispersed for this range of 𝛽1 . For 𝛽1 ≤ log(2)∕2 the sequence is not monotonic so Corollary 1 cannot be applied. Anyway, numerical calculation show that when 𝛽1 is close to 0 the distribution is overdispersed. Of course the boundary case 𝛽1 = 0 and 𝛽0 < 0 corresponds to the geometric distribution that is overdispersed. Note that for this distribution the log-ratio log(pn+1 ∕pn ) is just a linear function of n. The following example considers service rates of the form f (n) = anebn , which increase with n much more than in the preceding example. According to (4) we calculate the probabilities, 𝜌n e𝛽0 n−𝛽1 n 2 pn = n(n+1) p0 = p0 , n ≥ 1, (7) n!e 2 b n! with, (∞ )−1 ∑ 𝛽0 i−𝛽1 i2 p0 = (e )∕i! . i=0 In this case the domain of the parameters is 𝛽1 > 0 and 𝛽0 ∈ R, with the boundary case 𝛽1 = 0 and 𝛽0 < 0, corresponding to the Poisson distribution. Now the sequence 𝜃r ∝ e−2𝛽1 r , r = 0, 1,…, is always decreasing except for the boundary case 𝛽1 = 0, so the distribution is underdispersed. This distribution has the property that log(𝜃n ) is a linear function of n. Figure 2 shows the empirical version of log(𝜃n ) against n, for the distributions of dicentrics described in Table 1. The corre- sponding least-squares regression lines point out their linear trends. Consequently, we have fitted the data using the probability function in (7), obtaining the maximum likelihood estimates and the chi-square statistics shown in Table 2. In order to avoid expected frequencies of less than one, the two last cells within doses 5 and 7 were collapsed in the chi-square computation. As usual, the degrees of freedom (df) are calculated as the number of cells minus 1, minus 2 (the number of estimated parameters). If the model does not fit with experimental data, the valid- ity of the assumed service rates is put into question. Observe that 𝛽̂0 increases and 𝛽̂1 decreases with dose. It is interesting to point out that the MLE of the population mean and variance are exactly equal to the sample mean and variance. It happens because the distribution (7) used in this example, and the distribution (6) by Gelfand and Dalal (1990), are two-parameter exponential families with sufficient statistic for both the sum and the sum of the squares of the observations. For characterizations of distributions where the MLE of the population mean is the sample mean see Puig (2003) and Puig and Valero (2006) and the references therein. To include the dose as a covariate, we would need to analyze more experiments with a wide range of doses to investigate the dose dependence of the coefficients 𝛽0 and 𝛽1 . To guarantee that 𝛽1 > 0, a suitable link func- tion could be used. Anyway, the interpretation of these results from the point of view of queues with state-dependent service rates is very meaningful. This is the summary: 1. When the cell is irradiated, the particles (x rays) arrive following a Poisson process, producing damages in the cell. 2. The cell repair mechanism (service) tries to repair the damages. 14679469, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sjos.12677 by Readcube (Labtiva Inc.), Wiley Online Library on [24/10/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 10 PUIG et al. 2 1 log(θn) 0 Dose 7 Dose 5 Dose 3 −1 −2 0 1 2 3 4 5 6 7 8 n F I G U R E 2 Empirical log(𝜃n ) of the distributions of dicentrics against n, with least-squares regression lines, for three X-ray doses. T A B L E 2 MLE of the parameters with standard errors (in brackets) and expected frequencies of the distribution of dicentrics, using the probability function (7); values of the 𝜒 2 statistic and df (in brackets) for several doses. MLE Number of dicentrics: expected frequencies Dose 𝝌2 (grays) ̂0 𝜷 ̂1 𝜷 0 1 2 3 4 5 6 7 8 9 (df ) 3 0.1431 0.1613 208.6 204.8 72.8 12.5 1.2 3.93 (0.1331) (0.0542) (2) 5 1.6691 0.1288 5.4 25.1 45.3 42.1 22.7 7.6 1.6 0.2 9.30 (0.2179) (0.0381) (4) 7 1.9529 0.0642 6.5 18.8 32.1 36.1 28.5 16.5 7.2 2.4 0.6 6.41 (0.1940) (0.0221) (5) 3. If the mechanism is able to repair a damage, the corresponding aberration leaves the queue. The chromosome aberrations that are finally observed (dicentrics) are the number of events in the queue. 4. It is assumed that we have reached a point in the process where the distributions of the number of events will no longer change with time. Then, underdispersion means that the repair mechanism is strongly accelerated by the num- ber of cases of damage. We conjecture that the cell-repair mechanism is then described through function 𝜇n = f (n). In-depth experiments are required to prove this conjecture. 14679469, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sjos.12677 by Readcube (Labtiva Inc.), Wiley Online Library on [24/10/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License PUIG et al. 11 3.2 Weighted Poisson distributions A count random variable is weighted-Poisson distributed with weight w(k; b) if its probability function has the form, w(k; b) e−𝜆 𝜆k P(X = k) = , k = 0, 1, 2,…, C(𝜆, b) k! where C(𝜆, b) is a normalizing constant and b is a parameter (or a vector of parameters). The weight w(k; b) can be interpreted as a kind of recording or sampling mechanism. In fact, it is immediate to see that any arbitrary count distribution can be expressed as a weighted Poisson distribution with a suitable weight. Weighted distributions go back a long way in the literature. See, for instance, Patil et al. (2004) and the references therein. Weighted Poisson distributions have been frequently used to model underdispersed data, using weights for which the normalizing constant C(𝜆, b) takes a closed form. For instance, Pujol et al. (2014) use the weight w(k; b) = 1 + bk2 leading to a normalizing constant C(𝜆, b) = 1 + b(𝜆 + 𝜆2 ), and Cameron and Johansson (1997) use amenable polynomial weights of the form ∑ w(k; b) = (1 + ri=1 bi ki )2 . del Castillo and Pérez-Casany (2005) and Kokonendji et al. (2008), using different methods and notation, state the following result: Proposition 5. Let the weight w(x; b) be a log-concave(convex) function in x ∈ R. Then, the corresponding weighted Poisson distribution is under(over)-dispersed. Interestingly, Weighted Poisson distributions can also be viewed as stationary distributions of birth-death processes, such that w(r + 1; b) 𝜃r = (r + 1)𝜆r ∕𝜇r+1 = 𝜆 . (8) w(r; b) Therefore, a new and simple proof of Proposition 5 arises from the fact that if w(x; b) is a log-concave(convex) function then w(x + 1; b)∕w(x; b) is decreasing(increasing), and the corre- sponding weighted Poisson distribution remains under(over)-dispersed, accordingly to Corol- lary 1. There is also a direct relation between Weighted Poisson distributions and queues with state-dependent service rates, evidenced by considering 𝜆0 = 𝜆1 = · · · = 𝜆 in (8). Then we obtain, w(n − 1; b)n 𝜇n = f (n) = . (9) w(n; b) In fact, any service rate function proportional to f (n) will lead to the same equilibrium distribution family. Kokonendji et al. (2009) presented a nonparametric, kernel-based approach for estimat- ing the weight function. It is important to note that this method would additionally provide a nonparametric way to estimate the 𝜇n service-rate function using expression (9). As an example of how these concepts are related, let us consider the zero-modified Poisson distribution (Johnson et al., 2005), with probability function, p0 = 𝛾, e−𝜆 𝜆k pk = (1 − 𝛾) , k = 1, 2,…. (10) k!(1 − e−𝜆 ) 14679469, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sjos.12677 by Readcube (Labtiva Inc.), Wiley Online Library on [24/10/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 12 PUIG et al. This distribution is widely used to model both underdispersed and overdispersed data. When 𝛾 = e−𝜆 , that is p0 is equal to the proportion of zeros of the Poisson distribution, (10) is just the Pois- son probability function. Now the sequence (5) is 𝜃0 = 𝜙𝜆, 𝜃1 = 𝜃2 = · · · = 𝜃r = · · · = 𝜆, where 𝜙 = (1 − 𝛾)e−𝜆 ∕(𝛾(1 − e−𝜆 )). When 𝛾 < (>)e−𝜆 it means zero-deflation(inflation), 𝜙 > (<)1 and the sequence is decreasing(increasing), and the distribution is under(over)dispersed accordingly to Corollary 1. Note that the zero-modified Poisson distribution is a weighted-Poisson distribution with weight function 𝜔(0) = 1∕𝜙 and 𝜔(k) = 1, k = 1, 2,…. Then, from (9), we find that 𝜇1 = 1∕𝜙 and 𝜇n = n, n = 2, 3,…. There is an interesting interpretation: Poisson distribution arises as the equilibrium distribution of the number of customers in a queue, in a process with constant-rate exponential arrival times and exponential service-time rates proportional to the number of indi- viduals in the queue. The alteration in the service-time rate when there is only 1 individual in the queue (n = 1) has significant consequences. When the service-time is reduced(augmented) 1∕𝜙 < (>)1, that is the rate is augmented(reduced), the equilibrium distribution is Poisson zero-deflated(inflated), being under(over)dispersed. In the following section we are going to present other mechanisms leading to underdispersion that are not directly related with arrival or birth-death processes. 4 THINNING MECHANISMS Many observations in experimental sciences are related to the problem of thinning a random variable X. In this case X itself is not observed, rather a proportion 𝛼 or sampled version of X, X𝛼 , is observed. The standard binomial thinning is defined as follows, Definition 1. Let X be a count random variable, and let 𝜉1 , 𝜉2 , … be iid Bernoulli random variables with probability of success 𝛼 ∈ (0, 1], all of them independent of X. The count random variable, ∑ X X𝛼 = 𝜉i (X𝛼 = 0 if X = 0), i=1 is called a binomial 𝛼-thinning of X. Note that X𝛼 condition on the value of X is binomially distributed, that is, X𝛼 |X = x ∼ Bin(x, 𝛼). Poisson and Mixed-Poisson distributions are closed under independent binomial 𝛼-thinning, in the sense that if X is Poisson or Mixed-Poisson distributed, then X𝛼 belongs to the same family. Several properties regarding closure under binomial thinning are shown in Puig and Valero (2007). The expectation and variance of X𝛼 are, E(X𝛼 ) = 𝛼𝜇, Var(X𝛼 ) = 𝛼 2 𝜎 2 + 𝜇𝛼(1 − 𝛼), (11) where E(X) = 𝜇 and Var(X) = 𝜎 2 . It is plain to see that, ( 2 ) Var(X𝛼 ) 𝜎 =𝛼 − 1 + 1, ∀𝛼 ∈ (0, 1] E(X𝛼 ) 𝜇 and consequently if X is underdispersed(overdispersed) so is X𝛼 . Moreover, the dispersion index of X𝛼 shrinks toward 1. 14679469, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sjos.12677 by Readcube (Labtiva Inc.), Wiley Online Library on [24/10/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License PUIG et al. 13 4.1 Variable thinning Consider now a mechanism of sight such that the proportion of observed events, recorded in the random variable Y , changes with respect to the true value of the hid- den variable X. A simple model consists of assuming a variable binomial thinning as follows, Definition 2. Let X be a count random variable, and let 𝜶 = {𝛼1 , 𝛼2 ,…} be a sequence such that, 0 < 𝛼i < 1. The count random variable, Y , such that Y |X = n ∼ Bin(n, 𝛼n ), denoted as Y = X𝜶 is called a variable binomial 𝜶-thinning of X. Note that when 𝜶 = {𝛼1 , 𝛼1 , 𝛼1 ,…}, it corresponds to the standard binomial thinning X𝛼1 . The probability function of Y is, ∞ ( ) ∑ ∞ ∑ n P(Y = k) = P(Y = k|X = n)P(X = n) = 𝛼nk (1 − 𝛼n )n−k P(X = n). (12) n=0 n=k k The probability generating function (pgf) can be directly derived obtaining, ∑ ∞ ∑ ∞ 𝜙Y (t) = E(t ) = Y E(t |X = n)P(X = n) = Y (1 + 𝛼n (t − 1))n P(X = n). (13) n=0 n=0 The following result connects the behavior of 𝛼n with the under(over)dispersion of the observed events recorded by Y . Proposition 6. Let X be a random variable following a Poisson distribution with expectation E(X) = 𝜆, and Y = X𝜶 . 1. Let 𝛼n be increasing in n. Then Y is overdispersed. 2. Let 𝛼n be decreasing, and n𝛼n increasing in n. Then Y is underdispersed. Proof. Stein-Chen’s identity will be used in the proof. It states that if X follows a Pois- son distribution, and f is any function for which the expected value E(f (X)) exists, then E(Xf (X)) = E(X)E(f (X + 1)). Evaluating the derivatives of the pgf (13) at t = 1, we obtain, ∑ ∞ e−𝜆 𝜆n 𝜇Y = 𝜙′Y (1) = n𝛼n = E(X𝛼X ) = E(X)E(𝛼X+1 ), n=0 n! where X is a Poisson-𝜆 random variable. Similarly, ∑∞ e−𝜆 𝜆n 𝜎Y2 + 𝜇Y2 − 𝜇Y = 𝜙′′Y (1) = n(n − 1)𝛼n2 = E(X(X − 1)𝛼X2 ) = E(X)E(X𝛼X+1 2 ). n=0 n! Because 𝜎Y2 − 𝜇Y = 𝜙′′Y (1) − (𝜙′Y (1))2 , overdispersion or underdispersion depends of 2 the sign of E(X𝛼X+1 ) − E(X)E2 (𝛼X+1 ). Direct algebra shows that 2 E(X𝛼X+1 ) − E(X)E2 (𝛼X+1 ) = cov(X𝛼X+1 , 𝛼X+1 ) + E(𝛼X+1 )cov(X, 𝛼X+1 ). (14) 14679469, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sjos.12677 by Readcube (Labtiva Inc.), Wiley Online Library on [24/10/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 14 PUIG et al. T A B L E 3 MLE of the parameters with standard errors (in brackets) and expected frequencies of the distribution of dicentrics, using the 𝜶-thinning model with X ∼Poisson(𝜆) and logit(𝛼n ) = b∕n; values of the 𝜒 2 statistic and df (in brackets) for several doses. MLE Number of dicentrics: expected frequencies Dose 𝝌2 (grays) b̂ 𝝀̂ 0 1 2 3 4 5 6 7 8 9 (df ) 3 4.1487 0.8685 213.8 197.6 72.5 14.1 1.8 4.54 (1.2628) (0.0643) (2) 5 4.1246 3.4392 5.9 27.5 45.2 38.0 20.9 8.6 2.8 0.8 7.78 (1.1611) (0.3443) (4) 7 5.5545 5.6134 6.0 20.3 33.8 35.3 26.5 15.5 7.4 3.0 1.1 3.95 (1.7332) (0.6002) (6) Consequently, if 𝛼n is increasing 𝛼n+1 and n𝛼n+1 are also increasing and both covari- ances in (14) are positive implying overdispersion. On the other hand, if 𝛼n is decreasing and n𝛼n is increasing, 𝛼n+1 is decreasing and n𝛼n+1 is increasing. Therefore, both covariances in (14) are negative implying underdispersion. ▪ The condition for underdispersion shown in Proposition 6 is very meaningful. A decreasing 𝛼n means, in some sense, that it is more difficult to count a large number of objects than a few. In addition, an increasing n𝛼n means that, although the sight mechanism penalizes counts of a large number of objects, this penalization is not very hard, thus allowing the conditional average E(Y |X = n) = 𝛼n n to be increasing. In other words, it means that as the true number of objects X increases, the number of observed objects also increases on average. For example, consider a model with 𝛼n = a + bn−c , where a, b > 0, a + b < 1 and 0 < c < 1. It can immediately be checked that 𝛼n is decreasing and n𝛼n increasing leading to underdispersion. It is also useful to describe models in terms of logit(𝛼n ). For instance, the model with logit(𝛼n ) = b∕n, b > 0, also satisfies that 𝛼n is decreasing and n𝛼n increasing. This property is also shared using other link functions like probit(𝛼n ) or cloglog(𝛼n ). As an example of application of this distribution, we reanalyze the distributions of dicentrics described in Table 1. The maximum likelihood estimates, expected frequencies, and chi-square statistics are shown in Table 3. It shows a good fit, similar to that obtained in Section 3.1 using a service-rate-model distribution (Table 2). The interpretation is also meaningful: 1. Ionizing particles arrive to the cell accordingly to a Poisson process. It is assumed that the number of damages X is Poisson distributed. The estimated mean of the number of damages increases with the dose of radiation, as expected. 2. Cell repair mechanism acts, and the final number of aberrations is Y = X𝜶 , a variable thinning of the original number of damages X. The action of the repair mechanism, that is 𝛼n , is very similar for doses 3 and 5 grays, and slightly less intense (less efficient) for 7 grays (b̂ = 5.5545). 3. Cell repair mechanism is strongly stressed with the number of damages so that 𝛼n decreases with n. However, given a number of damages n, the expected number of aberrations, n𝛼n , increases. 14679469, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sjos.12677 by Readcube (Labtiva Inc.), Wiley Online Library on [24/10/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License PUIG et al. 15 FIGURE 3 Profile of the dispersion index 𝛿 as a function of 𝛼 and 𝜆. Another interesting model is that obtained from 𝛼n = 𝛽 + 𝛼 n , 𝛽 ≥ 0 and 0 < 𝛼 < 1. It is clear that 𝛼n is decreasing, but n𝛼n is increasing only when 𝛽 > 𝛼 −2∕ log(𝛼) . In particular, when 𝛽 = 0, this leads to a sight mechanism described by 𝛼n = 𝛼 n producing a two-parameter count distribution with pgf, ∑∞ e−𝜆 𝜆n 𝜙Y (t) = (1 + 𝛼 n (t − 1))n . n=0 n! In this case, the expectation and dispersion index take a simple form: 𝜎Y2 𝜇Y = 𝛼𝜆e𝜆(𝛼−1) 𝛿= = 𝛼 3 𝜆e𝜆𝛼(𝛼−1) − 𝜇Y + 1. 𝜇Y Figure 3 shows the profile of the dispersion index as a function of 𝛼 and 𝜆. Note that this model presents over-equi and underdispersion, being an example that the condition of increasing n𝛼n cannot be removed from Proposition 6 to ensure underdispersion. To show an example application of this distribution we analyze one of the classical fertil- ization datasets described in Morgan (1982) and reported in Table 4. This corresponds to an experiment conducted on sea urchins, where a batch of eggs was exposed to sperm at a fixed sperm concentration and a fixed temperature, with a dilute solution of nicotine. Table 4 shows the counts of the number of sperms that had penetrated the eggs at various time intervals for samples of 100 eggs. The row corresponding to time equal to 5 seconds has been removed from the original data here because it is not informative. This dataset was also analyzed by Ridout and Besbeas (2004). We can observe a linear trend profile if we plot the sample means of the number of fertil- izations shown in Table 4 against log-time. Moreover, 𝜇Y ≈ 𝜆 when 𝛼 is close to 1. Therefore, we have fitted the data considering parameter 𝜆 to be time-dependent according to the linear 14679469, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sjos.12677 by Readcube (Labtiva Inc.), Wiley Online Library on [24/10/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 16 PUIG et al. T A B L E 4 Distributions of sperm over eggs. Number of fertilizations Time (seconds) 0 1 2 3 4 5 6 7 x s2 12 82 17 1 0.19 0.1757 20 33 56 9 2 0.80 0.4646 40 17 68 12 2 1 1.02 0.4642 80 9 56 25 7 2 1 1.40 0.8081 180 1 59 27 9 2 1 1 1.60 0.9899 1,8 1,2 1,6 1 1,4 1,2 0,8 Variance 1 Mean 0,6 0,8 0,6 0,4 0,4 0,2 0,2 0 0 0 50 100 150 200 0 50 100 150 200 Time Time FIGURE 4 Empirical and fitted mean (left) and variance (right). relationship, 𝜆(t) = 𝛽0 + 𝛽1 log(t). The likelihood can be directly computed from the expression of the probabilities (12), and using the numerical maximization algorithm provided by func- tion nlm in R we obtain the maximum likelihood estimates: 𝛼̂ = 0.814, se = 0.012; 𝛽̂0 = −3.228, se = 0.460; 𝛽̂1 = 1.410, se = 0.175. The interpretation of this model is meaningful: 1. The number of sperm X contacting an egg follows a Poisson distribution with an increasing mean over time. 2. The number of fertilizations is Y = X𝜶 , a variable thinning of the number of sperm X. Because 𝛼n = 𝛼 n , 𝛼 can be interpreted as the rate of efficiency for the first fertilization. The rate of efficiency for the second fertilization is 𝛼 2 and so on. In other words, the eggs have blocking mechanisms limiting the occurrence of polyspermy and this is the cause of the underdispersion. Figure 4 shows the empirical mean and variance of the number of fertilizations for each time, together with their estimated values, exhibiting good performance. 4.2 Variable thinning and negative correlation As commented earlier, thinning mechanisms are experimentally meaningful. Sometimes, the 𝜶-thinning random variable Y = X𝜶 and the “remainder” Z = X − X𝜶 are both meaningful, lead- ing to a bivariate count random variable (Y , Z), which is a bivariate stochastic decomposition of X since X = Y + Z. 14679469, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sjos.12677 by Readcube (Labtiva Inc.), Wiley Online Library on [24/10/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License PUIG et al. 17 The probability function of (Y , Z) is, P(Y = r, Z = s) = P(Y = r|X = r + s)P(X = r + s) = ( ) r+s r = 𝛼r+s (1 − 𝛼r+s )s P(X = r + s). (15) r Most bivariate count distributions have overdispersed marginals that are positively correlated. However, the following result shows the relation between the dispersion of X and the sign of the correlation of Y and Z, for the classical binomial thinning (𝛼n = 𝛼, ∀n). Proposition 7. Let X be a count random variable with dispersion index 𝛿X = 𝜎X2 ∕𝜇X , Y = X𝜶 , 𝛼n = 𝛼, ∀n, and the remainder Z = X − Y . Then, the correlation coefficient is, 𝛼(1 − 𝛼)(𝛿X − 1) 𝜌YZ = √ . (16) (𝛼 (𝛿X − 1) + 𝛼)((1 − 𝛼)2 (𝛿X − 1) + 1 − 𝛼) 2 Proof. The proof is a direct consequence of (11) and the fact that Z = X1−𝜶 . Then, 𝛿Y = 1 + 𝛼(𝛿X − 1) and 𝛿Z = 1 + (1 − 𝛼)(𝛿X − 1). Moreover, cov(Y , Z) = 𝛼(1 − 𝛼)𝜇(𝛿X − 1), and after some straightforward calculations, the proof is concluded. ▪ Accordingly from expression (16), when X is overdispersed then Y and Z are both overdis- persed, and their correlation is positive. Similarly, when X is underdispersed Y and Z are both underdispersed and their correlation is negative. Note that, if X is under(overdispersed), then so are Y and Z, but both Y and Z are less under(overdispersed) than X. Direct calculations show that the maximum (minimum) correlation is attained at 𝛼 = 1∕2, obtaining |𝛿X − 1| max |𝜌YZ | = . 𝛿X + 1 Commonly used bivariate count models accommodate only a positive correlation between counts. The flexibility for allowing negative correlations for bivariate count distributions generated by the variable thinning stochastic decomposition is shown in the following result. Proposition 8. Let X be Poisson distributed and Y = X𝜶 , such that 𝛼n is decreasing in n. Consider the count random variable Z = X − Y , so that we obtain the decomposition X = Y + Z. Then, if E(𝛼X ) < 1∕2 the correlation is negative: Cor(Y , Z) < 0. Proof. Note that Z is also a variable thinning of X with probabilities 𝛼n∗ = 1 − 𝛼n . Then, 𝛼n∗ is increasing when 𝛼n is decreasing and vice versa. Let us now determine the sign of the covariance, cov(X𝜶 , X − X𝜶 ) = E(X𝜶 (X − X𝜶 )) − E(X𝜶 )E(X − X𝜶 ). Then, E(X𝜶 (X − X𝜶 )) = E(X)2 E(𝛼X+2 (1 − 𝛼X+2 )), according to the reasoning used in the proof of Proposition 6. Note that, E(X𝜶 ) = E(X)E(𝛼X+1 ) and E(X − X𝜶 ) = E(X)E(1 − 𝛼X+1 ). Therefore, cov(Y , Z) = E(X)2 (E(𝛼X+2 (1 − 𝛼X+2 )) − E(𝛼X+1 )E(1 − 𝛼X+1 )), 14679469, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sjos.12677 by Readcube (Labtiva Inc.), Wiley Online Library on [24/10/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 18 PUIG et al. and this is equal to, E(X)2 (cov(𝛼X+2 , (1 − 𝛼X+2 )) + E(𝛼X+2 )E(1 − 𝛼X+2 ) − E(𝛼X+1 )E(1 − 𝛼X+1 )). Note that cov(𝛼X+2 , (1 − 𝛼X+2 )) is negative because 𝛼n is decreasing, so that 1 − 𝛼n is increasing in n. The remaining part is also negative because it can be written as the product, (E(𝛼X+2 ) − E(𝛼X+1 ))(1 − (E(𝛼X+2 ) + E(𝛼X+1 ))), where the first term is negative because 𝛼n is decreasing, and the second term is positive because E(𝛼X+2 ) + E(𝛼X+1 ) < 2E(𝛼X ) < 1. ▪ For instance, condition E(𝛼X ) < 1∕2 can easily be checked for the model 𝛼n = 𝛼 n , 0 < 𝛼 < 1, presented in Section 4.1. It can immediately be seen that in this case E(𝛼 X ) is just the pgf of a Poisson-𝜆 distribution evaluated at 𝛼, so condition E(𝛼X ) < 1∕2 is equivalent to e𝜆(𝛼−1) < 1∕2, or 𝛼 < 1 − log(2)∕𝜆. Useful expressions for computing the expectations, variances and the covariance of Y and Z are also obtained from the proofs of Propositions 6 and 8: 𝜇Y = 𝜆E(𝛼X+1 ) , 𝜇Z = 𝜆 − 𝜇Y 𝜎Y2 = 𝜆E(X𝛼X+1 2 ) + 𝜇Y − 𝜇Y2 (17) 𝜎Z2 = 𝜆E(X(1 − 𝛼X+1 )2 ) + 𝜇Z − 𝜇Z2 cov(Y , Z) = 𝜆2 E(𝛼X+2 − 𝛼X+2 2 ) − 𝜇Y 𝜇Z Propositions 6 and 8 are useful to model bivariate count patterns, as we shall see in the following application. The genome of coronavirus SARS-CoV-2 is about 30,000 “letters” long. The sequence of nucleotides can be downloaded from the National Centre for Biotechnology Information (NCBI) (https://www.ncbi.nlm.nih.gov/). We have worked with the sequence collected in Spain in 2020-03-15, accession number MT359865. Sliding window analyses are a common approach to study the local variation within a genome sequence. We want to explore the presence of pairs of specific trinucleotides in nonoverlapping chunks of size 200 (window length). Consider for instance the number of trinucleotides atc (Y ) and the number of tta (Z) along chunks, and the sum of both X = Y + Z. Here the number of observations for each vari- able is m = 149 (the number of chunks). The sample means and variances are x = 8.0268, y = 2.2550 and z = 5.7718, and s2X = 7.4722, s2Y = 1.9615 and s2Z = 6.8935. The dispersion indexes are dX = 0.9309, dY = 0.8698, dZ = 1.1943, and the sample correlation coefficient between Y and Z is rY ,Z = −0.1880. The negative sign of rY ,Z is an indicator of a kind of competition between both trinucleotides. The Poisson assumption on X has been checked using the bat- tery of 13 goodness of fit tests described by Puig and Weiss (2020). This battery includes some classical tests, like Fisher’s dispersion test, and new tests designed for a large family of overdis- persed alternatives called the LC-class (Puig and Kokonendji (2018)). They explore different aspects of Poissonness. It is reasonable to assume that X follows a Poisson distribution because none of the 13 test rejected the Poisson assumption (the minimum p-value = 0.1211 was for the test Δ∞ ). 14679469, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sjos.12677 by Readcube (Labtiva Inc.), Wiley Online Library on [24/10/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License PUIG et al. 19 0.7 0.6 0.5 Proportion of atc 0.4 0.3 0.2 0.1 2 4 6 8 10 12 14 Total atc + tta FIGURE 5 Observed and fitted values of 𝛼n for the counts of atc. The values of 𝛼n = E(Y |X = n)∕n can be nonparametrically estimated averaging the counts of Y for each observed value of X, n, dividing the average by n. The dots in Figure 5 represent these estimates. To model the bivariate distribution of (Y , Z) we assume that the odds of 𝛼n is a power function of the total number of counts (X), that is, log(𝛼n ∕(1 − 𝛼n )) = a + b log(n). Interestingly, when the slope b is negative, 𝛼n is decreasing with n and n𝛼n is increasing with n. Therefore, Proposition 6 ensures that Y is underdispersed and Z is overdispersed. The likelihood function has the expression, m ( )( )yi ( ) zi ∏ yi + zi 1 e−a−b log(yi +zi ) e−𝜆 𝜆yi +zi L(Y , Z) = . (18) i=1 yi 1+ e−a−b log(yi +zi ) 1 + e−a−b log(yi +zi ) (yi + zi )! Note that parameter 𝜆 can be estimated separately because its MLE is just that of the Poisson distribution obtained with the observations of X (xi = yi + zi ), that is 𝜆̂ = x = 8.0268. The part of the likelihood function corresponding to parameters a and b is that of a logistic regression model, with logit link function, with one covariate (the values of log(X)). Then, standard soft- ware for Generalized linear models can be used, obtaining the MLE, â = 0.4318, se = 0.4117 and b̂ = −0.6434, se = 0.1918. Although the intercept is not significant, note the negative sign of b̂ as expected. Figure 5 also shows the values of 𝛼̂ n = 1∕(1 + e−â−b log(n) ) (red line) reflecting a ̂ ∑∞ −𝜆̂ ̂ n good performance. The estimated value of E(𝛼X ) is n=0 𝛼̂ n e 𝜆 ∕n! = 0.2994 < 1∕2, in agreement with Proposition 8, ensuring that both (Y , Z) are negatively correlated. Expectations, variances, and the correlation coefficient can be estimated using (17) obtaining, 𝜇̂ Y = 2.2559, 𝜇̂ Z = 5.7710, 𝜎̂ 2Y = 1.8142, 𝜎̂ 2Z = 7.2872 and 𝜌Y ,Z = −0.1478. There are globally 2016 different pairs of combinations of trinucleotides. In SARS-CoV-2, only 316 pairs are significantly negatively correlated (𝜌 < −.15) for chunks of size 200. Some of them can be analyzed using the methodology of the previous example (atc and tta), but the Poisson assumption on X is rejected for some other pairs. In these cases, although Propositions 6 and 8 14679469, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sjos.12677 by Readcube (Labtiva Inc.), Wiley Online Library on [24/10/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 20 PUIG et al. could not be illustrated or applied, a similar analysis could be done by using suitable alternatives to the Poisson distribution in the likelihood function (18). 5 DISCUSSION The main message of this paper might be that: if your data are underdispersed try to find out why this is happening, as the answer could be very enriching. We have pointed out two major families of mechanisms leading to underdispersed count data distributions: stochastic processes (arrival and birth-death) and binomial thinning (constant and variable). Renewal processes, such that the interarrival times follow an IFR distribution or have a coef- ficient of variation lower than 1, are an important source of underdispersed distributions of the number of arrivals. In Section 2.1 we reviewed an important result by Barlow and Proschan (1965), pointing out the true meaning of the result by Winkelmann (1995). Pure birth processes are another meaningful source, producing underdispersed distributions when the sequence of birth-rates is decreasing (Ball (1995)). On the other hand, any arbitrary count distribution admits a representation as the number of arrivals (births) of a specific pure birth process (Faddy’s representation). We have revisited these results in Section 2.2 pointing out the importance of Faddy’s representation as a tool for modeling. The steady state distributions of birth-death processes are undoubtedly important and meaningful sources of underdispersed distributions. We saw in Section 3 that decreasing sequences of ratio-plot elements 𝜃n provide underdispersed distributions. In particular, queues with state-dependent service rates provide interesting models leading to distributions like COM-Poisson and that introduced by Gelfand and Dalal (1990). It is worth noting that represent- ing the COM-Poisson distribution as a queuing model with service rate 𝜇n = anc provides for a simple proof, using Corollary 1, that the distribution is underdispersed (overdispersed) for c > 1 (c < 1). The theory of the steady state distributions of queues with state-dependent service rates, and Corollary 1, allows for construction of new underdispersed distributions in a simple way, for instance considering log(𝜃n ) as a linear function of n. This new distribution fits the frequencies of dicentrics very well, providing a biological explanation for why they are underdispersed. Several authors have applied Weighted Poisson distributions for dealing with underdispersion. In Section 3.2 we showed that Weighted Poisson distributions can be interpreted as the steady state distribution of a birth-death process. This makes it possible to prove some theoretical results, also using Corollary 1. Binomial thinning (or subsampling) is a meaningful operation that models sight or selection mechanisms in which an underlying count random variable X is partially observed or selected. We showed in Section 4 that if X is underdispersed, then so is any 𝛼-thinning X𝛼 . However, the variable thinning mechanism defined in Section 4.2 is more flexible and realistic. We show that if the underlying variable X is Poisson distributed, and the sequence of thinnings 𝛼n is decreasing but n𝛼n is increasing, then the variable 𝜶-thinning X𝜶 is underdispersed. The choice of 𝛼n is a source of new underdispersed statistical models. We have reanalyzed the frequencies of dicentrics with the new distribution generated by logit(𝛼n ) = b∕n providing meaningful results. Another new variable thinning model (𝛼n = 𝛼 n ) makes it possible to fit the classic fertilization datasets (Morgan (1982)) providing an interesting interpretation. Two models, one using the service-rate distribution (7) and the other using the Poisson-variable-thinning distribution with logit(𝛼n ) = b∕n, have been used to fit the frequencies of dicentrics. Which model is preferable? Looking at the results presented in Tables 2 and 3, we 14679469, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sjos.12677 by Readcube (Labtiva Inc.), Wiley Online Library on [24/10/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License PUIG et al. 21 can conclude that, statistically speaking, the Poisson-variable-thinning model is slightly better because the value of its global 𝜒 2 statistic is 16.27 (12 df), which is lower than the value for the service-rate model, 19.64 (11 df). However, because both models attempt to explain the cell-repair mechanism, extensive biological experiments are required to determine which model is the most accurate. Finally, Binomial thinning operators provide a natural bivariate stochastic decomposition of the form (X𝜶 , X − X𝜶 ). In Section 4.2, we explored some relationships between the underdis- persion of X𝜶 and the sometimes negative correlation between X𝜶 and X − X𝜶 . In particular, Propositions 6 and 8 allow to construct new bivariate count models, where one marginal is underdispersed, the other is overdispersed and the correlation is negative. We have illustrated some of the results of this section exploring the presence of a specific pair of trinucleotides in nonoverlapping chunks of size 200 of the genome of coronavirus SARS-CoV-2. We have used a new bivariate count distribution generated by a variable thinning where the odds of 𝛼n is a power function of n. Interestingly, this model can be fitted with standard software for Generalized linear models. Many topics remain open to further research. One of these is how to implement regression models to new distributions, like the one presented in (7). One issue is that its parameters can- not be directly interpreted as mean and variance. The COM-Poisson distribution has the same problem. A possible option is to parametrize the distribution in terms of the mean and a suit- able dispersion parameter, albeit this may need the numerical solution of a nonlinear system of equations. In general, each distribution should be investigated to see how its parameters may be affected by the covariates. It can be interesting to investigate order statistics of count distributions as a potential source of underdispersion. A recent article by Badiella et al. (2023), for instance, demonstrates how the Poisson order statistics are underdispersed. Another interesting aspect is that most count distributions linked to capture–recapture prob- lems are overdispersed and belong to the LC-class, a wide family of distributions endowed with a log-convex pgf (Puig and Kokonendji (2018)). Similarly, it can be directly checked that all count distributions with a log-concave pgf are underdispersed, but this property has no practical inter- pretation. What kind of mechanism is behind a log-concave pgf? This challenge is a motivation for further research. ACKNOWLEDGEMENTS The authors thank the two anonymous referees for their valuable comments and suggestions. This work was partially funded by the grant PID2022-137414OB-I00 from the Spanish Ministry of Science, Innovation and Universities and by the Spanish State Research Agency, through the Severo Ochoa and María de Maeztu Program for Centers and Units of Excellence in R&D (CEX2020–001084-M). ORCID Pedro Puig https://orcid.org/0000-0002-6607-9642 REFERENCES Badiella, L., del Castillo, J., & Puig, P. (2023). Ultra log-concavity of discrete order statistics. Statistics & Probability Letters, 201, 109900. https://doi.org/10.1016/j.spl.2023.109900 Ball, F. (1995). A note on variation in birth processes. Mathematical Scientist, 20, 50–55. 14679469, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sjos.12677 by Readcube (Labtiva Inc.), Wiley Online Library on [24/10/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 22 PUIG et al. Barlow, R. E., & Proschan, F. (1965). Mathematical theory of reliability, with contributions by Larry C. Hunter. Reprint of the 1965 original. Classics in Applied Mathematics (Vol. 17). Society for Industrial and Applied Mathematics (SIAM), 1996. Böhning, D. (2016). Ratio plot and ratio regression with applications to social and medical sciences. Statistical Science, 31(2), 205–218. Cameron, A. C., & Johansson, P. (1997). Count data regression using series expansions: With applications. Journal of Applied Econometrics, 12(3), 203–223. Conway, R. W., & Maxwell, W. L. (1961). A queuing model with state dependent service rate. Journal of Industrial Engineering, 12(2), 132–136. Cox, D. R., & Isham, V. (1980). Point processes. In Chapman & Hall/CRC Monographs on Statistics and Applied Probability. Chapman & Hall/CRC. Crawford, F. W., Minin, V. N., & Suchard, M. A. (2014). Estimation for general birth-death processes. Journal of the American Statistical Association, 109(506), 730–747. del Castillo, J., & Pérez-Casany, M. (2005). Overdispersed and underdispersed Poisson generalizations. Journal of Statistical Planning and Inference, 134, 486–500. Faddy, M. J. (1994). On variation in Poisson processes. Mathematical Scientist, 19, 47–51. Faddy, M. J. (1997). Extended Poisson process modelling and analysis of count data. Biometrical Journal, 4, 431–440. Gelfand, A. E., & Dalal, S. R. (1990). A note on overdispersed exponential families. Biometrika, 77(1), 55–64. IAEA. (2011). Cytogenetic Dosimetry: Applications in Preparedness for and Response to Radiation Emergencies. International Atomic Energy Agency. Johnson, N. L., Kemp, A. W., & Kotz, S. (2005). Univariate discrete distributions. In Wiley Series in probability and statistics. Wiley. Kokonendji, C. C., Mizère, D., & Balakrishnan, N. (2008). Connections of the Poisson weight function to overdis- persion and underdispersion. Journal of Statistical Planning and Inference, 138, 1287–1296. Kokonendji, C. C., & Puig, P. (2018). Fisher dispersion index for multivariate count distributions: A review and a new proposal. Journal of Multivariate Analysis, 165, 180–193. Kokonendji, C. C., Senga Kiessé, T., & Balakrishnan, N. (2009). Semiparametric estimation for count data through weighted distributions. Journal of Statistical Planning and Inference, 139, 3625–3638. Levine, M. W. (1991). The distribution of intervals between neural impulses in the maintained discharges of chemically isolated ganglion cells in goldfish retina. Biological Cybernetics, 65, 459–467. McShane, B., Adrian, M., Bradlow, E. T., & Peter, S. (2008). Count models based on Weibull interarrival times. Journal of Business and Economic Statistics, 26(3), 369–378. Morgan, B. J. T. (1982). Modelling polyspermy. Biometrics, 38(4), 885–898. Moriña, D., Serra, I., Puig, P., & Corral, A. (2019). Probability estimation of a Carrington-like geomagnetic storm. Scientific Reports, 9(1), 2393. https://doi.org/10.1038/s41598-019-38918-8 Patil, G. P., Rao, C. R., & Zelen, M. (2004). Weighted distributions. In Encyclopedia of statistical sciences. John Wiley & Sons Inc. Podlich, H. M., Faddy, M. J., & Smyth, G. K. (2004). Semi-parametric extended Poisson process models for count data. Statistics and Computing, 14(4), 311–321. Puig, P. (2003). Characterizing additively closed discrete models by a property of their MLEs, with an application to generalized Hermite distributions. Journal of the American Statistical Association, 98(463), 687–692. Puig, P., & Kokonendji, C. C. (2018). Non-parametric estimation of the number of zeros in truncated count distributions. Scandinavian Journal of Statistics, 45, 347–365. Puig, P., & Valero, J. (2006). Count data distributions: Some characterizations with applications. Journal of the American Statistical Association, 101(473), 332–340. Puig, P., & Valero, J. (2007). Characterization of count data distributions involving additivity and binomial subsampling. Bernoulli, 13(2), 544–555. Puig, P., & Weiss, C. (2020). Some goodness-of-fit tests for the Poisson distribution with applications in Biodosime- try. Computational Statistics & Data Analysis, 144, 1–12. https://doi.org/10.1016/j.csda.2019.106878 Pujol, M., Barquinero, J. F., Puig, P., Puig, R., Caballín, M. R., & Barrios, L. (2014). A new model of biodosimetry to integrate low and high doses. PLOS One, 9(12), e114137. https://doi.org/10.1371/journal.pone.0114137 14679469, 0, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/sjos.12677 by Readcube (Labtiva Inc.), Wiley Online Library on [24/10/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License PUIG et al. 23 Ridout, M. S., & Besbeas, P. (2004). An empirical model for underdispersed count data. Statistical Modelling, 4(1), 77–89. Ross, S. M. (1995). Stochastic processes. In Wiley series in probability and statistics: Probability and statistics (2nd ed.). Wiley. Sellers, K. F., & Morris, D. S. (2017). Underdispersion models: Models that are “under the radar”. Communications in Statistics—Theory and Methods, 46(24), 12075–12086. Seshadri, V. (1998). The inverse Gaussian distribution: Statistical theory and applications. Springer. Shmueli, G., Minka, T. P., Kadane, J. B., Borle, S., & Boatwright, P. (2005). A useful distribution for fitting discrete data: revival of the Conway–Maxwell–Poisson distribution. Journal of the Royal Statistical Society: Series C, 54(1), 127–142. Winkelmann, R. (1995). Duration dependence and dispersion in count data models. Journal of Business and Economic Statistics, 13, 467–474. Wise, J. (1962). The relationship between the mean and variance of a stationary birth-death process, and its economic application. Biometrika, 49, 253–255. Zeviani, W. M., Ribeiro, P. J., Bonat, W. H., Shimakura, S. E., & Muniz, J. A. (2014). The Gamma-count distribution in the analysis of experimental underdispersed data. Journal of Applied Statistics, 41, 2616–2626. How to cite this article: Puig, P., Valero, J., & Fernández-Fontelo, A. (2023). Some mechanisms leading to underdispersion: Old and new proposals. Scandinavian Journal of Statistics, 1–23. https://doi.org/10.1111/sjos.12677