Inference and sample size calculation for clinical trials with incomplete observations of paired binary outcomes

We introduce the following notations. Let Y it be a binary outcome (1 for Yes and 0 for No) from the i th ( i = 1, ··· , n ) subject under treatment t ( t = 0, 1). Here we generically use treatment to denote research categories: case/control, MR/EC, pre-/post-interventions, etc. For paired outcomes ( Y i 0 , Y i 1 ), we define joint probabilities h u v = P ( Y i 0 = u , Y i 1 = v ) for u = 0 , 1 and v = 0 , 1 . Then we have marginal probabilities p 0 = P ( Y i 0 = 1) = h 10 + h 11 and p 1 = P ( Y i 1 = 1) = h 01 + h 11 . The relationship between joint and marginal probabilities are illustrated in Table 1 . We also define ρ = Corr( Y i 0 , Y i 1 ) to be the within-subject correlation coefficient. It can be shown that (1) ρ = h 11 − p 0 p 1 p 0 ( 1 − p 0 ) p 1 ( 1 − p 1 ) . For binary variables, there is a range for the valid values of within-subject correlation, ρ ∈ [ ρ L , ρ U ], determined by the marginal probabilities: ρ L = max { − p 0 p 1 ( 1 − p 0 ) ( 1 − p 1 ) , − ( 1 − p 0 ) ( 1 − p 1 ) p 0 p 1 } , ρ U = min { p 0 ( 1 − p 1 ) p 1 ( 1 − p 0 ) , p 1 ( 1 − p 0 ) p 0 ( 1 − p 1 ) } . Details of the derivation can be found in Chaganty and Joe [ 15 ]. Equation (1) also suggests that there is a one-to-one correspondence between ( p 0 , p 1 , h 11 ) and ( p 0 , p 1 , ρ ), and they fully characterize the statistical properties of ( Y i 0 , Y i 1 ). In the rest of this paper we will use ( p 0 , p 1 , ρ ) and ( p 0 , p 1 , h 11 ) exchangeably in the description of paired outcomes. The treatment effect is denoted by Δ = p 1 – p 0 = h 01 – h 10 . We are interested in testing the null hypothesis H 0 : p 1 – p 0 = 0. Without loss of generality, we define n = m 1 + m 2 + m 3 and number the subjects (from 1 to n ) in such a fashion that the first m 1 subjects only contribute outcomes under treatment 0, the next m 2 subjects contribute outcomes under both treatments, and the last m 3 subjects only contribute outcomes under treatment 1. We use Y 1 = { Y i 0 , i = 1, ··· , m 1 }, Y 2 = {( Y i 0 , Y i 1 ), i = m 1 + 1, ··· , m 1 + m 2 }, and Y 3 = { Y i 1 , i = m 1 + m 2 + 1, ··· , n } to denote the collection of outcomes from the three groups. We also define r = m 1 / n and s = m 3 / n to be the proportions of subjects with partial outcomes. In the following discussion, we assume a missing completely at random mechanism. That is, the occurrence of missing values is independent of the observed or unobserved outcomes. Note that the setting described above is equivalent to that assumed by Tang et al. [ 16 ], which stated that paired observations ( Y 2 ) follow a multinomial distribution with probabilities ( h 00 , h 01 , h 10 , h 11 ), and incomplete observations ( Y 1 and Y 3 ) follow independent binomial distributions with probabilities p 0 and p 1 . The general idea of existing approaches to estimating the treatment effect proceeds as follows. 1) Using only complete pairs of observations from Y 2 to obtain estimators of p 0 and p 1 , denoted by p ^ 0 ( P ) and p ^ 1 ( P ) . Then use unpaired observations from Y 1 and Y 3 to obtain estimators p ^ 0 ( U ) and p ^ 1 ( U ) , respectively. Here superscripts ( P ) and ( U ) indicate whether the estimators are based on paired or unpaired observations. 2) Because p ^ t ( P ) and p ^ t ( U ) are both unbiased estimators of the true response rates p t ( t = 0, 1), a class of hybrid estimators can be constructed as the weighted average of p ^ t ( P ) and p ^ t ( U ) , p ^ t ( H ) = ϕ t p ^ t ( P ) + ( 1 − ϕ t ) p ^ t ( U ) , with weight 0 ≤ ϕ t ≤ 1. Here the superscript ( H ) indicates the hybrid nature of the estimator. 3) A class of unbiased estimator for treatment effect is Δ ^ 1 = p ^ 1 ( H ) − p ^ 0 ( H ) . Specifically, the aforementioned estimators of p 0 and p 1 are p ^ t ( P ) = ∑ i = m 1 + 1 m 1 + m 2 Y i t ∕ m 2 based on Y 2 for t = 0 , 1 ; p ^ 0 ( U ) = ∑ i = 1 m 1 Y i 0 ∕ m 1 based on Y 1 ; p ^ 1 ( U ) = ∑ i = m 1 + m 2 + 1 n Y i 1 ∕ m 3 based on Y 3 . A traditionally adopted weighting scheme [ 6 , 7 , 8 ] is to set ϕ 0 = m 2 /( m 1 + m 2 ) and ϕ 1 = m 2 /( m 2 + m 3 ). Under this weighting scheme, the hybrid estimator can be simplified as (2) p ^ 0 ( H ∗ ) = ∑ i = 1 m 1 + m 2 Y i 0 m 1 + m 2 , p ^ 1 ( H ∗ ) = ∑ i = m 1 + 1 n Y i 1 m 2 + m 3 , and Δ ^ 1 ∗ = p ^ 1 ( H ∗ ) − p ^ 0 ( H ∗ ) . Under weights ϕ 0 = m 2 /( m 1 + m 2 ) and ϕ 1 = m 2 /( m 2 + m 3 ), n ( Δ ^ 1 ∗ − Δ ) approximately follows a normal distribution with mean 0 and variance (3) σ 1 2 = p 0 ( 1 − p 0 ) 1 − s + p 1 ( 1 − p 1 ) 1 − r − 2 ( 1 − r − s ) ( h 11 − p 0 p 1 ) ( 1 − s ) ( 1 − r ) . The expression of σ 1 2 has been investigated by many researchers [ 6 , 7 , 8 ]. For readers’ reference, we include the derivation details in Appendix A . The pre- and post-intervention outcomes ( Y i 0 , Y i 1 ) are observed as pairs in Y 2 . In constructing the hybrid estimator in (2) , however, they are used as if they were independently observed. i.e., Y i 0 and Y i 1 from Y 2 are used separately in p ^ 0 ( H ∗ ) and p ^ 1 ( H ∗ ) . As a result, although Δ ^ 1 ∗ remains an unbiased estimator, it might lose efficiency. In the following, we propose a different hybrid estimator which directly estimates the treatment effect Δ (instead of estimating p 0 and p 1 first and then taking the difference) based on paired and unpaired observations, denoted by Δ ^ ( P ) and Δ ^ ( U ) . The advantage of Δ ^ ( P ) is that it accounts for the pairing of observations in Y 2 . Then a new class of hybrid estimator can be defined as the weighted sum of Δ ^ ( P ) and Δ ^ ( U ) , Δ ^ 2 = w Δ ^ ( P ) + ( 1 − w ) Δ ^ ( U ) , with weight 0 ≤ w ≤ 1. We further investigate the optimal weighting scheme that minimizes the variance of Δ ^ 2 . (4) τ 1 2 = p 1 ( 1 − p 1 ) ∕ s + p 0 ( 1 − p 0 ) ∕ r . The estimator of the treatment effect based on independent observations ( Y 1 and Y 3 ) is Δ ^ ( U ) = p ^ 1 ( U ) − p ^ 0 ( U ) . The large sample theory suggests that n ( Δ ^ ( U ) − Δ ) is approximately normally distributed with mean 0 and variance The estimator of the treatment effect based on paired observations ( Y 2 ) is constructed using within-subject difference d i = Y i 1 − Y i 0 , Δ ^ ( P ) = ∑ i = m 1 + 1 m 1 + m 2 d i ∕ m 2 , which takes into account potential within-subject correlation. [ 4 ] has shown that Var( d i ) = h 01 + h 10 – ( h 01 – h 10 ) 2 . Thus the approximate distribution of n ( Δ ^ ( P ) − Δ ) is normal with mean 0 and variance (5) τ 2 2 = h 01 + h 10 − ( h 01 − h 10 ) 2 1 − r − s . Because both Δ ^ ( U ) and Δ ^ ( P ) are unbiased estimators of the treatment effect, we can define a new class of unbiased Because both estimators as Δ ^ 2 = w Δ ^ ( P ) + ( 1 − w ) Δ ^ ( U ) with 0 ≤ w ≤ 1. The approximate distribution of n ( Δ ^ 2 − Δ ) is normal with mean 0 and variance Var ( n Δ ^ 2 ) = w 2 τ 2 2 + ( 1 − w ) 2 τ 1 2 . The optimal weight, denoted by w *, which minimizes Var ( n Δ ^ 2 ) , is solved by setting ∂ Var ( n Δ ^ 2 ) ∕ ∂ w = 0 . It is easy to obtain solution w ∗ = r 1 2 ∕ ( τ 1 2 + τ 2 2 ) . Thus the new hybrid estimator is Δ ^ 2 ∗ = w ∗ Δ ^ ( P ) + ( 1 − w ∗ ) Δ ^ ( U ) . It combines inference from two data sources: unpaired outcomes from independent cohorts ( Y 1 , Y 3 ) and paired outcomes from a single cohort ( Y 2 ). It has the smallest variance among all weighted averages of Δ ^ ( P ) and Δ ^ ( U ) . Under the optimal weight w ∗ , n ( Δ ^ 2 ∗ − Δ ) is approximately normal with mean 0 and variance σ 2 2 = τ 1 2 τ 2 2 τ 1 2 + τ 2 2 = ( 1 ∕ τ 1 2 + 1 ∕ τ 2 2 ) − 1 . Note that σ 2 2 is the harmonic mean of τ 1 2 and τ 2 2 . It is always true that σ 2 2 ≤ min { τ 1 2 , τ 2 2 } . i.e., the hybrid estimator Δ ^ 2 ∗ is more efficient than the original estimators Δ ^ ( P ) and Δ ^ ( U ) . Furthermore, by rewriting w * as w ∗ = 1 ∕ τ 2 2 1 ∕ τ 1 2 + 1 ∕ τ 2 2 , it is clear that between Δ ^ ( P ) and Δ ^ ( U ) , the one with the better precision is assigned the larger weight in the hybrid estimator Δ ^ 2 ∗ . Calculating the optimal weight w * requires the values of ( p 0 , p 1 , h 11 ). In practice, their true values are usually unknown, and are replaced by observed sample proportions. Due to the independence between Y 1 and Y 2 , as well as between Y 2 and Y 3 , the weights in Δ ^ 1 ∗ , ϕ 0 = m 2 /( m 1 + m 2 ) and ϕ 1 = m 2 /( m 2 + m 3 ), are optimal in minimizing the variances of p ^ 0 ( H ) and p ^ 1 ( H ) , respectively. However, when the paired outcomes from the same subject are correlated, minimizing the individual variances of p ^ 0 ( H ) and p ^ 1 ( H ) does not guarantee minimizing the variance of their difference Δ ^ 1 = p ^ 1 ( H ) − p ^ 0 ( H ) . That is, Δ ^ 1 ∗ might not be optimal among the class of estimators Δ ^ 1 . The proposed estimator Δ ^ 2 ∗ is advantageous in that it accounts for the pairing of observations in Y 2 and the optimal weight w * directly minimizes the variance of the treatment effect estimator Δ ^ 2 = w Δ ^ ( P ) + ( 1 − w ) Δ ^ ( U ) . It should be pointed out, however, that the optimality of Δ ^ 2 ∗ is achieved within the class of estimators Δ ^ 2 , which is a subclass of Δ ^ 1 under the constraint ϕ 0 = ϕ 1 . Specifically, Δ ^ 2 = w Δ ^ ( P ) + ( 1 − w ) Δ ^ ( U ) = w ( p ^ 1 ( P ) − p ^ 0 ) ( P ) ) + ( 1 − w ) ( p ^ 1 ( U ) − p ^ 0 ( U ) ) = [ w p ^ 1 ( P ) + ( 1 − w ) p ^ 1 ( U ) ] − [ w p ^ 0 ( P ) + ( 1 − w ) p ^ 0 ( U ) ] . We have the following result concerning the relative efficiency between Δ ^ 1 ∗ and Δ ^ 2 ∗ : If r = s , then σ 2 2 ≤ σ 1 2 . The “=” sign holds only when ρ = 0. See Appendix B . Theorem 1 offers assurance that when the proportions of subjects with missing outcomes are balanced between the two treatments ( r = s ), the proposed hybrid estimator Δ ^ 2 ∗ will always be more efficient than the traditional estimator Δ ^ 1 ∗ . There is no definite conclusion where r ≠ s . In Figure 1 we plot the relative efficiency between Δ ^ 1 ∗ and Δ ^ 2 ∗ , defined as R = σ 1 2 ∕ σ 2 2 , against correlation coefficient ρ , under different combinations of ( p 0 ; p 1 ; r , s ). We declare Δ ^ 2 ∗ to be more efficient than Δ ^ 1 ∗ , if R > 1, and vice versa. In each graph, the title indicates the combination of ( p 0 ; p 1 ) under investigation. The horizontal axis shows the value of ρ . Note that the valid range of ρ are different depending on ( p 0 ; p 1 ) in each graph. The vertical axis represents the relative efficiency R . Equal efficiency R = 1 is marked by a horizontal dashed line. We explore four missing patterns ( r = 0.1, s = 0.1), ( r = 0.2, s = 0.2), ( r = 0.1, s = 0.3), and ( r = 0.3, s = 0.1), represented by colors black, red, blue, and, green, respectively. Figure 1 empirially confirms the conclusion of Theorem 1, where the black and red curves (both under r = s ) are always above the line of equal efficiency. When r ≠ s , however, we observe that Δ ^ 1 ∗ tends to be more efficient than Δ ^ 2 ∗ if the within-subject correlation is negative or weakly positive. Such scenarios are relatively rare in paired observations from pre-post studies. The proposed estimator Δ ^ 2 ∗ becomes superior when there is a moderate to strong positive correlation. Such scenarios are more frequely the case in practice. We further observe that although Δ ^ 1 ∗ might be more efficient than Δ ^ 2 ∗ , it is only by a small margin when it does. Overall all the scenarios we have explored, the smallest relative efficiency we have observed is R = 0.949. That is, Δ ^ 1 ∗ is at most 5% more efficient that Δ ^ 2 ∗ . On the other hand, the largest relative efficiency we observed is R = 2.55, at ( p 0 = 0.4, p 1 = 0.45, r = 0.3, r = 0.1) and ρ = 0.9. That is, Δ ^ 2 ∗ can be substantially more efficient than Δ ^ 1 ∗ when there is a strong positive correlation between paired outcomes. Finally, with the formulae of σ 1 2 and σ 2 2 provided, it is straightforward for practitioners to numerically examine the relative efficiency of Δ ^ 1 ∗ verse Δ ^ 2 ∗ based on assumed ( p 0 , p 1 , ρ , r , s ) at the design stage. Under the proposed estimator Δ ^ 2 ∗ , in order to test null hypothesis H 0 : p 0 – p 1 = 0, we calculate test statistics (6) Z = ∣ Δ ^ 2 ∗ ∣ σ 2 2 ∕ n = n ∣ w ∗ Δ ^ ( P ) + ( 1 − w ∗ ) Δ ^ ( U ) ∣ w ∗ 2 τ 2 2 + ( 1 − w ∗ ) 2 τ 1 2 , and reject H 0 if Z > z 1– α /2 . Here z 1– α /2 is the 100(1 – α )th percentile of the standard normal distribution. Note that based on ( Y 1 , Y 3 ) only, the test statistic Z 1 = Δ ^ ( U ) ∕ τ 1 2 ∕ n is widely used to compare proportions between two independent groups. Based on Y 2 only, the test statistic Z 2 = Δ ^ ( P ) ∕ τ 2 2 ∕ n can be shown to be equivalent to the McNemar's test statistic. Equation (6) presents Z as a hybrid test statistic obtained from Z 1 and Z 2 . As for sample size determination, under the alternative hypothesis H α : p 1 – p 0 = Δ, achieving a testing power of 1 – γ implies that (7) P ( Z > z 1 − α ∕ 2 ∣ H α ) = 1 − γ . Here γ is the type II error. From (7) we have a closed-form sample size solution n 2 = σ 2 2 ( z 1 − α ∕ 2 + z 1 − γ ) 2 Δ 2 = τ 1 2 τ 2 2 ( z 1 − α ∕ 2 + z 1 − γ ) 2 ( τ 1 2 + τ 2 2 ) Δ 2 , where τ 1 2 and τ 2 2 are presented in (4) and (5) . Through similar steps, we can derive a close-form sample size formula based on estimator Δ ^ 1 ∗ . We denote this sample size by n 1 = σ 1 2 ( z 1 − α ∕ 2 + z 1 − γ ) 2 Δ 2 , where the expression of σ 1 2 is (3) .

