PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Commun Stat Simul Comput. Author manuscript; available in PMC 2018 January 1.
Published in final edited form as:
Commun Stat Simul Comput. 2017; 46(2): 807–814.
Published online 2016 November 11. doi:  10.1080/03610918.2014.960091
PMCID: PMC5736152
NIHMSID: NIHMS900341

The Validation of a Beta-Binomial Model for Overdispersed Binomial Data

Abstract

The beta-binomial model has been widely used as an analytically tractable alternative that captures the overdispersion of an intra-correlated, binomial random variable, X. However, the model validation for X has been rarely investigated. As a beta-binomial mass function takes on a few different shapes, the model validation is examined for each of the classified shapes in this paper. Further, the mean square error (MSE) is illustrated for each shape by the maximum likelihood estimator (MLE) based on a beta-binomial model approach and the method of moments estimator (MME) in order to gauge when and how much the MLE is biased.

Key words and phrases: beta-binomial distribution, intra-correlated binary data, model adequacy, overdispersion

1. Introduction

The beta-binomial model has been widely used as an analytically tractable alternative that captures the overdispersion of a binomial random variable, X, which is a sum of Bernoulli random variables with success probability p and intraclass correlation ρ, since Skellam (1948). When p is assumed to follow a beta distribution, p ~ Beta (α, β) and then unconditional on p, the X follows a beta-binomial distribution (BBD). Here, we denote this as Z. Similar to a beta distribution, a probability mass function of Z takes on a few different shapes: for example, Bell-shape, J-shape, inverse J-shape, and U-shape/Bimodal. Although the success probability and intraclass correlation of Z are identical to those of X, there may exist significant discrepancy in distribution between them as the bimodal shape is more likely to reflect the intraclass correlation between Bernoulli random variables within X. In such cases, inference based on the beta-binomial model approach is biased.

A few methods have been developed for the goodness-of-fit (GOF) test of a beta-binomial model (Brooks et. al, 1997 and Garren et. al, 2000). However, they have been rarely used to investigate the adequacy of the beta-binomial model for X, which may or may not follow a BBD. In this paper the model validation for X is examined. Specifically, we investigate the model validation for each of the classified shapes through the simulation studies in Section 2, along with a real application for screening mammography data. The mean square error (MSE) is illustrated for each shape by the maximum likelihood estimator (MLE) based on a beta-binomial model approach and the method of moments estimator (MME) in order to gauge when and how much the MLE is biased in Section 3. The paper concludes with final comments in Section 4.

2. Simulation Studies: Validation of the beta-binomial model

2.1. Parameter set-up

Consider a discrete random variable Xi [set membership] {0, 1, …, ni} which is a sum of intra-correlated Bernoulli random variables Xij in a sample i, i = 1, …, m, j = 1, …, ni, with success probability p and intraclass correlation, Corr(Xiu, Xiv) = ρ, uv. The mean of Xi is identical to that of a binomial random variable but if ρ > 0, Xi is referred to as an overdispersed binomial random variable because the variance, Var(Xi) = nip(1 − p){1 + ρ(ni − 1)}, exceeds the binomial.

Suppose that conditional on P = p, a discrete random variable Zi [set membership] {0, 1, …, ni} follows a binomial distribution with parameter p and ni and that P follows a beta distribution with two shape parameters, α and β > 0. Unconditional on P, Zi follows a BBD with the probability mass function

fz(zi|ni,α,β)=(nizi)B(α+zi,nizi+β)B(α,β),

where B(α, β) is a beta function. The mean and variance of Zi are

E(Zi)=niαα+β=niμandVar(Zi)=niμ(1μ){1+ϕ(ni1)},

where μ = α/(α + β) and ϕ = 1/(α + β + 1). Letting μ = p and ϕ = ρ, a beta-binomial model has been used as an analytically tractable alternative to the binomial that captures the overdispersion of Xi.

