Homogeneity test of proportions for combined unilateral and bilateral data via GEE and MLE approaches

In comparative clinical trials, binary bilateral observations often consist of paired organ (e.g., eyes, ears, and kidneys) records in dichotomous form, with '1' denoting a cured or affected status, and '0' otherwise. In practice, such datasets often comprise a mixture of bilateral and unilateral observations, since some subjects contribute paired organ records, while others only a single record, for reasons such as prior surgical removal, congenital absence, or missing measurement for one organ. It is therefore essential to analyze the combined bilateral and unilateral data together, rather than discarding unilateral observations, to avoid loss of valuable information and potential bias in inference. Furthermore, appropriate modeling of the intra-subject correlation inherent in paired bilateral outcomes is necessary to accurately capture the dependence structure and ensure valid statistical conclusions. Rosner^[1] first proposed the 'constant R model' to test whether the proportion of affected eyes is equal across all groups, assuming equal dependence between the two eyes of the same individual. Dallal^[2] critiqued Rosner's model for its poor fit when the binary status (e.g., cured or affected) is present bilaterally, yet shows substantially different prevalences across groups. To address this, he proposed a model where the conditional probability of a response in one organ, given a response in the other, is a constant. Donner^[3] introduced an alternative approach assuming a constant intra-subject correlation across all individuals, which Thompson^[4] subsequently demonstrated to be robust through simulation studies.

Recent literature (2015−2025) demonstrates that these three models remain central to modern analyses of bilateral and combined uniliateral-bilateral data. Independent research groups continue to apply Rosner's, Dallal's, and Donner's methods to a variety of clinical contexts, including ophthalmology, otolaryngology, and other paired organ investigations. Several recent studies employ Rosner's model for the analysis of correlated bilateral outcomes^[5−11], while others rely on Dallal's model in applications involving correlated bilateral data^[12−16], or Donner's constant correlation framework to evaluate dependence in bilateral measurements^[17−20]. These works collectively reaffirm that bilateral and combined unilateral-bilateral data continue to represent an active statistical research area, and that classical dependence models are still widely used in modern biomedical statistics.

Testing the homogeneity of proportions addresses a fundamental inferential question in comparative studies, namely whether multiple groups share a common marginal probability of response while accounting for within-subject dependence in correlated outcomes. Methods based on maximum likelihood estimates (MLEs), and likelihood-based tests under the foregoing models have been developed in numerous subsequent studies. For example, Tang et al.^[21] proposed an asymptotic method for testing the equality of proportions between two groups under Rosner's model for correlated binary data. Ma et al.^[5], and Ma & Liu^[17] investigated three tests (likelihood ratio, Wald-type, and score) for testing homogeneity among g ≥ 2 groups under Rosner's and Donner's models, respectively. Ma & Wang K.^[7], and Ma & Wang H.^[19] further examined several test procedures to test the homogeneity of general g ≥ 2 proportions for combined unilateral and bilateral data under Rosner's and Donner's models, respectively.

Alternatively, regression models for correlated paired organ data can be used to perform hypothesis testing^[22,23], with the generalized linear mixed model (GLMM)^[24−27], and generalized estimating equations (GEE)^[28−30] being the most common approaches. Both methods are implemented in standard statistical software such as SAS, R, and Stata. The GEE method, introduced by Liang & Zeger for longitudinal data^[28−30], and later extended to clustered data^[31], generalizes the generalized liner model (GLM) framework to account for within-cluster dependence via estimating equations. It has been widely adopted for longitudinal and clustered data, and is often preferred over GLMM because it yields consistent parameter estimates and robust (co)variance estimates even when the 'working' correlation structure is misspecified^[32−35].

