{'fraction_non_alphanumeric': 0.02917922354192519, 'fraction_numerical': 0.016060982095373162, 'mean_word_length': 5.0, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Performance of the Beta-Binomial Model for Clustered Binary Responses: Comparison with Generalized Estimating Equations

This study examined performance of the beta-binomial model in comparison with GEE using clustered binary responses resulting in non-normal outcomes. Monte Carlo simulations were performed under varying intracluster correlations and sample sizes. The results showed that the beta-binomial model performed better for small sample, while GEE performed well under large sample.


Introduction
A series of item responses coded as 1 or 0 is often observed in various research studies in the social and behavioral sciences. Binary responses usually indicate pass/fail, agree/disagree, yes/no, and correct/incorrect. Subjects' traits measured by these items are often summarized by averaging or summing up the binary values. In educational research, the normal linear models (e.g., linear regression, ANOVA, mixed linear model) have been commonly accepted to analyze the summed scores. They can be reasonable choices when the central limit theorem holds. However, this standard approach may be limited in some situations where the normal models cannot properly approximate binomial probabilities of summary scores, often with skewed distributions (Warton & Hui, 2011). In addition, the fact that the summed test scores bound within a lower-to-upper limit is contradictory to the characteristic of a normal distribution having an infinite range, possibly leading to predicted SEONGAH IM 3 values that are out of range. This could be an important consideration in educational research especially when predicted scores and estimates from a statistical model are utilized for making decisions on teachers and students. If normal linear modeling using summed or aggregate scores is not appropriate, models based on the binomial distribution can be an alternative that can address concerns related to non-normality and limited range (Ferrari & Comelli, 2016;Jaeger, 2008).
When item responses constitute total scores or binomial counts in binomial modeling, it should be noted that repeated responses within an individual respondent are more alike than those from different respondents, thus leading to extra-binomial variation, a feature known as overdispersion. When an analysis ignores dependence within a subject or cluster (i.e., intracluster correlation), the results can produce biased estimates and distort standard errors, leading to invalid inferences (Jaeger, 2008). Thus, a study that proposes appropriate statistical approaches to properly deal with correlated binary responses resulting in skewed outcomes would be necessary and important. Among the statistical approaches that are capable of dealing with the concerns mentioned above, this study employed the beta-binomial model (Williams, 1982) and examined its performance using Monte Carlo simulation.

Studies Using the Beta-Binomial Model
The beta-binomial model is a commonly used approach for overdispersion that is often a manifestation of intracluster correlation. The beta distribution allows a binomial parameter to be heterogeneous and enables the beta-binomial model to explain various shapes of distributions with negative to positive skewness. Empirical studies using the beta-binomial regression with clustered binary data can be found in some fields including public health, economics, and ecology. Among them, some relevant examples are studies in health sciences that used survey data to fit varying shapes of health-status survey scores from different subdomains of health (e.g., Arostegui et al., 2007;Khan & Morris, 2014;Lamu & Olsen, 2018). In these studies, the beta-binomial regression worked well for the aggregate scores that had various non-normal distributions, compared to the normal linear models or some nonparametric approaches such as those with extended estimating equations (Basu & Rathouz, 2005) and quantile regressions. Recently, Najera-Zuloaga et al. (2017) suggested the extend use of the beat-binomial model in the context of the longitudinal design to identify risk factor for chronic disease patients and found that the application leaded to clinically relevant results. The main interest of the beta-binomial analysis is on the marginal distributions of aggregate counts or proportions, rather than modeling of initial binary observations at data collection. When the binary responses from the initial data collection are available and the marginal scores in relation to covariates are of interest, a population averaged or marginal model such as the generalized estimating equations (GEE) of Liang and Zeger (1986) would be desirable because taking the summary scores in each cluster means information regarding the trials or items is lost in the aggregation process. GEE analyzes the mean response and the within-cluster dependence, assuming the primary interest is in the former and regarding the latter as a nuisance that must be taken into account for valid inference. The parameter estimates of GEE are known to be consistent even when the working correlation matrix reflecting dependence within clusters is misspecified (Wang & Carey, 2003).
Often researchers who employ the GEE contrasted it with generalized linear mixed models (GLMM). The studies noted that GEE produced slightly overestimated parameters and GLMM reduced power. Type 1 errors were elevated with small sample sizes in both GEE and GLMM, while either of the two methods clearly outperformed (Hallgren et al., 2016;Hubbard et al., 2010). While both approaches can adeptly handle clustered outcomes, they address dependency differently. The way to interpret the regression coefficients of the two models are different, and the size of the coefficients also substantially differ depending on the extents of correlation within clusters (Hu et al., 1998;Lee & Nelder, 2004). Rather the GEE regression estimates are similar to those from the beta-binomial model and other quasi-likelihood binomial models with beta-binomial type variance (Agresti, 2013). The nuisance parameter, i.e., the within-cluster correlation, can be compared to intracluster correlation estimated from the beta-binomial model.
The beta-binomial model was also compared to GLMM with aggregate data in a few studies. Harrison (2015) found that the beta-binomial model performed well even with data sets generated under random intercepts GLMM with experimental factors in ecology. Ferrari and Comelli (2016) also found that the beta-binomial model was powerful and outperformed GLMM in a range of experimental conditions. Yet, studies that systematically compare the betabinomial model and GEE have been rare. In an empirical study in health service, Dilba and Aerts (2004) compared the beta-binomial model and GEE for clustered data and found that Type 1 error rates of the two methods were close under an exchangeable correlation matrix, yet the generalizability of the finding was questionable as the study was performed using an empirical data set.

5
In education, early researchers (e.g., Carlin & Rubin, 1991;Gross & Shulman, 1980;Huynh, 1976) used the beta-binomial distribution to depict varying distributions from criterion referenced testing. Later, Fox (2008) applied the model to a multivariate randomized response data. These studies focusing on the depiction of the non-normal outcomes, however, did not aim to examine regression modeling of aggregate or clustered scores in relation to a set of covariates. Knowing that the beta-binomial model could be suitable for modeling non-normal test scores from dependent binary responses, yet the systematic examination compared with a wellknown marginal model has been rare, this study investigated the performance of the beta-binomial model in comparison with GEE using correlated binary responses simulated under three experimental conditions that are relevant in real world testing situations. In the next section, GEE with logit link and the beta-binomial logistic model are described, followed by the simulation design, and results of analyses from the two models.

Two Models Generalized Estimating Equations
The Generalized Estimating Equations of Liang and Zeger (1986) was designed to develop a marginal or population-averaged model that tests influence of covariates on exponentially distributed (e.g., Poisson, Binomial, Gamma) response variables, and facilitates analysis of data collected in repeated, longitudinal, and panel designs. GEE is an extension of generalized linear model (GLM: McCullagh & Nelder, 1989), in which a specific type of correlation structure is incorporated into the variance function. The dependence among response variables across trials is specified as a working correlation matrix that accounts for the form of withinsubject correlation.
For a case of binary responses, yij denotes a response of subject i at repeated trials j with either taking the value of 1 for correct or positive endorsement or 0 for incorrect or negative endorsement, where i = 1,…, N and j = 1,…, m. The GEE with logit link is (1) As likelihood estimation methods require that each observation in the model is independent of the others, GEE uses a quasi-likelihood method for estimating parameters using the first two moments and produces parameter estimates that are asymptotically normal (Liang & Zeger, 1986). A consistent estimator of the asymptotic covariance matrix of β, the solution of the equation (1), is obtained by the robust (also often called sandwich) estimator. The GEE parameter estimates are more efficient and consistent when the number of clusters increase even under misspecification of working correlation matrix. In this study, exchangeable correlation matrix was used, assuming the correlation among responses between any two observations is a constant because the correlation structure is appropriate to use with level-1 clustered data (Hilbe, 2009, p. 447).

Beta-Binomial Model
The beta-binomial model (Williams, 1982) is a well-known conjugate mixture model for clustered binary data. It applies with totals or aggregates from binary items or trials. Suppose 1 ; 1, 2, , ; 1, 2, , where the Bernoulli random variable yij takes value 1 with a probability p and 0 with a probability q = 1p, and Yi = 0, 1,…, mi. Under conventional binomial model, Yi ~ Binomial(mi, p) and the correlation coefficient between yij and yik is 0 for every j ≠ k. When the independence assumption is violated, some alternative approaches to deal with the dependence should be considered. The beta-binomial model can be legitimately chosen to revamping extra-binomial variation or overdispersion due to intracluster dependence. The beta-binomial assumes that a success probability pi for cluster i (i.e., a subject i) is beta-distributed with two shape parameters α > 0 and β > 0. The appropriateness of the beta-binomial distribution 7 is dependent upon how well the beta distribution can represent the population of pis. The versatility of the beta distributions with different combinations of α and β allows one to handle various shapes of distributions. Its probability mass function of beta-binomial distribution is where B(., .) is the beta function. Let πi = αi / (αi + βi) satisfy a logistic relationship with a predictor vector X T . Then Yi follows a beta-binomial distribution with the mean and the variance as follows: ( ) where intracluster correlation ρi = (αi + βi + 1) -1 . As shown in equation (3), the beta-binomial form of expectation is the summation of individual means that is the same as that for binomial distribution, but the variance (4) is different from the binomial variance, πi(1 -πi). The multiplicative component of [1 + ρi(mi -1)] in (4) specifically designates the role of statistical dependence among trials or items (Agresti, 2013;McCullagh & Nelder, 1989). The mean equation for the binomial logistic regression model with the grouped or summed outcome variable is where X T is a transposed vector of explanatory variables and β is a vector of regression parameter estimates. The parameter ρi can be formulated as a function of predictors, but more typically it is taken to be a constant and thereby restricting αi + βi to be constant for all i (Simonoff, 2003). In this study, the intracluster correlation was modeled as a constant, assuming correlations among responses are exchangeable, so as to compare results of beta-binomial model to those from GEE.
THE BETA-BINOMIAL MODEL FOR CLUSTERED BINARY RESPONSES 8

Method Experimental Factors of Simulation Design
The simulation design of this study included two types of non-normal marginal distributions, negative and positive skewed distributions. Three experimental factors were manipulated for each type of distributions, totaling 36 experiments. Data generation in each experiment was performed with 1000 repetitions.

Intracluster Correlation.
The impact of clustering on statistical results crucially depends on the strength of intracluster correlation. Specifically, if the intracluster correlation is relatively strong, the failure to take clustering into account is likely to have more profound effect on the outcomes of statistical analysis (Galbraith et al., 2010). This study set three different exchangeable correlations, r = 0.2, 0.4, and 0.6. The three correlations chosen reflect a range of magnitudes of correlations among items that are often found in surveys, tests, and questionnaires in education.

Number of Items.
Three different numbers for items or trials were used, i.e., m = 8, 15, and 30. As this study tried to mimic real testing situations, those numbers were chosen based on hypothetical reliability coefficients calculated based on the Spearman-Brown formula (Spearman, 1910), which predicts a psychometric reliability for a lengthened test. All three conditions of 8, 15, and 30 resulted in predicted reliability around or above 0.8 (a satisfactory condition for reliability) for each of working correlation coefficients of 0.2, 0.4, and 0.6. An exception was made for the case of the working correlation 0.2 with m of 8, which resulted in 0.6 that is considered an acceptable cutoff.
Sample Size. Two different sample or cluster sizes were chosen, N = 50 and 200. The two conditions were chosen based on Pan's recommendation (2001b) using varying conditions of the GEE simulations with exchangeable correlations of 0.2, 0.4, and 0.6. The study showed that when N = 50, in most cases the power was less than or around 0.8 in most of the experimental conditions, and N = 200 always led to the power above 0.9.

Data Generation
The use of the GEE approach to generate correlated binary responses was advocated by Prentice (1988) because fully parametric approaches to estimating marginal response probabilities from a series of binary responses could be intractable except paired binary cases. Emrich and Piedmonte (1991) developed an algorithm reflecting Prentice's suggestion. The first step in the algorithm is to obtain correlated binary responses by applying threshold approaches to correlated continuous responses satisfying a desired marginal logit model. Each set of correlated binary observations or Bernoulli responses is treated like a block fulfilling the marginal condition. Thus, the 0 or 1 pattern within a block does not necessarily need to be the same as the patterns in other blocks as long as it meets the marginal condition, indicating the binary response patterns can vary across different blocks or individuals. Then, the dependence structure within the responses is parameterized in term of the given correlation matrix of the latent continuous responses (Emrich & Piedmonte, 1991;Touloumis, 2016).
An R package SimCorMultRes (Touloumis, 2017) adopting the suggested steps was used to simulate data sets with a multilevel structure (e.g., items at level 1 and subjects at level 2 in this study) given varying exchangeable correlations, sample sizes, and item numbers as mentioned before. The simulated multilevel data sets were analyzed using GEE. In the second step, all of the simulated data sets were made level-1 aggregated or clustered, which is the same as the process of summing up the binary responses to produce test total scores in practice. Two predictors included in data generation were X1~ N(0, 1) and X2 ~ Po(3). Both of the two covariates X1 and X2 were the level-2 predictors because the beta-binomial model cannot accommodate any level-1 predictor as the model utilizes the aggregate outcomes only and the two models can be only comparable under the same conditions.
For each type of distribution (i.e., positively and negatively skewed), a set of true parameters for the coefficients (β0, β1, β2) were submitted. For the negatively skewed distributions, the three true parameters were 0.7, 0.3, and 0.2, respectively. In the positively skewed distributions, each parameter was fixed at -0.7, 0.3, and -0.2 correspondingly. Those true values and the distributions of the two covariates reflected results of an empirical analysis performed with an educational test. For GEE modeling, the R package geepack (Højsgaard et al., 2016) and another package, aod (Lesnoff & Lancelot, 2012) for beta-binomial modeling with maximum likelihood estimation, were used.

Evaluation Criteria
Performance of GEE and the beta-binomial model was evaluated using six evaluation criteria. They reported different ways to appreciate accuracy and THE BETA-BINOMIAL MODEL FOR CLUSTERED BINARY RESPONSES 10 consistency of estimates and inferences from the two models. All of the calculations under each experimental condition were based on 1000 replications. To evaluate the parameter estimates, the bias and mean squared error (MSE) were reported. In addition, the coverage probability (CP) of 95% confidence intervals for each model was calculated. To evaluate inference results, Type 1 error and power rates were reported. Empirical Type 1 error rates were computed as the ratio that a test statistic produced a p-value less than 0.05 under the null hypothesis. Power was calculated as the ratio of a p-value below 0.05 where each coefficient was different from 0. Additionally, agreement rates of the inferential decisions on each coefficient at the nominal significance level of 0.05 were reported.

Errors in the GEE and the Beta-Binomial Estimates
This section reports the errors in the parameter estimates of the two models. The bias, mean squared error (MSE), and coverage probability (CP) of 95% confidence intervals for the coefficient estimates are reported in Table 1. As the positively and negatively skewed shapes did not seem to influence the sizes and patterns of the errors, the MSE values averaged over the two types of shapes are displayed in Figure 1, depending on the degree of intracluster correlations, the number of items, and sample size.

13
The MSE values of the beta-binomial estimates were found to be smaller than those of the GEE estimates across all of the experimental conditions. Regardless of the models, the higher correlations resulted in the higher MSE values. The higher the sample size, the lower the MSE values. The added number of items had smallest effect on the MSE values. Table 1 shows that the coverage probabilities demonstrated almost nominal 95% coverage in both models across most of the experimental conditions, although they seem to be slightly liberal. Recalling that β1 = 0.3 for both L and J shapes, and β2 = 0.2 for J-shapes and -0.2 for L-shapes, signs of biases as shown in Table 1 indicated that the GEE coefficients were overestimated, while the beta-binomial coefficients were underestimated in most of the experimental conditions (see also Table 2).

The GEE and Beta-Binomial Estimates
The parameter estimates of GEE and the beta-binomial model are presented in Table 2. Figure 2 summarizes the estimates under different conditions of sample size and intracluster correlation. When the sample size was small, the beta-binomial estimates were closer to the true parameters than the GEE estimates in most of the experimental conditions. In contrast, when the sample size got larger, it was noteworthy that the GEE estimates were consistent, staying close to the true parameters even with the highest correlation of 0.6, but the beta-binomial estimates were diverged with the added correlations. The differences between the estimates of the two models increased with the added correlations. Figure 3 compares the density plots of the GEE and beta-binomial estimates along with histograms of the aggregate outcomes in three experimental conditions that were chosen to contrast the small, medium, and large mean differences in the parameter estimates of the two models. One should note that the ranges of the estimates at the bottom plot were much larger than those shown in the upper two plots, although they do not look much different at first glance. In the uppermost plot where the mean differences were small, the density lines of the GEE and betabinomial estimates look very close. In the plot at the bottom with the large differences between the estimates, the lines were quite farther apart. The differences grew larger with higher intracluster correlations.    The last two columns under each model in Table 2 show that the estimates of intracluster correlation stayed very close across all of the experimental conditions. The sizes of the estimated intracluster correlations were slightly different, but the differences were fringe. For the simulated data sets with r = 0.2, the averaged correlation estimates were 0.099 for GEE and 0.100 for the beta-binomial model. When r = 0.4, the estimates were 0.214 and 0.217, respectively. Under r = 0.6, the estimates were 0.352 and 0.360 for each model correspondingly. The smaller estimates than the true correlations used for the data generation indicated that the two hypothetical predictors explained some amount of overdispersion due to dependent responses within clusters. In a follow-up analysis, the estimated correlations from the intercept only model without the two predictors were much closer to the true parameters, though they were about 0.04 to 0.1 lower than the true values on average for both.

Inference on the Beta-Binomial and GEE Estimates
The power and empirical Type 1 error rates and percent agreement of statistical decision (AD) on the estimates are reported in Table 3. The Type 1 error rates in both models ranged from 0.031 to 0.075 and were slightly larger than the nominal 0.05 level of significance on average. But they remained fairly stable across all of the conditions of correlations, sample sizes, and numbers of items.
The two models agreed upon rejection or acceptance for the coefficients at α = 0.05. On average, the agreement rate was about 95% for the large sample and dropped to 88% for the small sample. The inference on the correlation estimates (not included in Table 3) showed higher agreement, 92% for the small sample condition and 99.98% for the large sample condition on average. The power rates of the GEE and beta-binomial estimates were averaged over the shapes and the number of items in Figure 4. Notably, the power of the estimates drastically changed from the small sample condition to the large sample condition. The higher intracluster correlations diminished the power of the estimates.

Discussion and Conclusions
This study examined performance of the beta-binomial model for aggregate outcomes from clustered binary responses leading to non-normal marginal distributions. Its performance was compared to the results of GEE. The simulated data sets used for each model had different structures. In GEE, the binary observations signifying correct or incorrect responses to each simulated item were the units of analysis, and the binary responses were nested within subjects. In the beta-binomial model, the units of analysis were the aggregate outcomes of correlated binary responses. This study found that GEE tended to overestimate the regression parameters, consistent with previous findings from the Hallgren et al. (2016) study. GEE under the large sample size performed consistently and was only slightly affected by the different extents of intracluster correlation. However, under small sample size, the estimates were prone to bias with higher correlations. Though it is known that GEE is less dependent on the sample size (Kenward et al., 1994;Muth et al., 2015), some researchers (e.g., Gunsolley et al., 1995) expressed concern about the performance of GEE when used for small sample. In the current study, GEE estimates under small sample or cluster size were found to be vulnerable for a particular case of binomial data with high intracluster correlation. In such a case, use of a modified variance estimator (Pan, 2001a) or different tests such as the robust score test (Guo et al., 2005) can be preferable.
It was found that the beta-binomial model can be a good alternative for GEE when sample size is small. Its estimates were closer to the true parameters than those of GEE. Indeed, the strength of the beta-binomial model would be valuable because many researchers and practitioners employ small sample or cluster size conditions when collecting data sets for their pilot studies or in experimental design settings. The beta-binomial model could be readily applied when sample size of clustered binary data is not large enough.
Under the large sample size condition, the beta-binomial model generated more biased estimates. The consistency of the beta-binomial estimates was questionable especially with high correlations, opposite to the results of GEE that showed stable performance with the large sample. Hence, GEE and the betabinomial model can complement each other.
This study also found that in both small and large sample size conditions, the beta-binomial model was more efficient than GEE. The absolute parameter estimates of the beta-binomial model were underestimated throughout experimental conditions, which was consistent with the result of Harrison (2015). The two models agreed upon rejection or acceptance for the coefficients at the significance level of 0.05. The power of the parameter estimates was low for the small sample, yet very high at around 0.9 for the large sample size, which was similar to the result of Pan's (2001b) study. The Type 1 error probabilities of the two models were very close across the different conditions, confirming the empirical finding of Dilba and Aerts (2004). Though the values were slightly larger than the nominal significance level of 0.05, they remained fairly stable across the conditions in both models. Among the three experimental factors, the number of items or trials had the smallest effect on the parameter estimates, akin to the finding regarding the beta-binomial estimates in Harrison (2015). In the current study, the extent of intracluster correlation had profound effects on the estimation as well as on the power of the estimates for both models as found in Galbraith et al. (2010).
As this study examined cases mimicking educational tests where aggregate or clustered outcomes have non-normal and skewed distributions from correlated binary responses, generalizability of the conclusions within a broader context of aggregate outcomes might be limited. To overcome this limitation, a future study can include other aggregate or clustered designs such as a cross-sectional design where the intraclass correlation could be much smaller. Findings of this study warrant further investigation on the effect of varying degrees of intracluster correlation under more refined sample size conditions. Different data generation processes can be considered especially to further investigate the biased results of the two models under different sample size conditions. One can also ponder a systematic study including other statistical approaches that take into account overdispersion differently such as GLMM. This might answer the reason for somewhat underestimated correlation estimates even when the two predictors were absent in both models because other possible source of overdispersion such as interindividual variability could not be captured by the two models.