Depending on α and β, the probability mass function of Zi can be largely classified into Bell-shape, J-shape, inverse J-shape, and U-shape/Bimodal as depicted in Figure 1. In simulation studies, 40 sets of parameters (ϕ, μ) of Zi, ϕ = 0.01, 0.025, 0.05, 0.075, 0.1, 0.2, 0.3, 0.5 and μ = 0.5, …, 0.9, are selected to validate the model adequacy; 28 cases from Bell-shape with small to large intraclass correlation (region I), 7 cases from J-shape (region II), and 5 cases from U-shape (region III). We consider the case of having the equal number of the intra-correlated Bernoulli random variables per sample for the subsequent simulation studies (i.e., n = ni for all i).

Figure 1
Selected sets of parameters, (ϕ, μ), ϕ = 0.01, 0.025, 0.05, 0.075, 0.1, 0.2, 0.3, 0.5 and μ = 0.5, …, 0.9, for simulation studies; 28 sets of parameters from Bell-shape with small to large intraclass correlation ...

2.2. Beta-binomial model validation

For simulation studies, the Oman and Zucker method (2001) is employed to generate intra-correlated binomial data, Xi, because their method doesn't require any assumption on a distribution for p, and thus Xi may or may not follow a BBD. For given parameters, μ = p and ϕ = ρ, selected in Section 2.1, the discrepancy between Xi and a BBD is examined through the Pearson's goodness-of-fit (GOF) test (Garren et al., 2000). Note that p and ρ are the parameters of Xi while μ and ϕ are the corresponding parameters of a BBD throughout the paper.

Table 1 shows the descriptive statistics of p-values by the Pearson's GOF test using a Monte Carlo run of 100,000 simulations for ρ = 0.01 ~ 0.3, m = n = 50. The results of ρ = 0.5 are similar to that of ρ = 0.3 and are omitted in Table 1. Q1, Q3, and “PP” indicate 25th and 75th percentile of p-values and the proportion of p-values < 0.05, respectively. The p-values are not uniformly distributed except when ρ is very small, and the proportion of p-values < 0.05 increases as (i) ρ increases for a fixed p and (ii) p departs from 0.5 for a fixed ρ. The results in Table 1 and Supplemental Figure A reveal that the discrepancy between Xi and a BBD with μ = p and ϕ = ρ is manifest as the set of parameters, (ρ, p), falls into or gets close to region II (J-/Inverse J-shape) or III (U-shape/Bimodal shape) in Figure 1.

Table 1
Descriptive statistics of p-values by the Pearson's GOF test using a Monte Carlo run of 100,000 simulations when m = n = 50.

With n = 50, four sets of parameters are selected for visual evaluation, depending on shape of the BBD; (A) Bell-shape with small correlation (μ = 0.5, ϕ = 0.01), (B) Bell-shape with moderate correlation (μ = 0.5, ϕ = 0.1), (C) J-shape (μ = 0.75, ϕ = 1/3), (D) U-shape (μ = 0.5, ϕ = 0.4). Figure 2 depicts the comparisons of the distribution of Xi to the binomial and the beta-binomial mass functions. The beta-binomial model fits well for the small value of correlation (Figure 2 (A)), while for the other three cases (Figure 2 (B) – (D)) substantial discrepancy is observed. The shapes of Xi for (B) – (D) are bimodal while the shapes of the beta-binomial mass functions are Bell-, J-, and U-shape, respectively. Thus, it appears that the bimodal shape occurs when an intraclass correlation is moderate to large, so for these cases inference based on the beta-binomial model approach may be biased.

Figure 2
Comparisons of the distribution of Xi to the beta-binomial and the binomial mass functions by the shape of the beta-binomial mass function; (A) Bell-shape with small correlation (μ = 0.5, ϕ = 0.01), (B) Bell-shape with moderate correlation ...

2.3. Real Example: Screening Mammogram