Bilateral observations from the same subject can be regarded as repeated measurements on that individual, making the GEE framework a natural choice for analyzing such correlated binary data. In this paper, the GEE method is employed to analyze combined unilateral and bilateral data, and compare its performance with established likelihood-based procedures, including the likelihood ratio, Wald-type, and score tests. Specifically, the homogeneity of proportions for combined unilateral and bilateral data under both Rosner's and Donner's models are investigated. A recent study by Zhang & Ma^[36] compared the score test with the GEE method for the equal proportion test in the combined data framework under Rosner's model. The present work extends this line of research by conducting a systematic comparison between the GEE approach, and the three likelihood-based methods under both Rosner's and Donner's models, providing a broader evaluation in this context.

The rest of the paper is organized as follows. The Methods section presents the MLE approaches using likelihood ratio, Wald-type and score statistics, as well as the GEE method for correlated binary outcomes and its implementation in SAS. The Results section reports the numerical results, including a simulation study evaluating the empirical type I error rates and powers of these methods under different models, and two real-world applications in otolaryngologic and ophthalmologic studies. The paper concludes with a discussion of the findings and their implications.

Consider a study involving the combined unilateral and bilateral date, where in the i-th group (i = 1, …, g), there are m_+i subjects who contribute data from both paired organs (bilateral), and n_+i subjects who contribute data from one of the paired organs (unilateral), respectively. Let m_ri be the number of bilateral subjects who have r (r = 0, 1, 2) organs cured or affected, and n_ri be the number of unilateral subjects who have r (r = 0, 1) organs cured or affected, such that:

where, Z_ijk denotes the response (1 cured or affected; 0 otherwise) of the k-th paired organ (k = 1, 2) of the j-th subject in the i-th group. The data structure is summarized in Table 1, where a subscript '+' indicates the summation over the corresponding index.

The proportion of cured or affected organs in the i-th group is assumed to be Pr (Z_ijk = 1) = π_i. Given m_+i and n_+i in the i-th group, (m_0i, m_1i, m_2i) follows trinomial distribution and (n_0i, n_1i) follows binomial distribution, i.e.,

where, the joint probabilities p_ri's (r = 0, 1, 2) read

with Corr (Z_ij1, Z_ij2) being the intra-subject correlation between the two responses from the j-th subject in the i-th group.

In what follows, two parametric models proposed by Rosner^[1] and Donner^[3] are considered to address the intra-subject correlation for the binary data. Under Rosner's model, the conditional probability is specified as Pr (Z_ijk = 1 | Z_ij,_3−k = 1) = Rπ_i, where R is a scalar parameter. Based on this specification, the intra-subject correlation becomes Corr (Z_ij1, Z_ij2) = (R − 1)π_i/(1 − π_i). Under Donner's model, a constant correlation ρ is assumed across all g groups, i.e., Corr (Z_ij1, Z_ij2) = ρ. Under either model, the joint probabilities in Eq. (2) can be written in terms of π_i and the nuisance parameter κ (κ = R under Rosner's and κ = ρ under Donner's model).

MLE approach

Let β = (π₁, …, π_g, κ)^T be the vector of parameters. For given observation:

the log-likelihood function reads:

where, the term 'const' denotes a constant depending on (m, n).

Our interest is the homogeneity test of proportions across the g groups. Thus, the hypotheses are:

The MLEs $ \hat{ {\boldsymbol{\beta}}}_0 $ under H₀ can be solved analytically, while the MLEs $ \hat{ {\boldsymbol{\beta}}}_1 $ under H₁ have no closed-form solutions and must be computed using an iterative method. Detailed procedures for obtaining $ \hat{ {\boldsymbol{\beta}}}_0 $ and $ \hat{ {\boldsymbol{\beta}}}_1 $ can be found in the works of Ma & Wang K.^[7] and Ma & Wang H.^[19], corresponding to Rosner's and Donner's models, respectively. Based on these MLEs, three likelihood-based test statistics are considered and defined as follows:

Likelihood Ratio (LR) test:

Wald-type test:

where, C^T is a (g − 1) × (g + 1) hypothesis matrix with (C^T)_ii = 1, (C^T)_i,i+1 = −1 for i = 1, …, g − 1, and (C^T)_ij = 0 otherwise; I⁻¹ (β) is the inverse of the (g + 1) × (g + 1) Fisher's information matrix whose explicit forms under Rosner's and Donner's models are given in Appendices in the works of Ma & Wang K.^[7], and Ma & Wang H.^[19], respectively.

Score test:

where, U = (∂l/∂π₁, …, ∂l/∂π_g, 0) is the score function, and I⁻¹ (β) is the inverse of Fisher's information matrix I (β).

It can be shown that the above three test statistics asymptotically follow the chi-square distribution with degrees of freedom df = g −1 under H₀, i.e., Q_LR, Q_W, Q_S $ \stackrel{d}{\to}\chi^2_{g-1} $.

GEE approach

Treating the bilateral observations as repeated measurements on the same subject, the GEE method can be employed to perform an equal proportion test with the mean model $ \mathbb{E}$(Z)_ij = $\mu $_ij = π_i${\rm 1\mskip-4.5mu l}_{k_{ij}} $ ($ {\rm 1\mskip-4.5mu l}_{k_{ij}} $ denotes a k_ij × 1 column vector of ones), and the estimating equations S(β) = $\sum_{j=1}^{n_i} $∂${\mu}^{T}_{ij} $/∂π_iV_ij⁻¹(Z_ij − ${\mu}_{ij} $) = 0, where Z_ij denotes the response(s) of the j-th subject in the i-th group, which is a 2 × 1 vector for bilateral data and a scalar for unilateral data, k_ij = 2 for bilateral data and k_ij = 1 for unilateral data, and V_ij = ${{A}}_{ij}^{1/2} {{R}}_{ij} {A}_{ij}^{1/2} $ is the "working covariance" matrix. Here, R_ij is the "working correlation" and A_ij is the diagonal matrix of marginal variances of Z_ij. The working correlation matrix R_ij generally depends on unknown parameters that must be estimated. In the case of paired organ data, R_ij contains only a single parameter. Various working correlation structures can be specified, such as independent, exchangeable, and unstructured correlations. Estimates for the parameters of both the mean model and the working covariance are obtained through an iterative procedure based on the estimating equations (see, e.g., the GENMOD procedure in the SAS user's guide^[37]).

Several statistical software systems can implement the GEE approach for analyzing combined unilateral and bilateral data. In this study, the SAS GENMOD procedure is used. Table 2 illustrates the structure of the combined unilateral and bilateral data for the case g = 2 in a single simulation.

An example of SAS code invoking the GENMOD procedure for g = 2 is presented in Table 3. In this example, the link function is the identity link specified with link = identity, and the unstructured working correlation matrix is specified with type = un. In addition, the by replicate statement repeats the analysis for each simulation case, and the contrast statement enables pairwise evaluation of differences, corresponding to the equal proportion test. Note that with repeated statement, the default statistic provided by the contrast statement is a generalized score statistic Q_GS computed based on the generalized score function S (β)^[37]. A generalized Wald-type statistic can be obtained with the wald option in the contrast statement. However, the generalized Wald-type statistic is known to exhibit poorer finite-sample performance and less accurate control of nominal significance levels compared with the generalized score statistic^[37,38].

A simulation study was performed to assess the performance of the three likelihood-based test statistics along with the GEE-based generalized score test, by investigating the empirical type I error and powers, respectively, under Rosner's and Donner's model. Equal and unequal sample size for g = 2, 4, 8 are considered. Specifically, for equal sample size design, we set:

and for unequal sample size design, we set:

Empirical type I error

To compute the empirical type I error rates, the dataset is generated according to the distributions in Eq. (1) with the above sample size designs under H₀ : π₁ = … = π_g = π₀ with a pre-specified π₀ = 0.2, 0.5. Additionally, the intra-subject correlation is set to be ρ₀ = 0.1, 0.5, 0.7. Thus, in data generation, R = R₀ = (1 − π₀)ρ₀/π₀ + 1 is set under Rosner's model, and ρ = ρ₀ under Donner's model.

After generating the dataset with each sample size design and parameter configuration, the three likelihood-based test statistics Q_LR, Q_W, Q_S in Eqs. (5)−(7), and the GEE-based generalized score test statistic denoted by Q_GS are computed accordingly. The null hypothesis H₀ : π₁ = … = π_g is rejected if Test Statistic > $ \chi_{1-\alpha,\; g-1}^2 $, where $ \chi_{1-\alpha,\; g-1}^2 $ is the quantile function of $ \chi_{g-1}^2 $ valued at 1 − α. The above simulation is replicated for N = 100,000 times. Then the empirical type I error rate is calculated as:

where, $ Q_i^{H_0} $ denotes the test statistic resulting from the dataset generated under H₀, and the subscript 'i' represents the type of tests for i = LR, W, S, GS.

The results for the empirical type I error with each sample size design and parameter configuration are presented in Table 4 under Rosner's model, and in Table 5 under Donner's model, respectively, where the liberal results are highlighted in boldface (The results are classified as liberal if $ \widehat{{\rm{TIE}}} $ > 0.06, as robust if 0.04 ≤ $ \widehat{{\rm{TIE}}}$ ≤ 0.06, and as conservative if $ \widehat{{\rm{TIE}}}$ < 0.04). For simplicity, the designs of the sample size (m₊₁, …, m_+g), and (n₊₁, …, n_+g) are referred to as (m, n) in Tables 4 and 5. The two equal sample size designs, where m₊₁ = … = m_+g = n₊₁ = … = n_+g = 20, 40, are denoted as E₁ and E₂, respectively. The unequal sample size designs are labeled as U_g for g = 2, 4, 8. For instance, U₂ refers to the design where (m₊₁, m₊₂) = (n₊₁, n₊₂) = (20, 40). Note that these notations also apply to the tables presented for power calculation.

Under Rosner's model, as shown in Table 4, when g = 2, the LR and Wald-type tests produce a few liberal results for the sample size design E₁, E₂ and E₁, U₂, respectively. As g = 4, the number of liberal results from the Wald-type tests grows substantially across all sample size designs (E₁, E₂ and U₄), whereas the LR tests show a similar number of the liberal results as in the g = 2 case. When g = 8, the majority of empirical type I error rates from the Wald-type tests are liberal, with several being considerably inflated (rate > 0.10) across all sample size designs. Overall, as the number of groups increases, the LR tests become mildly liberal, while the Wald-type tests exhibit a dramatic increase in both the frequency and amplitude of the liberal behavior. In addition, no clear association is observed between the empirical type I error rates, and either the null proportion π₀ or the intra-subject correlation ρ₀ for the LR or Wald-type tests. In contrast, the score tests consistently yield non-liberal empirical type I error rates across all simulation settings. The GEE tests show robust performance for all sample size designs when g = 2 and g = 4, and only a few mildly liberal results when g = 8. Overall, the results from the score tests closely align with those from the GEE approach and remain closer to the nominal level, as compared to the LR and Wald-type tests.

Under Donner's model, as shown in Table 5, the Wald-type tests begin to produce liberal results when g = 4 for sample size designs E₁, E₂, and U₄. As g = 8, most empirical type I error rates from the Wald-type tests are liberal across all designs. Compared with the results under Rosner's model, however, the liberal behavior under Donner's model is less severe, with empirical rates generally remaining below approximately 0.08. In contrast, the LR tests yield only a couple of liberal results as g = 2 for sample size design E₁ and U₂, when both π₀ and ρ₀ are small. The Score and GEE tests exhibit behavior similar to that observed under Rosner's model.

It should be noted that when the sample size is small (e.g., (m, n) = E₁) and both π₀ and ρ₀ are low, some score tests tend to be conservative, whereas the GEE tests may exhibit mild liberal behavior, particularly under Rosner's model. These deviations are likely attributable to sparse cell counts in the resulting contingency tables, which can undermine the validity of asymptotic approximations. As the sample size increases (e.g., (m, n) = E₂), however, both the score and GEE tests exhibit improved and stable type I error control.

Powers

The powers are calculated in a similar way as the empirical type I error. The difference is that instead of generating the dataset under H₀, they are generated under certain alternative hypothesis. There are two alternatives, H_1A and H_1B used for power calculation, which are shown below.

for g = 2, 4, 8. Tables 6 and 7 present the results for powers under under H_1A and H_1B within Rosner's (with R = 1.2, 1.5, 1.8), and Donner's (with ρ = 0.1, 0.5, 0.7) model, respectively.

It can be seen that the powers under H_1A are generally lower than those under H_1B. This is because the discrepancy in proportions specified under H_1A is smaller compared to that under H_1B. Within each power table, the Wald and LR tests are generally the most and second most powerful tests, respectively, for all the sample size designs and parameter configurations. The score tests and GEE tests are less powerful, but their associated powers are comparable to those of Wald and LR tests. In addition, for each given number of groups g and sample design, the powers obtained by the GEE tests decrease as the R (under Rosner's model) or ρ (under Donner's model) increases. The same pattern is observed in the other three likelihood-based tests under Donner's model.

Therefore, based on the simulation results, we conclude that the likelihood-based score test and the GEE test are robust under both Rosner's and Donner's models, with their resulting type I error rates under control and powers that are generally expected, though slightly underestimated. This echoes the findings for the score test in the work of Ma & Wang K.^[7] and Ma & Wang H.^[19], and in the meantime, recommends the GEE as an alternative approach.

It should be noted that the similar empirical performance observed between the likelihood-based score test Q_S and the GEE-based generalized score test Q_GS does not imply that the MLE and GEE approaches are theoretically equivalent. Rather, this similarity arises because both procedures are score-type tests targeting the same null hypothesis and are asymptotically equivalent under correctly specified marginal mean models. In the simulation settings considered here, the estimating equations underlying the GEE generalized score test closely resemble the score equations derived from the likelihood under the assumed dependence structures, leading to comparable finite-sample behavior. Nevertheless, important distinctions remain between the two approaches. The MLE-based score test relies on a fully specified likelihood and explicit modeling assumptions for the within-subject dependence, whereas the GEE-based generalized score test is quasi-likelihood-based and retains validity even when the working correlation structure is misspecified. Therefore, the observed numerical similarity in these simulations should be interpreted as a consequence of the specific study design rather than as evidence of general equivalence between GEE and likelihood-based methods.

Real world example

Two real-world examples are studied to illustrate the equal proportion test. Since there is more than one model available for analysis, a model selection process was first conducted to identify a more suitable model to describe the dataset, following the goodness-of-fit test procedure^[39].

The first example consists of a subset of 214 children who were admitted because of acute otitis media with effusion (OME), and were randomized into two groups, respectively treated with cefaclor and amoxicillin^[40]. Table 8 shows the number of cured ears in 173 children at 42 d.

Based on the goodness-of-fit test procedure by Zhou & Ma^[39], it shows that both Rosner's and Donner's models are suitable for this dataset (all p-values for goodness-of-fit test $ \gtrsim $ 0.9), with Rosner's model yielding slightly lower AIC (274.1305 under Rosner's model vs 274.1406 under Donner's model). Therefore, Rosner's model was selected to perform the equal proportion test. The constrained and unconstrained MLEs, along with the statistics and p-values under Rosner's model, can be found in Table 9. It can be seen that all the p-values from the three likelihood-based tests, and from the GEE test are of similar magnitude ($ \gtrsim $ 0.8), and much larger than the nominal level of α = 0.05. Therefore, we fail to reject the null hypothesis H₀ : π₁ = π₂, indicating no significant difference in the proportion of the cured ears between the two treatments (cefaclor versus amoxicillin). This conclusion aligns with the findings in the original study^[40], which reported that 'by 42 d after entry, the percentage of children whose ears were without effusion or "improved" was equal in both treatment groups (68.9% in the cefaclor group and 67.5% in the amoxicillin group)', despite overlooking the intra-subject correlation. However, by taking the intra-subject correlation into account in the likelihood-based tests and the GEE test, stronger evidence for the equal proportion of cured ears across treatment groups is observed, as reflected by very large p-values.

The second example involves the ophthalmologic study conducted in the Massachusetts Eye and Ear Infirmary between 1970 and 1979, where data were collected from an outpatient population of patients with retinitis pigmentosa (RP), and their normal relatives^[41]. Patients were classified into four types of genetic groups: (1) autosomal dominant RP (DOM); (2) autosomal recessive RP (AR); (3) sex-linked RP (SL); and (4) isolate RP (ISO). A subset of 218 patients aged 20−39 was originally analyzed by Rosner using the constant R model^[1], where an eye was considered affected if the best corrected Snellen visual acuity (VA) was 20/50 or worse, and normal if VA was 20/40 or better. Table 10 presents the distribution of the number of affected eyes for 216 patients in the four genetic groups who had complete information for VA on both eyes. This example can be considered as a special case in the combined data framework where unilateral observations are absent.

Using the goodness-of-fit test results presented in Zhou & Ma^[39], it can be seen that Donner's model is much superior to Rosner's model as reflected by one magnitude larger p-values, and smaller AIC (443.7967 under Donner's model vs 449.9490 under Rosner's model). Therefore, Donner's model was selected to perform the equal proportion test. The constrained and unconstrained MLEs, along with the statistics and p-values under Donner's model, are presented in Table 11. As can be seen, all the p-values are smaller than 0.05, with the p-values from the score and GEE tests being close to each other. Therefore, the null hypothesis H₀: π₁ = π₂ = π₃ = π₄ is rejected, and it is concluded that there is an overall difference in proportions of the affected eyes across the four genetic groups. This finding is consistent with the original paper by Rosner^[1], which also noted that the overall difference was completely attributed to the differences between the SL group, and each of the other three groups.

In the present study, the GEE approach is considered as an alternative to three likelihood-based methods (likelihood ratio, Wald-type, and score statistics) for testing homogeneity of proportions across g ≥ 2 groups for the combined unilateral and bilateral data. The likelihood-based methods are revisited under two statistical models: i) Rosner's constant R model; and ii) Donner's equal correlation ρ model, accounting for intra-subject correlations in the bilateral portion, and introduce the GEE method along with its implementation in SAS GENMOD procedure for the combined organ data.