We perform simulation studies to assess the performance of the proposed sample size formula n 2 . Combinations of ( p 0 , p 1 , r , s ) similar to those in Figure 1 are explored. We investigate four non-negative values of within-subject correlation, ρ = 0, 0.1, 0.25, and 0.5, because in biomedical studies we seldom encounter a negative correlation and the above values all fall within the valid range of ρ under the various pairs of ( p 0 , p 1 ) considered. We set the nominal levels of power at 1 – γ = 0.8 and type I error at α = 0.05. For each combination of ( p 0 , p 1 , ρ , r , s ), we first calculate sample size n 2 . Then 10,000 data sets are generated, each contains n 2 pairs of binary outcomes following ( p 0 , p 1 , ρ ), denoted by {( Z i 0 , Z i 1 ), i = 1, ··· , n 2 . To impose missing values, we generate a random vector of indicators g = ( g 1 , ··· , g n 2 ), where g i are independent and follow a multinomial distribution, taking values 1, 2, and 3, with probabilities r , (1 – r – s ), and s , respectively. We set Z i 1 to be missing if g i = 1, and set Z i 0 to be missing if g i = 3. The final simulated dataset is defined by Y 1 = { Z i 0 : g i = 1}, Y 2 = {( Z i 0 , Z i 1 ) : g i = 2}, and Y 3 = { Z i 1 : g i = 1}. Then hypothesis testing based on (6) is conducted on data set. Empirical power and type I error are computed as the proportion of times H 0 is rejected under the alternative and null hypothesis, respectively. A similar simulation study is performed for n 1 . In Tables 2 and 3 we show the simulation results for n 1 and n 2 . In the last three columns, we present the calculated sample size together with empirical power 1 − γ ^ and empirical type I error α ^ under different design configurations. In the column named n 0 , we list the sample sizes based on the McNemar's test [ 4 ] where all subject contribute complete pairs of observations ( r = s = 0). Over the scenarios that we have explored, the calculated sample sizes have a wide range, from 120 to 2308, and the empirical powers and type I errors are generally close to their nominal levels, indicating good performance of the sample size formulas. We also observe that the strength of within-subject correlation has a profound impact on sample size requirement. For example, under marginal probabilities ( p 0 = 0.4, p 1 = 0.45) and missing pattern ( r = 0.1, s = 0.1), if we assume no within-subject correlation ( ρ = 0), the required sample size is 1701. If we assume a correlation of ρ = 0.5, the sample size drops to 901, a 47% reduction. Thus at the design stage of experiments with paired outcomes, it is critical to conduct literature review or preliminary analysis to evaluate within-subject correlation. Furthermore, the simulations demonstrate that utilizing incomplete data does lead to saving in sample size. For example, under ( p 0 = 0.25, p 1 = 0.35, ρ = 0.25), if we exclude incomplete data from analysis, then in order to obtain n 0 = 251 complete pairs of outcomes to achieve 80% power at 5% significance level, we need to enroll 314 (251/0.8) and 419 (251/0.6) unique subjects under missing patterns ( r = 0.1, s = 0.1) and ( r = 0.2, s = 0.2), respectively. The proposed sample size n 2 indicates that, if we appropriately incorporate incomplete data into statistical inference, 280 and 327 unique subjects are actually needed to achieve the desired power and type I error, both substantially smaller than those required when ignoring incomplete data. Similar saving in sample size is observed on n 1 . Finally, comparing Tables 2 and 3 , under balanced missing patterns, ( r = 0.1, s = 0.1) and ( r = 0.2, s = 0.2), the sample sizes ( n 2 ) in Table 2 are always smaller than their counter parts ( n 1 ) in Table 3 . Under the unbalanced missing pattern, ( r = 0.1, s = 0.3), the sample sizes in Table 3 can be smaller. However, when within-subject correlation is moderately strong ( ρ = 0.5), the sample sizes in Table 2 are all smaller than those in Table 3 . This observation empirically confirms the conclusion about relative efficiency in Theorem 1. Due to space limit, we omit the simulation results under missing pattern ( r = 0.3, s = 0.1).

Evans et al. [ 17 ] performed a cohort study on mothers’ mood through pregnancy and after childbirth. The Edinburgh postnatal depression scale was implemented on each woman at 18 and 32 weeks of pregnancy, as well as 8 weeks and 8 months postpartum. 13,799 women were identified as eligible, among whom 12,059 completed at least one of the four questionnaires and 9,028 completed all four. The authors performed analysis based on the 9,028 women with complete data. The outcome of primary interest is whether a woman had depressive disorder at the four stages, indicated by a score above 12. It was reported that the numbers (%) of women with depressive disorder are 1062 (11.8%) at 18 weeks of pregnancy, 1222 (13.5%) at 32 weeks of pregnancy, 821 (9.1%) at 8 weeks postpartum, and 738 (8.2%) at 8 months postpartum. The joint probabilities of depression at the four stages were not reported. Suppose we would like to conduct a similar study to determine whether there is a significant difference in the risk of depression for women between the stages at 18 weeks of pregnancy and 8 weeks postpartum. During these two stages women undergo dramatical changes both physically and psychologically. We assume that the true incidence rates of depression at these two stages are p 0 =11.8% and p 1 =9.1%, respectively. Because in the previous study 9028/12059=75% of subjects contributed complete data, we specify 1 – r – s = 75%. We assume the incomplete data to be evenly distributed between the two stages, thus r = s = 12.5%. To detect a difference of Δ = p 1 – p 0 = −2.7% with 90% power and 5% type I error, the sample sizes calculated based on n 2 3077, 2538, 1966, 1667, 1359, and 713, under within-subject correlation ρ =0, 0.2, 0.4, 0.5, 0.6, and 0.8, respectively. The sample sizes calculated based on n 1 under the same configurations are 3077, 2553, 2029, 1767, 1505, 981. For illustration purposes, in the above calculation we have assumed MCAR. In studies of mood it is likely that the probabilities of missing data might depend on observed or unobserved data. In such cases an additional model for the missing mechanism is needed to achieve valid inference, and simulation studies are needed to assess sample size requirement.

Paired experimental design is frequently adopted in scientific studies. In many cases, the outcome of primary interest is binary. For example, in radiology, researchers perform two imaging procedures on the same group of patients to compare the rates of correct diagnosis between procedures. In social and behavioral sciences, pre-post studies are conducted such that responses from the same group of subjects before and after an intervention are compared to evaluate the effectiveness of an intervention. In the above two examples, each experimental unit contributes a pair of outcome measurements. Paired outcome measurements can also arise from matched case-control studies or matched randomized trials, where subjects in two comparison groups are matched by a set of characteristics potentially related to the outcome. In practice, researchers conducting paired experiments often encounter the missing data problem due to patient refusal, lost to follow-up, equipment failure or unavailability, etc. In such scenarios researchers are presented with a “mixture” type of data: paired (dependent) outcomes from subjects who contribute complete data and unpaired (independent) outcomes from those who contribute incomplete data. For example, non-invasive diagnostic tests of ultrasonography (US), magnetic resonance (MR) and echocolonoscopy (EC) were evaluated for their ability to provide correct diagnosis of endometriosis, but a significant portion of patients have missing values for MR or EC due to unavailability of the corresponding equipments at the occasion of evaluation [ 1 ]. Missing data can also be caused by protocol changes during the trial. Campbell et al. [ 2 ] described a controlled clinical trial on the safety and effectiveness of a cardiovascular device against the standard medical therapy, where in the middle of the trial, a protocol change in the assessment of primary outcome leads to missing baseline measurements in 60% of patients. Sometimes the “mixture” type of data arises from the experimental design. In a pre-post survey on 1000 subjects to evaluate the impact of a 3-month radio campaign on people's attitude toward colorectal cancer screening. Expecting that 20% of the subjects will be lost to follow-up at the end of the campaign, researchers plan to include 200 new subjects in the post-campaign survey. Thus a total of 1200 unique subjects will be involved, 200 pre-campaign only, 200 post-campaign only, and 800 both. If every subject contributes a complete pair of observations, then the McNemar's test is a straightforward choice to detect difference. On the other hand, if the outcomes are obtained from two independent groups of subjects, then the Chi-square test will be appropriate. The sample size calculation for the McNemar's test or the Chi-square test has been well established [ 3 , 4 , 5 ]. Recently there has been some development in statistical inference based on paired binary outcomes with missing data. Based on the large sample theory, various estimators were proposed along with the confidence intervals for the intervention effect [ 6 , 7 , 8 ]. Shih [ 9 ] proposed maximum likelihood estimation and likelihood ratio test with incomplete pairs. They employed an EM algorithm in computation. Choi and Stablein [ 10 ] investigated the effect of non-random missing data mechanisms. Tang and Tang [ 11 ] proposed an exact-test approach that compares two paired proportions with incomplete data. Tang et al. [ 12 ] investigated exact and approximate unconditional confidence intervals for proportional difference in the presence of incomplete data. Li et al. [ 13 ] developed a hypothesis testing method for a simple ordering in incomplete contingency tables and Li et al. [ 14 ] constructed confidence intervals for rate ratio in matched-pair studies with incomplete data. Despite the extensive application of paired experimental design, there has been limited research on sample size determination for such studies accounting for the impact of missing data. The rest of the paper is arranged as follows. In Section 2 we review the existing methods and present a new test statistic which combines information from paired and unpaired observations. Furthermore, closed-form sample size formulas are presented for paired experimental design with incomplete binary outcomes. In Section 3 we conduct simulation studies to examine the performance of sample size formulas. In Section 4 we demonstrate the sample size method using a real application example. Finally, we discuss limitations and possible extensions in Section 5.

Paired binary outcomes with potentially incomplete data are frequently collected and analyzed in biomedical and social behavioral studies. In this paper we propose a new estimator for the treatment effect which appropriately accounts for the mixed nature of the data and directly minimizes variance for the estimator of interest. We theoretically prove that under balanced missing patterns ( r = s ), the proposed estimator is always more efficient than the traditional estimator. We demonstrate through simulation that under unbalanced missing patterns ( r ≠ s ) the proposed estimator will be superior when there is a moderate to strong positive correlation between paired outcomes. We also derive close-form sample size formulas to help researchers determine how many unique subjects need to be enrolled in this type of studies. The sample size formulas are shown to preserve the nominal levels of power and type I error under various design configurations. One limitation of the proposed estimator is that Δ ^ 2 is not well defined when outcomes are observed from all subjects under one of the treatments ( r = 0 or s = 0), which is not the case for the traditional estimator Δ ^ 1 . Both Δ ^ 1 and Δ ^ 2 are developed based on normal approximation under the large sample theory. Thus the performance of the estimators and the corresponding sample size formulas depends on the quality of normal approximation. It is widely known that the normal approximation is less satisfactory under extreme values of probabilities ( p 0 and p 1 close to 0 or 1) or small sample sizes. Here we need to check the numbers of observations for both paired and unpaired outcomes. When normal approximation is unsatisfactory, nonparametric exact methods should be implemented, which will be one of our future research topics. The proposed sample size formula is derived based on the MCAR assumption. When the MCAR assumption does not hold, an additional model accounting for the missing mechanism is required to achieve valid inference. Because the missing mechanism is usually unknown and unique in each study, it is unlikely to derive a universally applicable sample size formula for studies without the MCAR assumption. In such cases, sample size requirement should be evaluated using customized simulation schemes designed for particular studies, which is beyond the scope of this paper. It is generally expected, however, that a more complicated model requires a larger sample size to maintain desired power and type I error.