The model validation for the data of a screening mammogram study (Beam et al., 2003) is considered as a real example. In their study, 64 out of 148 films were used for estimating sensitivity while 84 were used for specificity. One hundred eight radiologists each read the 148 mammography films (i.e., m=108, n=64 for sensitivity and n=84 for specificity). The distributions of the total number of positive mammography readings out of confirmed cancer cases and negative readings out of non-cases are depicted and compared to the binomial and the beta-binomial mass functions in Figure 3 (A) and (B), respectively. The parameters of the BBD are estimated by the MLE (denoted as [mu] and [phi with hat]) developed by Kim and Lee (2013). It is shown that the data are overdispersed, and the beta-binomial model provides a superior fit to the data, compared to the binomial model (the p-values by the Pearson's GOF test are < 0.001). The estimates of sensitivity and intraclass correlation, [mu] and [phi with hat], are 0.874 and 0.054, respectively, and the model fits well (p-value by the GOF test = 0.659). For specificity, [mu] and [phi with hat] are 0.786 and 0.086, and no significant difference between the distribution of negatives to non-cases and a BBD is found (p-value = 0.394) but it appears that the shape of the specificity data is slightly bimodal. Thus, these results support the simulation study results presented in Table 1 and Figure 2.

Figure 3
Distributions of Positives to Cases (A) and Negatives to Non-Cases (B) with comparison to the binomial and the beta-binomial mass functions.

3. Comparisons of the MLE to the MME

Based on the first and second moments of Xi, a method of moments estimator (MME; denoted as [p with tilde] and [rho with tilde]) for p and ρ,

p=1mi=1mXi/n,ρ=1n1[nSm12p(1p)1],Sm12=1m1i=1m(Xi/np)2,
(1)

can be considered as a distribution-free point estimator as no particular distribution for Xi is assumed. In simulation studies, the MLE is used to estimate μ and ϕ of Zi and is compared to the distribution-free MME, [p with tilde] and [rho with tilde], based on (1) in order to gauge when and how much the MLE is biased for each of the corresponding model shapes.

The comparisons between the MLE and the MME are implemented through a Monte Carlo run of 20,000 simulations with m = 25, 50, 100, and n = 25, 50, 100. Table 2 illustrates the mean and the mean squared error (MSE) for each estimator when m = n = 25, 50, and 100. The results of other n and m are similar to Table 2 (results are omitted). The MSE decreases with the number of samples, m, and the number of Bernoulli variables in each sample, n. As Kleinman (1973) pointed out, the MLE is more efficient than the MME when the MLE is unbiased or very close to the true parameter value. The simulation study results show that (i) the distribution-free MME are more reliable and robust across all regions in Figure 1, (ii) the MLE for success probability, [mu], is slightly biased for (p, ρ) = (0.7, 0.5) and (0.9, 0.5), (iii) the MLE for intraclass correlation, [phi with hat], is significantly biased except when ρ is very small. Specifically, [phi with hat] is underestimated for parameters which fall into or get close to regions II (J-/Inverse J-shape) or III (Bell-/U-shape) in Figure 1. A 15% or more biased estimate, compared to the true parameter value (i.e., (estimate - true value)/ (true value) > 0.15) is highlighted with boldface in Table 2. Supplemental Figures B1 and B2 illustrate the histogram of the MLE, [phi with hat], for ρ = 0.05, 0.1, 0.3, and 0.5 using a Monte Carlo run of 20,000 when m = n = 100 and p = 0.5, 0.7, and 0.9.

Table 2
Comparisons of the MLE to the MME for p and ρ using a Monte Carlo run of 20,000 simulations when m=n=25, 50, and 100.

4. Discussion

In this paper, the beta-binomial model validation is examined using the random samples Xi which are generated by the Oman and Zucker method. The discrepancy between the distribution of Xi and a BBD with μ = p and ϕ = ρ is examined through the Pearson's GOF test for each of the classified shapes. The p-values by the Pearson's GOF test are not uniformly distributed except when ρ is small, and the proportion of p-values < 0.05 increases as (i) ρ increases for a fixed success probability, p and (ii) p departs from 0.5 for a fixed ρ. The simulation study results presented in Table 1 and and22 reveal that the discrepancy between them is manifest and the MLE for intraclass correlation, [phi with hat], is significantly underestimated when the values of true parameters (p, ρ) fall into or get close to region II (J-/Inverse J-shape) or III (U-shape/Bimodal). As the bimodal shape is more likely to reflect ρ within each sample, the beta-binomial model may not be appropriate for such cases. However, the beta-binomial model may provide a better fit to the data with small ρ (Bell-shape), the large zero frequency data (Inverse J-shape), or large n frequency data (J-shape). For example, the screening mammogram data discussed in Section 2.3 have small ρ for sensitivity and specificity and thus are well fitted to a BBD. Skellam (1948) and Chatfield and Goodhardt (1970) proposed the beta-binomial model fit to the traffic accident data (Xi = the number of traffic accidents) and the consumer purchasing behavior data (Xi = the number of weeks out of n in which the consumer makes at least one purchase), respectively. Both data have large zero frequency and decreasing pattern (inverse J-shape), and the beta-binomial model showed a better fit to the data, compared to the Poisson model or the Negative binomial model. However, no interpretation on ρ in those datasets was provided.

In summary, we would recommend that both the Pearson's GOF test and the distribution-free estimate be conducted prior to making parametric inference under the beta-binomial model. Specifically, if the distribution-free estimate for ρ, [rho with tilde], is greater than 0.1, the GOF test should be considered to investigate if there is a significant evidence that the beta-binomial model poorly fits the data, which could result in biased inference. When the beta-binomial model is not valid for a given data, one may consider an alternative model including the negative binomial model. However, finding a valid model should be differentiated for the given data along with solid model validation process.

Supplementary Material

Supplemental Figure A. Histograms of p-values by the Pearson's GOF test using a Monte Carlo run of 100,000 simulations when m = n = 50, p =0.5 (top) ~ 0.9 (bottom 3 panels), and ρ =0.01 (left), 0.05 (middle), and 0.075 (right 5 panels)

Supplemental Figure B1. Histogram of the MLE, [phi with hat], for ρ = 0.05 (left) and 0.1 (right panels) using a Monte Carlo run of 20,000 simulation when m = n = 100 and p =0.5 (top), 0.7 (middle), and 0.9 (bottom panels).

Supplemental Figure B2. Histogram of the MLE, [phi with hat], for ρ = 0.3 (left) and 0.5 (right panels) using a Monte Carlo run of 20,000 simulation when m = n = 100 and p =0.5 (top), 0.7 (middle), and 0.9 (bottom panels).

Acknowledgments

The authors expresses many thanks to Dr. Craig Beam for allowing the use of the screening mammogram data. This work was partially supported by NIH/NCI grant 5P30CA076292-16.

References

  • Beam C, Conant E, Sickles E. Association of volume and volume-independent factors with accuracy in screening mammogram interpretation. Journal of the National Cancer Institute. 2003;95:282–290. [PubMed]
  • Brooks S, Morgan B, Ridout M, Pack S. Finite mixture models for proportions. Biometrics. 1997;53:1097–1115. [PubMed]
  • Chatfield C, Goodhardt G. The Beta-Binomial Model for Consumer Purchasing Behavior. Applied Statistics. 1970;19:240–250.
  • Garren S, Simth R, Piegorsh W. On a likelihood-baed goodness-of-fit test of the beta-binomial model. Biometrics. 2000;56:947–950. [PubMed]
  • Kim J, Lee JH. Simultaneous confidence intervals for a success probability and intraclass correlation, with an application to screening mammography. Biometrical Journal. 2013;55(6):944–954. [PubMed]
  • Kleinman JC. Proportions with extraneous variance: single and independent sample. Journal of the American Statistical Association. 1973;68:46–54.
  • Oman SD, Zucker DM. Modelling and generating correlated binary variables. Biometrika. 2001;88:287–290.
  • Skellam JG. A probability distribution derived from the Binomial distribution by regarding the probability of success as variable between the sets of trials. Journal of the Royal Statistical Society, B. 1948;10:257–261.