Monte Carlo simulations evaluate the empirical type I error rates, and powers of each method under various sample sizes and parameter configurations. Results show that the GEE and score tests exhibit comparable performance, outperforming the likelihood ratio and Wald-type tests in controlling type I error. All tests demonstrate similar power, though the likelihood ratio, and Wald-type tests are slightly inflated. These findings confirm that the GEE test performs at least as well as the score test and revalidate previous results reported by Ma and others^[7,19]. Importantly, the GEE method offers additional flexibility, allowing incorporation of continuous covariates and more complex model structures, which provides an advantage over the score test.

Applications to two real datasets from otolaryngologic and ophthalmologic studies illustrate these methods. Goodness-of-fit tests guide model selection for the likelihood-based methods, showing that Rosner's model is preferred for the OME dataset, while Donner's model better fits the retinitis pigmentosa dataset. Inference from both the GEE method and the likelihood-based methods after model selection are consistent with the original study results.

In conclusion, the GEE and score tests are recommended for homogeneity testing for the combined unilateral and bilateral data. Both methods provide robust type I error control, and strong power across different scenarios. For extended analyses involving continuous covariates or multiple explanatory variables, the GEE approach is preferred. The score test is computationally efficient, particularly for large samples or numerous groups, and is well suited for exact tests, such as Fisher's exact or permutation tests, when data are sparse, which is less straightforward with standard software implementation of GEE.