A Bimodal Extension of the Beta-Binomial Distribution with Applications

Citation: Reyes, J.; Najera-Zuloaga, J.; Lee, D.-J.; Arrue, J.; Iriarte, Y.A. A Bimodal Extension of the Beta-Binomial Distribution with Applications. Axioms 2024,13, 662. https://doi.org/10.3390/ axioms13100662 Academic Editor: Simeon Reich Received: 24 August 2024 Revised: 17 September 2024 Accepted: 19 September 2024 Published: 25 September 2024 Copyright: © 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). axioms Article A Bimodal Extension of the Beta-Binomial Distribution with Applications Jimmy Reyes 1, Josu Najera-Zuloaga 2, Dae-Jin Lee 3, Jaime Arrué 1and Yuri A. Iriarte 1,* 1Departamento de Estadística y Ciencia de Datos, Facultad de Ciencias Básicas, Universidad de Antofagasta, Antofagasta 1270300, Chile; jimmy.r[email protected] (J.R.); [email protected] (J.A.) 2Department of Mathematics, University of the Basque Country UPV/EHU, 48940 Leioa, Spain; [email protected] 3School of Science and Technology, IE University, 28046 Madrid, Spain; [email protected] *Correspondence: [email protected] Abstract: In this paper, we propose an alternative distribution to model count data exhibiting uni/bimodality. It arises as a weighted version of the beta-binomial distribution, which is defined by a parametric weight function that admits up to two modes for the resulting probability mass function. Like the baseline beta-binomial distribution, the proposed distribution performs well in modeling overdispersed binomial data. Structural properties of the new distribution are studied. Raw moments are derived, which are used to describe the dispersion behavior relative to the mean and the skewness behavior. Parameter estimation is carried out using the maximum likelihood method. A simulation study is conducted in order to illustrate the behavior of the estimators. Finally, two applications illustrating the usefulness of the proposal are presented. Keywords: beta-binomial distribution; bimodality; count data; maximum likelihood; moments; overdispersion MSC: 62E10; 62F10 1. Introduction Count data represents the number of times a particular event occurs in an interval of time, space, or other unit of measurement. This type of data is commonly found in various areas, such as medicine, economics, and engineering, to name a few. For example, Böhning et al. [1] analyzed count data from a dental epidemiological study under the situation of additional zeros. Salman et al. [2] analyzed bankruptcy count data from Swedish small manufacturing firms with the aim of investigating the business failure risk factors of small manufacturing firms. Calabria et al. [3] analyzed the reliability of repairable systems from in-service failure count data. There are many real-world scenarios where the probability of success in binomial experiments cannot be considered constant. For example, the probability of consuming alcohol across the 7 days of a particular week varies from one individual to another (see Alanko and Lemmens [4] ). Considering a beta distribution for the probability of success in a binomial distribution (which gives rise to the beta-binomial distribution) is not overly restrictive since the beta distribution is very flexible in terms of the shapes of its probability density function. A random variable X follows the beta-binomial distribution, denoted X∼BB(n , α , β) , if its probability mass function (p.m.f.) is given by P(X=x) = n xB(x+α,n−x+β) B(α,β),x=0, 1, 2, . . . , n,α,β>0, (1) Axioms 2024,13, 662. https://doi.org/10.3390/axioms13100662 https://www.mdpi.com/journal/axioms Axioms 2024,13, 662 2 of 16 where B(a,b) = R1 0ua−1(1−u)b−1du,a,b>0, is the beta function. In Bayesian inference, the beta-binomial distribution is used to make predictions about the number of successes in future trials, taking into account the uncertainty in the estimate of the probability of success. In classical inference, the beta-binomial distribution can be used to model data with overdispersion in binomial experiments, i.e., when the observed variability is greater than that expected under a standard binomial distribution. A review of the applicability and extensions of the beta-binomial distribution can be found in Wilcox [5] . The use of the beta-binomial distribution in the context of regression is discussed in Crowder [6] . Details on the estimation of the parameters of the beta-binomial distribution can be found in Tripathi et al. [7]. Regarding more recent applications of the beta-binomial distribution, several studies can be found in the literature. To name a few, Palm et al. [8] use the beta-binomial distribution in the formulation of the BBARMA (Beta-Binomial Autoregressive Moving Average) model, which can capture the temporal dynamics and autoregressive structure in count data. Chen et al. [9] use the beta-binomial distribution to propose a GARCH model that captures the variation in the number of new cases of cryptosporidiosis infection, obtaining a useful model for time series data that present bounded counts and high volatility. Jansen and Holling [10] , under a Bayesian approach, use the beta-binomial distribution in the meta-analysis of rare events. Although the beta-binomial distribution is applied in various real-world settings, its performance is not good when empirical distributions exhibit bimodality, i.e., when there are two modes or peaks in the empirical distributions. The presence of bimodality can be explained by the existence of two groups or subpopulations with unique characteristics or by the existence of latent variables that significantly influence the distribution of the population. A very popular methodology in the literature to incorporate flexibility in terms of asymmetry and multimodality is related to the definition of weighted distributions proposed by Fisher [11] and Rao [12] . Suppose that X is a random variable with probability function f(x). The weighted random variable Xwhas PDF fXw(x) = w(x)f(x) µw, (2) where w(·)is a nonnegative weight function and µw=E[W(X)] <∞. A particularly salient case of (2) is obtained when w(x) = x , which defines a lengthbiased distribution. These distributions arise naturally in applied fields, such as reliability and survival analysis, when individuals or mechanical units are sampled with unequal probability due to the experimental design or the existing unequal probability of detection. On the other hand, it is possible to find in the literature weight functions that can lead to multimodality for the weighted distributions resulting from (2). For example, if w(x) = 1 +1−α(x−µ) σ2 , α∈R , and f(x) is the pdf of the normal distribution with mean µ∈R and variance σ2> 0, then (2) reduces to the family of bimodal distributions called the alpha-skew-normal distribution, see Elal-Olivero [13] . Based on the same weight function, Gómez-Déniz et al. [14] introduces a bimodal version of the Poisson distribution. Cortés et al. [15] propose a parametric weight function that involves a power function of exponent 4, which can lead to a probability function with up to three modes. In this paper, we propose an extension of the beta-binomial distribution appropriate to fit overdispersed binomial data that may exhibit both unimodality and bimodality. The proposal arises from (2), using the weight function proposed by Elal-Olivero [13] under a beta-binomial baseline distribution. In this way, the new distribution is aimed at expanding the use of beta-binomial distributions to real-world scenarios where empirical distributions exhibit bimodality. The remainder of the paper is organized as follows. In Section 2, we define the bimodal beta-binomial random variable and study some of its properties, such as the probability Axioms 2024,13, 662 3 of 16 mass function, cumulative distribution function, and the raw moments. The latter are used to describe the behavior of the relative dispersion with respect to the mean and the skewness behavior of the distribution. In Section 3, parameter estimation for the new distribution using the maximum likelihood method is discussed. A simulation study is carried out to evaluate the behavior of the estimators. In Section 4, two application examples with real data are presented to illustrate the usefulness of the proposed distribution. Finally, concluding remarks are presented in Section 5. 2. Bimodal Beta-Binomial Distribution In this section, we derive the new distribution and study some of its main properties. 2.1. Bimodal Beta-Binomial Random Variable The following proposition presents the p.m.f. of the new distribution. Proposition 1. Let X ∼BB(n,α,β)and w(·)be a parametric function given by w(x) = 1+1−q(x−µ) σ2 ,x=0, 1, . . . , n, where µ=nα α+βand σ2=nαβ(α+β+n) (α+β)2(α+β+1), are the mean and variance of X, respectively. Then, the p.m.f. of the weighted random variable Xwis fXw(x;α,β,q) = P(Xw=x) =1 2+q2(1+1−q(x−µ) σ2)n xB(x+α,n−x+β) B(α,β),x=0, 1, . . . , n, (3) such that α,β>0, q ∈Rand B(·,·)is the beta function. Proof. First, we observe that fXw(x)> 0 for all x= 0, 1, 2, . . . , n when α , β> 0 and q∈R . Second, it can be seen that n ∑ x=0 fXw(x) = 1 (2+q2) n ∑ x=0 1+1−q(x−µ) σ2!fX(x) =1 2+q2 2−2q σ n ∑ x=0 (x−µ)fX(x) + q2 σ2 n ∑ x=0 (x−µ)2fX(x)! =1. In consequence, it is concluded that (3) is a valid p.m.f. Definition 1. Let Xw be a random variable with p.m.f. given in (3), then we say that Xw follows a bimodal beta-binomial distribution. We denote this as Xw∼BBB(n,α,β,q). The name given in Definition 1to refer to the new distribution is based on the bimodal behavior that the p.m.f. can present. Figure 1shows some plots of the p.m.f. of the bimodal beta-binomial distribution for different values of its parameters. In the figure, it can be seen that the BBB p.m.f can present a great variety of shapes depending on its parameters: monotonic shape, symmetric/asymmetric unimodal shape, bathtub shape, or asymmetric bimodal shape. A function in the R programming language [ 16 ] for computing (3) is provided in Appendix A. Axioms 2024,13, 662 4 of 16 5 10 15 20 0.00 0.10 x p.m.f. BBB(α = 4, β = 4, q=0) 5 10 15 20 0.00 0.10 x p.m.f. BBB(α = 4, β = 4, q=3) 5 10 15 20 0.00 0.10 x p.m.f. BBB(α = 4, β = 4, q= −3) 5 10 15 20 0.00 0.15 x p.m.f. BBB(α = 1, β = 5, q=0) 5 10 15 20 0.00 0.06 x p.m.f. BBB(α = 8, β = 1, q=0) 5 10 15 20 0.00 0.10 0.20 x p.m.f. BBB(α = 3, β = 8, q=0) Figure 1. Plots of the p.m.f. of the bimodal beta-binomial distribution with n= 20 and different values of α,βand q. 2.2. Two Related Distributions Corollary 1. Let Xw∼BBB(n,α,β,q). Then, 1. fXw(x ; α , β , q= 0 ) = n xB(x+α,n−x+β) B(α,β) , x= 0, 1, 2, . . . , n , α , β> 0, which is the p.m.f. of the beta-binomial distribution. Axioms 2024,13, 662 5 of 16 2. If n =1, then fXw(x;α,β,q) = θx(1−θ)1−x, x =0, 1, such that θ=   β (2+q2)(α+β)h1+1+qpα/β2i,if x =0, α (2+q2)(α+β)h1+1+qpβ/α2i,if x =1. Corollary 1is a direct consequence of (3) considering fixed values for q and n . Part 1 shows that the beta-binomial distribution is a special case of the bimodal beta-binomial distribution obtained when q= 0. The second part shows that the bimodal beta-binomial distribution reduces to the Bernoulli distribution with parameter θ , where θ is a function of the parameters α,β, and q. 2.3. Cumulative Distribution Function The cumulative distribution function (c.d.f.) of Xw∼BBB(n , α , β , q) can be obtained straightforwardly from Proposition 1. Corollary 2. Let Xw∼BBB(n , α , β , q) . Then, the cumulative distribution function (c.d.f.) of Xw is given by FXw(x) = P(Xw≤x) =           0, if ⌊x⌋<0, 1 2+q2 ⌊x⌋ ∑ t=0(1+1−q(t−µ) σ2)n tB(t+α,n−t+β) B(α,β),if 0≤ ⌊x⌋<n, 1, if ⌊x⌋ ≥ n, (4) where ⌊x⌋=max{k∈Z|k<x}, x ∈R. Figure 2shows some plots of the c.d.f. of the bimodal beta-binomial distribution for different values of α , β , and q . As expected, the figure shows that the frequencies are not decreasing as x increases. However, two sharp increases in frequency can be observed in two different intervals of x , which is explained by the bimodal behavior of the corresponding p.m.f. A function in the R programming language for computing (4) is provided in Appendix A. 2.4. Moments The following proposition derives the raw moments of the beta-binomial distribution. Essentially, these moments are expressed as a function of the raw moments of the betabinomial distribution. Proposition 2. Let Xw∼BBB(n,α,β,q). Then, the rth raw moment of Xwis given by E(Xwr)=aµr−bµr+1+cµr+2,r=1, 2, . . . , (5) where a=1 2+q22+2qµ σ+q2µ2 σ2,b=1 2+q22q σ+2q2µ σ2,c=q2 (2+q2)σ2, such that µj=EXj= n ∑ x=0 xjn xB(x+α,n−x+β) B(α,β),j=1, 2, . . . is the jth raw moment of the beta-binomial distribution. Axioms 2024,13, 662 6 of 16 5 10 15 20 0.0 0.4 0.8 x c.d.f. BBB(α = 4, β = 4, q=0) 5 10 15 20 0.0 0.4 0.8 x c.d.f. BBB(α = 4, β = 4, q=3) 5 10 15 20 0.0 0.4 0.8 x c.d.f. BBB(α = 4, β = 4, q= −3) Figure 2. Plots of the c.d.f. of the bimodal beta-binomial distribution with n= 20, α= 4, β= 4 and different values of q. Proof. By definition of expectation, we have that E(Xr w)= n ∑ x=0 xrfXw(x;α,β,q) =1 2+q2 n ∑ x=0 xr(1+1−q(x−µ) σ2)fX(x;α,β), (6) Axioms 2024,13, 662 7 of 16 where fX(x ; α , β) is the p.m.f. of the beta-binomial distribution. Therefore, after some algebra, we see that E(Xr w)=1 2+q22+2qµ σ+q2µ2 σ2n ∑ x=0 x fX(x;α,β)−1 2+q22q σ+2q2µ σ2n ∑ x=0 xr+1fX(x;α,β) +q2 (2+q2)σ2 n ∑ x=0 xr+2fX(x;α,β), and the result is obtained by recognizing the raw moments of the beta-binomial distribution in the above expression. Alternatively, in (6) we can write h1−q(x−µ) σi2=c21−q σcx2 , where c= 1 +qµ σ . Then, using the binomial theorem, we have 1−q(x−µ) σ2 =c22 ∑ k=0 vkxk, with vk=qk σkck. Thus, we can write (7) as E(Xwr)=1 2+q2"n ∑ x=0 xrfX(x;α,β) + c22 ∑ k=0 vk∑ x=0 xr+kfX(x;α,β)# =1 2+q2 µr+c22 ∑ k=0 vkµr+k!,r=1, 2, . . . , where µris the rth raw moment of the beta-binomial distribution. Corollary 3. Let Xw∼BBB(n , α , β , q) . Then, the coefficient of variation ( c . v . (Xw) ) and the Fisher’s skewness coefficient (pβ) of Xware given by c.v.(Xw) = qµ2−bµ3+cµ4−(aµ1−bµ2+cµ3)2 aµ1−bµ2+cµ4 and pβ=aµ3−bµ4+cµ5−3(aµ1−bµ2+cµ3)(aµ2−bµ3+cµ4)+2(aµ1−bµ2+cµ3)3 haµ2−bµ3+cµ4−(aµ1−bµ2+cµ3)2i3/2 , where µ1=Γ(α+1)Γ(α+β)n Γ(α+β+1)Γ(α), µ2=Γ(α+1)Γ(α+β)n(nα+n+β) (α+β+1)Γ(α+β+1)Γ(α), µ3=Γ(α+1)Γ(α+β)n3αβn+β2+3nβ+3n2α+n2α2+2n2−αβ (2+α2+3α+2αβ +β2+3β)Γ(α+β+1)Γ(α), µ4=Γ(α+1)Γ(α+β)A n (6+α3+12αβ +6β2+6α2+3αβ2+3βα2+11α+11β+β3)Γ(α+β+1)Γ(α), µ5=Γ(α+1)Γ(α+β)B n (α+β+3)(α+β+2)(α+β+1)(α+β+4)Γ(α+β+1)Γ(α). such that a, b and c are as in Proposition 2and Axioms 2024,13, 662 8 of 16 A=−αβ +6n3−4αβ2+7nβ2+n3α3+12n2β+β3+18αn2β+7β2nα−4βnα2−5βnα+βα2 +6βn2α2−nβ−β2+6n3α2+11αn3, B=15β3n+50n2β2+10n4α3+60n3β+50n4α−15nβ2+35n4α2+n4α4−30β2α2n−35αn2β +25β2n2α2+15β3nα−10βα3n2+5βα3n+10βn3α3−45β2nα−11αβ3−5βnα+5βα2 −35βn2α2+60n3α2β−α3β−10n2β+75n2αβ2+110n3αβ +β4+11β2α2+24n4−5β3. Figure 3shows some curves of the coefficient of variation and the coefficient of skewness of the bimodal beta-binomial distribution as a function of q under fixed values for α and β . In the figure, it can be seen that the bimodal beta-binomial distribution (depending on q ) can present a greater or lesser relative dispersion (and a greater or lesser skewness level) than the beta-binomial distribution (special case q=0). −4 −2 0 2 4 0.2 0.4 0.6 0.8 1.0 q Coefficient of variation α = 2, β = 2 α = 4, β = 2 α = 2, β = 4 −10 −5 0 5 10 −1.5 −1.0 −0.5 0.0 0.5 1.0 1.5 q Skewness α = 2, β = 2 α = 4, β = 2 α = 2, β = 4 Figure 3. Plots of the coefficient of variation and the skewness coefficient (as a function of q ) of the bimodal beta-binomial distribution with n=20 and different values of αand β. Functions in the R programming language for computing the r th moment (7) and for the coefficients of variation and the coefficient of skewness of Corollary 3are provided in Appendix A. 3. Parameter Estimation In this section, we discuss the maximum likelihood estimator and conduct a simulation study to evaluate the performance of the estimators. 3.1. Maximum Likelihood Estimation Given a random sample X1 , . . . , Xm of the random variable Xw∼BBB(n , α , β , q) , the log-likelihood function for θ= (α,β,q)can be written as Axioms 2024,13, 662 9 of 16 ℓ(θ;xi) = log m ∏ i=1 fXw(xi;α,β,q) =c+ m ∑ i=1 log n xi+ m ∑ i=1 log(x1i) + m ∑ i=1 log Γ(xi+α) + m ∑ i=1 log Γ(n−xi+β), (7) where c=−mlog2+q2+mlog Γ(α+β)−mlog Γ(α)−mlog Γ(β)−mlog Γ(α+β+n) , x1i=1+[1−q(xi−µ)/σ]2and Γ(a) = R∞ 0ua−1e−udu,a>0, is the gamma function. Then, the score functions are given by ∂ℓ(θ;xi) ∂α =c1−2q1+qµ σm ∑ i=1 x2i x1i +2q2 σ m ∑ i=1 xix2i x1i + m ∑ i=1 Ψ(xi+α), (8) ∂ℓ(θ;xi) ∂β =c2−2q1+qµ σm ∑ i=1 x3i x1i +2q2 σ m ∑ i=1 xix3i x1i + m ∑ i=1 Ψ(n−xi+β), (9) ∂ℓ(θ;xi) ∂q=c3+2µ σ21−q σm ∑ i=1 1 x1i −2 σ1−qµ σ+q σm ∑ i=1 xi x1i +2q σ2 m ∑ i=1 x2 i x1i , (10) where c1=−mΨ(α) + mΨ(α+β)−mΨ(α+β+n) , c2=−mΨ(β) + mΨ(α+β)−mΨ(α+ β+n),c3=−2mq/(2+q2),Ψ(a) = ∂log Γ(a)/∂a, with a>0, is the digamma function, x2i=∂ ∂αxi−µ σ=−k1i 2hnα3β(α+β+1)(α+β+n)3i−1/2 and x2i=∂ ∂β xi−µ σ=k2i 2hnα3β(α+β+1)(α+β+n)3i−1/2, such that k1i=−α3xi+α3n− 2 xiα2n−xiα2β+ 2 n2α2+ 2 nα2β+n2αβ +nαβ2+αn2+ nαβ −αxin−xiαnβ+xiαβ2+αxiβ+xinβ+xiβ2+xiβ3+xiβ2n and k2i=−α3xi+α3n+ n2α2−xiα2n−α2xi+nα2+ 2 nα2β−xiα2β+nαβ2−αxiβ−αxin+ 2 nαβ +αn2+xiαnβ+ xiαβ2+xinβ+2xiβ2n+xiβ3. Maximum likelihood (ML) estimator ˆ θ= (ˆ α , ˆ β , ˆ q) of θ= (α , β , q) can be obtained by setting (8)–(10) equal to zero and solving the resulting system of equations. However, due to the analytical complexity of these equations, estimates must be obtained using numerical methods. The standard errors of the ML estimators can be obtained as the square roots of the elements of the diagonal of the matrix K−1(ˆ θ) = −∂2ℓ(θ;xi) ∂θ∂θTθ=ˆ θ−1 , where ∂ℓ(θ;xi)/∂θ∂θTis the hessian matrix. Alternatively, ML estimates can be obtained by solving the optimization problem maxθℓ(θ ; xi) , subject to α> 0, β> 0 and q∈R , where ℓ(θ ; xi) is as in (7). For this, we recommend the use of the function stat:optim() of the R programming language, which also returns the numeric Hessian function. In particular, we consider the L-BFGS-B method [ 17 ], which allows the imposition of box constraints on the parameters. This means that it is possible to specify lower and upper bounds for each parameter, which is very valuable in optimization problems with high dimensions and specific constraints. An R function for computing (7) is provided in Appendix A. Axioms 2024,13, 662 16 of 16 19. Manoj, C.; Wijekoon, P.; Yapa, R.D. The McDonald generalized beta-binomial distribution: A new binomial mixture distribution and simulation based comparison with its nested distributions in handling overdispersion. Int. J. Stat. Probab. 2013,2, 24. [CrossRef] 20. Akaike, H. A new look at the statistical model identification. IEEE Trans. Autom. Control 1974,19, 716–723. [CrossRef] 21. Schwarz, G. Estimating the dimension of a model. Ann. Stat. 1978,6, 461–464. [CrossRef] 22. Ameijeiras-Alonso, J.; Crujeiras, R.M.; Rodríguez-Casal, A. Mode testing, critical bandwidth and excess mass. Test 2019,28, 900–919. [CrossRef] 23. Ameijeiras-Alonso, J.; Crujeiras, R.M.; Rodriguez-Casal, A. Multimode: An R package for mode assessment. arXiv 2018, arXiv:1803.00472. [CrossRef] 24. Xiaohu, L.; Yanyan, H.; Xueyan, Z. The Kumaraswamy binomial distribution. Chin. J. Appl. Probab. Stat. 2011,27, 511–521. 25. Rodríguez-Avi, J.; Conde-Sánchez, A.; Sáez-Castillo, A.; Olmo-Jiménez, M. A generalization of the beta–binomial distribution. J. R. Stat. Soc. Ser. C Appl. Stat. 2007,56, 51–61. [CrossRef] 26. Paul, S. A three-parameter generalization of the binomial distribution. Hist. Philos. Log. 1985,14, 1497–1506. [CrossRef] Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.