- Research Article
- Open access
- Published:

# Confidence intervals construction for difference of two means with incomplete correlated data

*BMC Medical Research Methodology*
**volume 16**, Article number: 31 (2016)

## Abstract

### Background

Incomplete data often arise in various clinical trials such as crossover trials, equivalence trials, and pre and post-test comparative studies. Various methods have been developed to construct confidence interval (CI) of risk difference or risk ratio for incomplete paired binary data. But, there is little works done on incomplete continuous correlated data. To this end, this manuscript aims to develop several approaches to construct CI of the difference of two means for incomplete continuous correlated data.

### Methods

Large sample method, hybrid method, simple Bootstrap-resampling method based on the maximum likelihood estimates (*B*
_{1}) and Ekbohm’s unbiased estimator (*B*
_{2}), and percentile Bootstrap-resampling method based on the maximum likelihood estimates (*B*
_{3}) and Ekbohm’s unbiased estimator (*B*
_{4}) are presented to construct CI of the difference of two means for incomplete continuous correlated data. Simulation studies are conducted to evaluate the performance of the proposed CIs in terms of empirical coverage probability, expected interval width, and mesial and distal non-coverage probabilities.

### Results

Empirical results show that the Bootstrap-resampling-based CIs *B*
_{1}, *B*
_{2}, *B*
_{4} behave satisfactorily for small to moderate sample sizes in the sense that their coverage probabilities could be well controlled around the pre-specified nominal confidence level and the ratio of their mesial non-coverage probabilities to the non-coverage probabilities could be well controlled in the interval [0.4, 0.6].

### Conclusions

If one would like a CI with the shortest interval width, the Bootstrap-resampling-based CIs *B*
_{1} is the optimal choice.

## Background

Incomplete data often arise in various research fields such as crossover trials, equivalence trials, and pre and post-test comparative studies. For instance, ([1] pp. 212) designed a crossover clinical trial to measure the onset of action of two doses of formoterol solution aerosol: 12 ug and 24 ug. In this study, twenty-four patients were randomly allocated in equal numbers to one of the six possible sequences of two treatments at a time. Each patient was received two aerosols at each of visits 2 and 4. After four weeks, researchers measured the forced expiratory volume of a second (FEV_{1}) indicators for twenty-four patients. Due to the fact that researches did not consider all possible combinations of three treatments (e.g., placebo, 12 ug and 24 ug aerosols), which indicates that the missing data mechanism is missing completely at random (MCAR) thus FEV_{1} was only observed for 7 patients under both treatments (e.g., 12 ug and 24 ug aerosols), 9 patients only for 12 ug aerosol, and 8 patients only for 24 ug aerosol. The resultant data are shown in Table 1, which consist of two parts: the complete observations and the incomplete observations.

For the above crossover clinical trial, our main interest is to test the equivalence between 12 ug and 24 ug formoterol solution aerosols with respect to the FEV_{1} value. To this end, we can construct a (1−*α*)100 % confidence interval for the difference of two FEV_{1} values. If the resultant confidence interval (CI) lies entirely in the interval (−*δ*
_{0},*δ*
_{0}) with *δ*
_{0}(>0) being some pre-specified clinical acceptable threshold, we thus could conclude the equivalence between two doses of formoterol solution aerosol at the *α* significance level. As a result, reliable CIs for the difference in the presence of incomplete data are necessary.

The problem of testing the equality and constructing CI for the difference of two correlated proportions in the presence of incomplete paired binary data has received considerable attention in past years. For example, ones can refer to [2–6] for the large sample method, and [7] for the corrected profile likelihood method. When sample size is small, [8] proposed the exact unconditional test procedure for testing equality of two correlated proportions with incomplete correlated data. Tang, Ling and Tian [9] developed the exact unconditional and approximate unconditional CIs for proportion difference in the presence of incomplete paired binary data. Lin et al. [10] presented a Bayesian method to test equality of two correlated proportions with incomplete correlated data. Li et al. [11] discussed the confidence interval construction for rate ratio in matched-pair studies with incomplete data. However, all the aforementioned methods were developed for incomplete paired binary data.

Statistical inference on the difference of two means with incomplete correlated data has received a limited attention. For example, [12] discussed the problem of testing the equality of two means with missing data on one response and recommended [13] statistic when the variances were not too different. Lin and Stivers [14] also gave a similar comparison. Lin and Stivers [15] and [12] suggested some test statistics for testing the equality of two means with incomplete data on both response. However, to our knowledge, little work has been done on CI construction for the difference of two means with incomplete correlated data under the MCAR assumption.

Inspired by [16–19], we develop several CIs for the difference of two means with incomplete correlated data under the MCAR assumption based on the large sample method, hybrid method and Bootstrap-resampling method. The presented Bootstrap-resampling CIs have not been considered in the literature related to missing observations.

The rest of this article is organized as follows. Several methods are presented to construct CIs for the difference of the two means with incomplete correlated data in Section “Methods”. Simulation studies and an example are conducted to evaluate the finite performance of the proposed CIs in terms of coverage probability, expected interval width, and mesial and distal non-coverage probabilities in Section “Results”. A brief discussion is given in Section “Discussion”. Some concluding remarks are given in Section “Conclusion”.

## Methods

Suppose that ** x**=(

*x*

_{1},

*x*

_{2})

^{′}is a 2×1 vector of random variables, and follows a distribution with mean

**and covariance matrix**

*μ***given by**

*Σ*respectively. Let {(*x*
_{1m
},*x*
_{2m
}):*m*=1,⋯,*n*} be *n* paired observations on *x*
_{1} and *x*
_{2}, \(\left \{x_{1,n+1}, \cdots, x_{1,n+n_{1}}\right \}\) be *n*
_{1} additional observations on *x*
_{1}, \(\left \{x_{2,n+1}, \cdots, x_{2,n+n_{2}}\right \}\) be *n*
_{2} additional observations on *x*
_{2}. Thus, there are *n*
_{1} missing observations on *x*
_{2}, and *n*
_{2} missing observations on *x*
_{1}. Without loss of generality, the data may be presented as follows:

where (*x*
_{1m
},*x*
_{2m
}) is referred to as a paired observation, while *x*
_{1,n+j
} and *x*
_{2,n+k
} are referred to as incomplete or unpaired observations. Similar to [20, 21], throughout this article, it is assumed that the missing data mechanism is MCAR (i.e., independent of treatment and outcome). Based on these observations, we here want to construct reliable explicit CIs for the difference of two means *δ*=*μ*
_{1}−*μ*
_{2} under MCAR assumption.

### Confidence interval based on the large sample method

To make a comparison with the following proposed methods, we assume that ** x** follows a bivariate normal distribution in this subsection. In this case, if only variable

*x*

_{1}or

*x*

_{2}is subject to missingness (i.e.,

*n*

_{1}=0 or

*n*

_{2}=0), one can obtain the closed forms of the maximum likelihood estimates (MLEs) of

**and**

*μ***[22]. However, there are no closed forms of the MLEs for**

*Σ***and**

*μ***when variables**

*Σ**x*

_{1}and

*x*

_{2}are simultaneously subject to missingness (i.e.,

*n*

_{1}≠0 and

*n*

_{2}≠0), though one can find the MLEs of

**and**

*μ***using an iterative algorithm [23]. To get the closed forms of MLEs for**

*Σ***and**

*μ***, [15] proposed the modified MLEs using a non-iterative procedure and provided several test statistics based on the obtained estimators of**

*Σ***and**

*μ***.**

*Σ*
*(i) Confidence interval based on Lin and Stivers’s test statistics*

Let \(\hat {\delta }=\hat {\mu }_{1}-\hat {\mu }_{2}\) be the MLE of *δ* under the bivariate normal assumption of ** x**. When

**is known, it follows from [15] that the MLE of**

*Σ**δ*is

and the asymptotic variance of \(\hat {\delta }\) can be expressed as

respectively, where \(\overline {x}_{1}^{(n)}=\frac {1}{n}\sum _{j=1}^{n}x_{1j}\), \(\overline {x}_{2}^{(n)}=\frac {1}{n}\sum _{j=1}^{n}x_{2j}\), \(\overline {x}_{1}^{(n_{1})}=\frac {1}{n_{1}}\sum _{j=1}^{n_{1}}x_{1,n+j}\), \(\overline {x}_{2}^{(n_{2})}=\frac {1}{n_{2}}\sum _{k=1}^{n_{2}}x_{2,n+k}\), *a*=*nh*(*n*+*n*
_{2}+*n*
_{1}
*β*
_{21}), *b*=*nh*(*n*+*n*
_{1}+*n*
_{2}
*β*
_{12}), *β*
_{21}=*ρσ*
_{2}/*σ*
_{1}, *β*
_{12}=*ρσ*
_{1}/*σ*
_{2}, *h*=1/{(*n*+*n*
_{1})(*n*+*n*
_{2})−*n*
_{1}
*n*
_{2}
*ρ*
^{2}}. An approximate 100(1−*α*) % CI of *δ* is given by \(\left (\hat {\delta }-\textit {z}_{\alpha /2}\sqrt {\text {Var}(\hat {\delta })}, \hat {\delta }+\textit {z}_{\alpha /2}\sqrt {\text {Var}(\hat {\delta })}\right)\), which is denoted as *T*
_{
w1}-CI.

Following [15], when ** Σ** is unknown, the statistic for testing

*H*

_{0}:

*δ*=

*δ*

_{0}versus

*H*

_{1}:

*δ*≠

*δ*

_{0}is given by

which is asymptotically distributed as *t*-distribution with *n* degrees of freedom under *H*
_{0}, where *V*
_{1}=[{*A*
^{2}/*n*+(1−*A*)^{2}/*n*
_{1} }*m*
_{1}+{ *B*
^{2}/*n*+(1−*B*)^{2}/*n*
_{2} } *m*
_{2}−2*ABm*
_{12}/*n*]/(*n*−1), *A*={*n*(*n*+*n*
_{2}+*n*
_{1}
*m*
_{12}/*m*
_{1}}/{ (*n*+*n*
_{1})(*n*+*n*
_{2})−*n*
_{1}
*n*
_{2}
*r*
^{2}}^{−1}, *B*={*n*(*n*+*n*
_{1}+*n*
_{2}
*m*
_{12}/ *m*
_{2}} /{ (*n* + *n*
_{1})(*n* + *n*
_{2})−*n*
_{1}
*n*
_{2}
*r*
^{2}}^{−1}, \(m_{1}=\sum _{j=1}^{n} \left (x_{1j}-\overline {x}_{1}^{(n)}\right)^{2}\), \(m_{2}=\sum _{j=1}^{n}\left (x_{2j}-\overline {x}_{2}^{(n)}\right)^{2}\), \(m_{12}=\sum _{j=1}^{n}\left (x_{1j}\,-\,\overline {x}_{1}^{(n)}\right)\left (x_{2j}-\overline {x}_{2}^{(n)}\right)\), \(r=m_{12}/\sqrt {m_{1}m_{2}}\). Therefore, the approximate 100(1−*α*) % CI on the basis of *T*
_{1} is given by (*L*, *U*), where \(L=A\left (\overline {x}_{1}^{(n)}-\overline {x}_{1}^{(n_{1})}\right)-B\left (\overline {x}_{2}^{(n)}-\overline {x}_{2}^{(n_{2})}\right)+\overline {x}_{1}^{(n_{1})}-\overline {x}_{2}^{(n_{2})}-t_{\alpha /2}(n)\sqrt {V_{1}}\), and \(U=A\left (\overline {x}_{1}^{(n)}-\overline {x}_{1}^{(n_{1})}\right)-B\left (\overline {x}_{2}^{(n)}-\overline {x}_{2}^{(n_{2})}\right)+\overline {x}_{1}^{(n_{1})}-\overline {x}_{2}^{(n_{2})}+t_{\alpha /2}(n)\sqrt {V_{1}}\), which is denoted as *T*
_{1}-CI.

Another test statistic defined by [15] for testing *H*
_{0}:*δ*=*δ*
_{0} versus *H*
_{1}:*δ*≠*δ*
_{0}, which is a generalization of [24] test statistic for two independent samples, is given by

which is asymptotically distributed as *t* distribution with degrees *ν* of freedom, where \(\bar {x}_{1}^{(n+n_{1})}=(n+n_{1})^{-1}\sum _{j=1}^{n+n_{1}}x_{1j}\), \(\bar {x}_{2}^{(n+n_{2})}=(n+n_{2})^{-1}\sum _{j=1}^{n+n_{2}}x_{2j}\), *h*
_{1}=*n*{(*n*+*n*
_{2})*m*
_{1}/(*n*+*n*
_{1})+(*n*+*n*
_{1})*m*
_{2}/(*n*+*n*
_{2})−2*m*
_{12}}/{(*n*−1)(*n*+*n*
_{1})(*n*+*n*
_{2})}, *h*
_{2}=*n*
_{1}
*b*
_{1}/{(*n*
_{1}−1)(*n*+*n*
_{1})^{2}}, *h*
_{3}=*n*
_{2}
*b*
_{2}/{(*n*
_{2}−1)(*n*+*n*
_{2})^{2}}, \(b_{1}=\sum _{j=n+1}^{n+n_{1}}\left (x_{1j}-\overline {x}_{1}^{(n_{1})}\right)^{2}\), \(b_{2}=\sum _{j=n+1}^{n+n_{2}}\left (x_{2j}-\overline {x}_{2}^{(n_{1})}\right)^{2}\), and \(\nu =\left (h_{1}+h_{2}+h_{3}\right)^{2}/\{{h_{1}^{2}}/(n-1)+{h_{2}^{2}}/(n_{1}-1)+{h_{3}^{2}}/(n_{2}-1)\}\). Therefore, the approximate 100(1−*α*) % CI of *δ* for statistic *T*
_{2} is denoted as *T*
_{2}-CI.

When *σ*
_{1}=*σ*
_{2}, it follows from [15] that the statistic for testing *H*
_{0}:*δ*=*δ*
_{0} versus *H*
_{1}:*δ*≠*δ*
_{0} can be expressed as

which is asymptotically distribution as *t*-distribution with degrees *n*+*n*
_{1}+*n*
_{2}−4 of freedom. Note that when *n*
_{2}>*n*
_{1}, *b*
_{1}+*c*
_{2} should be replaced by *b*
_{2}+*c*
_{1}. Thus, the approximate 100(1−*α*) % CI of *δ* for *T*
_{3} is denoted as *T*
_{3}-CI, where \(c_{1}=\sum _{j=1}^{n+n_{1}}\left (x_{1j}-{n+n_{1}}\sum _{j=1}^{n+n_{1}}x_{1j}\right)^{2}\), and \(c_{2}=\sum _{j=1}^{n+n_{2}}\left (x_{2j}-\frac {1}{n+n_{2}}\sum _{j=1}^{n+n_{2}}x_{2j}\right)^{2}\).

Also, [12] presented the similar but simpler test statistics for testing the mean difference *δ*=*μ*
_{1}−*μ*
_{2}, which are adopted to construct CIs of *δ* as follows.

*(ii) Confidence interval based on Ekbohm’s test statistics*

Following [12], an unbiased estimator of *δ* is given by \(\hat {\delta }=\bar {x}_{1}^{(n+n_{1})}-\bar {x}_{2}^{(n+n_{2})}\), and its variance is given by \(\text {Var}(\hat {\delta })= \text {Var}(\hat {\mu })=\left \{(n+n_{2}){\sigma _{1}^{2}}+(n+n_{2}){\sigma _{2}^{2}}-2n\rho \sigma _{1}\sigma _{2}\right \}/\left \{(n+n_{1})(n+n_{2})\right \}\). An approximate 100(1 − *α*) % CI of *δ* can be obtained by \(\left (\hat {\delta }-\textit {z}_{\alpha /2}\sqrt {\text {Var}(\hat {\delta })}, \hat {\delta }+\right.\)
\(\left.\textit {z}_{\alpha /2}\sqrt {\text {Var}(\hat {\delta })}\right),\) which is denoted as *T*
_{
w2}-CI.

When *σ*
_{1}=*σ*
_{2}, Ekbothm (1976) proposed the following statistic for testing *H*
_{0}: \(T_{4}=(\tilde {\delta }-\delta _{0})\sqrt {\!(n\,+\,\!n_{1})(n\,+\,\!n_{2})\,-\,n_{1}n_{2}\lambda ^{2}}/\)
\( \left \{\hat \sigma \!\sqrt {2n(1-\lambda)+(n_{1}+n_{2})(1-\lambda ^{2})}\right \}\), where \(\tilde {\delta }=\left [n\left (n+n_{2}+n_{1}\lambda \right)\overline {x}_{1}^{(n)}\!-n\left (n+n_{1}+n_{2}\lambda \right)\overline {x}_{2}^{(n)}+n_{1}\!\left \{n\,+\,n_{2}\!\left (1\,-\,\lambda ^{2}\right)-\!n\lambda \right \}\overline {x}_{1}^{(n_{1})}-n_{2}\left \{n\,+\,n_{1}\!\left (1\!\,-\,\!\lambda \right)^{2}\!\,-\,n\lambda \right \}\overline {x}_{2}^{(n_{2})}\right ]\!\big /\!\!\left \{(n\,+\,n_{1})(n+n_{2})\,-\,n_{1}n_{2}\lambda ^{2}\right \}\), \(\hat {\sigma }^{2}\,=\,\left \{m_{1}\,+\,m_{2}\,+\,(1+\lambda ^{2})(b_{1}\,+\,b_{2})\right \}/\left \{2(n-\!1)+\!(1\,+\,\lambda ^{2})(n_{1}\,+\,n_{2}\,-\,2)\right \}\), and *λ*=2*m*
_{12}/(*m*
_{1}+*m*
_{2}). Under *H*
_{0}, *T*
_{4} is asymptotically distributed as *t*-distribution with degrees *n* of freedom. Therefore, the approximate 100(1−*α*) % CI is denoted as *T*
_{4}-CI.

Following [12], when *σ*
_{1}=*σ*
_{2}, another statistic for testing *H*
_{0} can be expressed as \(T_{5}= \left (\bar {x}_{1}^{(n+n_{1})}\!-\bar {x}_{2}^{(n+n_{2})}\,-\,\delta _{0}\!\right)\!\sqrt {(n\,+\,n_{1})(n\,+\,n_{2})/(R_{1}\!\,+\,R_{2})}\), which is asymptotically distributed as *t* distribution with degrees *ν*
_{
σ
} of freedom under *H*
_{0}, where *R*
_{1} = *n*(*m*
_{1} + *m*
_{2} − 2*m*
_{12}) /(*n* − 1), *R*
_{2} =(*n*
_{1} + *n*
_{2})(*b*
_{1}+*b*
_{2})/(*n*
_{1}+*n*
_{2}−2), and \(\nu _{\sigma }=\left (R_{1}+R_{2}\right)^{2}\!\!\!~/\left \{{R_{1}^{2}}/(n+1)+{R_{2}^{2}}/(n_{1}+n_{2})\right \}-2\). Thus, an approximate 100(1−*α*) % CI of *δ* for *T*
_{5} is denoted as *T*
_{5}-CI.

### Confidence interval based on the generalized estimating equations(GEEs)

To relax the bivariate normality assumption of ** x**, the method of the generalized estimating equations (GEEs) with exchangeable working correlation structure (e.g., [25]) can be adopted to make statistical inference on

*δ*in the incomplete correlated data because the GEE approach have become one of the most widely used methods in dealing with correlated response data [26, 27]. Following [28], the GEEs with exchangeable working correlation structure can be used to estimate parameter vector

**; the so-called sandwich variance estimator can be used to consistently estimate the covariance matrix of**

*μ***; and the ML method under a bivariate normal assumption via available paired observations is used to estimate the correlation parameter. Thus, an approximate 100(1−**

*μ**α*) % CI of

*δ*based on GEE method is denoted as

*T*

_{ g }-CI.

### Confidence interval based on the hybrid method

When the distribution function of ** x** is unknown, a hybrid method is developed to construct CI of

*δ*in this subsection. We first introduce the general concept of hybrid method. Let

*θ*

_{1}and

*θ*

_{2}be two parameters of interest. Now our main interest is to construct a 100(1−

*α*) % two-sided CI (

*L*,

*U*) of

*θ*

_{1}−

*θ*

_{2}via hybrid method. Let \(\hat {\theta }_{1}\) and \(\hat {\theta }_{2}\) be two estimates of

*θ*

_{1}and

*θ*

_{2}, respectively; and let (

*l*

_{1},

*u*

_{1}) and (

*l*

_{2},

*u*

_{2}) denote two approximate 100(1−

*α*) % CIs for

*θ*

_{1}and

*θ*

_{2}, respectively. Under the dependent assumption on \(\hat \theta _{1}\) and \(\hat \theta _{2}\), it follows from the central limit theorem that the approximate two-sided 100(1−

*α*) % CI of

*θ*

_{1}−

*θ*

_{2}is given by (

*L*,

*U*), where

Because \(\text {Cov}(\hat \theta _{1}, \hat \theta _{2})=\text {corr}(\hat \theta _{1}, \hat \theta _{2})\left \{\text {Var}(\hat {\theta }_{1})\text {Var}(\hat {\theta }_{2})\right \}^{1/2}\), the lower limit *L* and the upper limit *U* can be rewritten as

respectively. Note that (*l*
_{1},*u*
_{1}) contains the plausible parameter values of *θ*
_{1}, and (*l*
_{2},*u*
_{2}) contains the plausible parameter values for *θ*
_{2}. Among these plausible values for *θ*
_{1} and *θ*
_{2}, the values closest to the minimum *L* and maximum *U* are respectively *l*
_{1}−*u*
_{2} and *u*
_{1}−*l*
_{2} in spirit of the score-type CI [29]. From the central limit theorem, the variance estimates can now be recovered from *θ*
_{1}=*l*
_{1} as \(\widehat {\text {Var}(\hat \theta _{1})}=(\hat \theta _{1}-l_{1})^{2}/z_{\alpha /2}^{2}\) and from *θ*
_{2}=*u*
_{2} as \(\widehat {\text {Var}(\hat \theta _{2})}=\left (u_{2}-\hat \theta _{2}\right)^{2}\!\!\!\big /z_{\alpha /2}^{2}\) for setting *L*. As a result, the lower limit *L* for *θ*
_{1}−*θ*
_{2} is

Similarly, we can obtain

To obtain the above presented approximate 100(1−*α*) % hybrid CI for *μ*
_{1}−*μ*
_{2}, one requires evaluating the (1−*α*) 100 % CIs of *θ*
_{1} = *μ*
_{1} (denoted as (*l*
_{1}, *u*
_{1})) and *θ*
_{2}=*μ*
_{2} (denoted as (*l*
_{2}, *u*
_{2})), and estimating the correlation coefficient \(\widehat {\text {corr}}(\hat {\theta }_{1}, \hat \theta _{2})\). For the former, following [19], we consider the following two methods for getting the confidence limits (*l*
_{1}, *u*
_{1}) and (*l*
_{2}, *u*
_{2}) of *θ*
_{1} and *θ*
_{2}.

(i) *The Wilson score method*

where *N*
_{
i
}=*n*+*n*
_{
i
} and \(\hat {\theta }_{i}=\frac {1}{N_{i}}\sum _{j=1}^{N_{i}}x_{ij}\) for *i*=1,2.

(ii)*The Agresti-coull method*

where *N*
_{
i
}=*n*+*n*
_{
i
} and \(\tilde {\theta }_{i}=\left (\sum _{j=1}^{N_{i}}x_{ij}+0.5z_{\alpha /2}^{2}\right)/\left (N_{i}+z_{\alpha /2}^{2}\right)\) for *i*=1,2.

To construct CI for *δ*=*μ*
_{1}−*μ*
_{2} via the above described hybrid method, we can simply set *θ*
_{1}=*μ*
_{1} and *θ*
_{2}=*μ*
_{2}. If *Σ* is known, the estimated correlation coefficient \(\widehat {\text {corr}}(\hat {\mu }_{1}, \hat {\mu }_{2})\) of \(\hat {\mu }_{1}\) and \(\hat {\mu }_{2}\) is given by \(\widehat {\text {corr}}(\hat {\mu }_{1}, \hat {\mu }_{2})=2n\rho /\sqrt {(n+n_{1})(n+n_{2})}\). If *Σ* is unknown, \(\widehat {\text {corr}}(\hat {\mu }_{1}, \hat {\mu }_{2})\) is given by \(\widehat {\text {corr}}(\hat {\mu }_{1}, \hat {\mu }_{2})=nr/\left \{(n+n_{1})(n+n_{2})-n_{1}n_{2}r^{2}\right \}\), where \(r=m_{12}/\sqrt {m_{1}m_{2}}\), \(m_{1}=\sum _{j=1}^{n}\left (x_{1j}-\overline {x}_{1}^{(n)}\right)^{2}\) and \(m_{2}=\sum _{j=1}^{n}\left (x_{2j}-\overline {x}_{2}^{(n)}\right)^{2}.\) Thus, using Eqs. (1) and (2) yields CIs of *δ*=*μ*
_{1}−*μ*
_{2}. When *l*
_{
i
} and *u*
_{
i
} are estimated by the Wilson score method, we denote the corresponding CI as *W*
_{
s
}-CI; when *l*
_{
i
} and *u*
_{
i
} are estimated by the Agresti-coull method, the corresponding CI is denoted as *W*
_{
a
}-CI.

#### Bootstrap-resampling-based confidence intervals

When the distribution of ** x** is known, one can obtain the approximate CIs of

*δ*based on the asymptotic distributions of the constructed test statistics under the null hypotheses

*H*

_{0}:

*δ*=

*δ*

_{0}. However, when the distribution of

**is unknown, the asymptotic distributions of the constructed test statistics may not be reliable, especially with small sample size. On the other hand, estimators of some nuisance parameters have not the closed-form solutions even if the approximate distribution is reliable, and they must be obtained by using some iterative algorithms, which are computationally intensive. In this case, the Bootstrap method is often adopted to construct CIs of parameter of interest. The Bootstrap CIs can be constructed via the following steps.**

*x*
*Step 1.* Given the paired observations and incomplete observations

we draw *n* paired observations \(\left \{(x_{1m}^{*},x_{2m}^{*}):\! m=1, \cdots,n\right \}\) with replacement from *n* paired observations {(*x*
_{11},*x*
_{21}),⋯,(*x*
_{1n
},*x*
_{2n
})}, generate *n*
_{1} observations \(\{x_{1,n+j}^{*}: j=1, \cdots, n_{1}\}\) with replacement from \(\left \{x_{1,n+1}, \cdots, x_{1,n+n_{1}}\right \}\), and sample *n*
_{2} observations \(\left \{x_{2,n+k}^{*}: k=1, \cdots, n_{2}\right \}\) with replacement from \(\left \{x_{2,n+1}, \cdots, x_{2,n+n_{1}}\right \}\). Thus, we obtain the following Bootstrap resampling sample

*Step 2.* For the above generated Bootstrap resampling sample \(D_{b}^{*}\), we first compute \(\hat {\mu }_{1}^{*}=(n+n_{1})^{-1}\sum _{j=1}^{n+n_{1}}x_{1j}^{*}\) and \(\hat {\mu }_{2}^{*}=(n+n_{2})^{-1}\sum _{j=1}^{n+n_{2}}x_{2j}^{*}\), and then calculate the estimated value \(\hat {\delta }^{*}\) of *δ* via \(\hat {\delta }^{*}=\hat {\mu }_{1}^{*}-\hat {\mu }_{2}^{*}\).

*Step 3.* Repeating the above steps 1 and 2 for a total of *G* times yields *G* Bootstrap estimates \(\left \{\hat {\delta }_{g}^{*}: g=1,2,\cdots,G\right \}\) of *δ*. Let \(\hat \delta _{(1)}^{*}<\hat \delta _{(2)}<\cdots <\hat \delta _{(G)}^{*}\) be the ordered values of \(\left \{\hat \delta _{g}^{*}: g=1,2,\cdots,G\right \}\).

*Step 4.* Based on the bootstrap estimates \(\left \{\hat {\delta }_{g}^{*}, g=1,2,\ldots,G\vphantom {\left \{\hat {\delta }_{g}^{*}, g=\right.}\right \}\), Bootstrap-resampling-based CIs for *δ* can be constructed as follows.

Generally, the standard error se \((\hat {\delta })\) of \(\hat {\delta }\) can be estimated by the sample standard deviation of the *G* replications, i.e., \(\hat {\text {se}}(\hat {\delta })=\sqrt {(G-1)^{-1}\sum _{g=1}^{G}\left (\hat {\delta }_{g}^{*}-\bar {\delta }_{B}^{*}\right)^{2}}\), where \(\bar \delta _{B}^{*}=\left (\hat {\delta }_{1}^{*}+\cdots +\hat {\delta }_{G}^{*}\right)/G\). If \(\left \{\hat {\delta }_{g}^{*}: g=1,\cdots, G\right \}\) is approximately normally distributed, an approximate 100(1−*α*) % Bootstrap CI for *δ* is given by \(\left (\hat {\delta }-z_{\alpha /2}\hat {\text {se}}(\hat {\delta }), \hat {\delta }+z_{\alpha /2}\hat {\text {se}}(\hat \delta)\right)\), where *z*
_{
α/2} is the upper *α*/2-percentile of the standard normal distribution, which is referred as the simple Bootstrap confidence interval. When \(\hat {\delta }=a\overline {x}_{1}^{(n)}+(1-a)\overline {x}_{1}^{(n_{1})}-b\overline {x}_{2}^{(n)}-(1-b)\overline {x}_{2}^{(n_{2})}\), the corresponding simple Bootstrap CI is denoted as *B*
_{1}. When \(\hat {\delta }=\bar {x}_{1}^{(n+n_{1})}-\bar {x}_{2}^{(n+n_{2})}\), the corresponding simple Bootstrap CI is denoted as *B*
_{2}.

Alternatively, if \(\left \{\hat {\delta }_{g}^{*}: g=1,\cdots, G\right \}\) is not normally distributed, it follows from ([16] p.132) that the approximate 100(1−*α*) % Bootstrap-resampling-based percentile CI for *δ* is \(\left (\hat \delta _{\left ([G\alpha /2]\right)}^{*},\hat {\delta }_{([G(1-\alpha /2)])}^{*}\right)\), where [ *a*] represents the integer part of *a*, which is referred as the percentile Bootstrap CI. When \(\hat \delta =a\overline {x}_{1}^{(n)}+(1-a)\overline {x}_{1}^{(n_{1})}-b\overline {x}_{2}^{(n)}-\left (1-b\right)\overline {x}_{2}^{(n_{2})}\), the corresponding percentile Bootstrap CI is denoted as *B*
_{3}. When \(\hat {\delta }=\bar {x}_{1}^{(n+n_{1})}-\bar {x}_{2}^{(n+n_{2})}\), the corresponding percentile Bootstrap CI is denoted as *B*
_{4}.

## Results

### Simulation studies

In this subsection, we investigate the finite performance of various CIs in terms of empirical coverage probability (ECP), empirical confidence widths (ECW), and distal and mesial non-coverage probabilities (DNP and MNP) in various parameter settings via Monte Carlo simulation studies. A summary of abbreviation for various confidence intervals is presented in Table 2.

In the first simulation study, we consider the following case that (*n*,*n*
_{1},*n*
_{2}) is set to be (5,2,2); *μ*
_{1}=0,1,2; *μ*
_{2}=0.25,1,1.5; *ρ*=−0.9,−0.5,−0.1,0,0.1,0.5,0.9; *δ*=*μ*
_{1}−*μ*
_{2}=−0.25,0,0.5; \({\sigma _{2}^{2}}=4\); \({\sigma _{1}^{2}}=1,8\) and *α*=0.05. For a given combination (*n*,*n*
_{1},*n*
_{2},*μ*
_{1},*μ*
_{2},*ρ*,*σ*
_{1},*σ*
_{2}), we generate *n*+*n*
_{1}+*n*
_{2} random samples of (*x*
_{1},*x*
_{2})^{′} from a bivariate normal distribution with ** μ**=(

*μ*

_{1},

*μ*

_{2})

^{′}and

Then, for the generated *n*+*n*
_{1}+*n*
_{2} random samples, the *n*
_{1} observations on *x*
_{2} are deleted randomly. For the remaining paired *n*+*n*
_{2} random samples, the *n*
_{2} observations on *x*
_{1} are deleted randomly. Thus, (*x*
_{1m
},*x*
_{2m
})^{′}(*m*=1,⋯,*n*) are *n* pairs observations on (*x*
_{1},*x*
_{2})^{′}; *x*
_{1,n+j
}(*j*=1,⋯,*n*
_{1}) are *n*
_{1} additional observations on *x*
_{1}; *x*
_{2,n+k
}(*k*=1,⋯,*n*
_{2}) are *n*
_{2} additional observations on *x*
_{2}. Based on the observation {(*x*
_{1j
},*x*
_{2j
}):*m*=1,⋯,*n*}, {*x*
_{1,n+j
}:*j*=1,⋯,*n*
_{1}}, {*x*
_{2,n+k
}:*k*=1,⋯,*n*
_{2}}, we can draw 5000 bootstrap resampling samples. Independently repeating the above process *M*=10000 times, we can compute their corresponding ECP, ECW, MNP and DNP values. The ECP, ECW, MNP and DNP are defined by

respectively, where \(I\{\delta \in \mathcal {A}\}\) is an indicator function, which is 1 if \(\delta \in \mathcal {A}\) and 0 otherwise. The ratio of the MNP to the non-coverage probability (NCP) is defined as

Results are presented in Tables 3, 4 and 5. Also, to investigate the performance of the proposed CIs under the assumption \({\sigma _{1}^{2}}={\sigma _{2}^{2}}=\sigma ^{2}\), we calculate the corresponding results for *T*
_{3}, *T*
_{4}, *T*
_{5}, hybrid CIs, Bootstrap-resampling-based CIs when *σ*
^{2}=4 and (*n*,*n*
_{1},*n*
_{2})=(5,5,2), which are given in Tables 9, 10 and 11.

Following [17, 30], an interval can be regarded as satisfactory if (i) its ECP is close to the pre-specified 95 % confidence level, (ii) it possesses shorter interval width, and (iii) its RNCP lies in the interval [0.4,0.6]; *too mesially located* if its RNCP is less than 0.4; and *too distally* if its RNCP is greater than 0.6.

In the second Monte Carlo simulation study, we assume that the random samples of bivariate variables *x*
_{1} and *x*
_{2} are generated from a bivariate *t*-distribution with five degrees of freedom, and mean ** μ** and scale parameter

*Σ*specified in the first simulation study. The corresponding results with (

*n*,

*n*

_{1},

*n*

_{2})=(5,5,5) are given in Tables 6, 7 and 8. Similarly, we calculate the corresponding results for

*T*

_{3},

*T*

_{4},

*T*

_{5}, hybrid CIs, Bootstrap-resampling-based CIs when

*σ*

^{2}=4 and (

*n*,

*n*

_{1},

*n*

_{2})=(5,5,2), which are given in Tables 9, 10 and 11.

To investigate powers for the proposed CIs, we calculated the power in both the first and second simulation study. The results are shown in Tables 12 and 13. There is very little power in both the first and second simulation study to exclude a difference of zero.

### Results of simulation studies

From Tables 3, 4, 5, 6, 7, 8, 9, 10, 11, 12 and 13, we have the following findings. First, when *Σ* is unknown, the CIs based on the the Bootstrap-resampling-based methods except for *B*
_{3} behave satisfactorily in the sense that their ECPs are close to the pre-specified confidence level 95 % (e.g., see Tables 3 and 6); the CI based on the Bootstrap-resampling-based method *B*
_{1} generally yielded shorter ECWs than others (e.g., see Tables 4 and 7); the CIs corresponding to bivariate *t*-distribution are generally wider than those corresponding to bivariate normal distribution; the ECWs decrease as the correlation coefficient *ρ* increases. Second, the RNCPs of all the considered CIs lie in the interval [0.4,0.6] (e.g., see Tables 5 and 8), which show that our derived CIs generally demonstrate symmetry. Third, when \({\sigma _{1}^{2}}={\sigma _{2}^{2}}\), the CIs based on statistics *T*
_{3}, *T*
_{4} and *T*
_{5} behave unsatisfactory (e.g., see Tables 9 and 10) because their corresponding ECPs are almost less than the pre-specified confidence level 95 %. Fourth, powers corresponding to *W*
_{
a
} and *B*
_{1} are larger than others (e.g., see Tables 12 and 13). From the above findings, we would recommend the usage of the Bootstrap-resampling-based CI (i.e., *B*
_{1}) because its coverage probability is generally close to the pre-chosen confidence level, it consistently yields the shortest interval width even when sample size is small, it usually guarantees its ratios of the MNCPs to the non-coverage probabilities lying in [0.4, 0.6], and its power is usually larger than others.

### An worked example

In this subsection, the data introduced in Section for the action of two doses of formoterol solution aerosol are used to illustrate the proposed methodologies. In this example, we are interested in CI construction of the difference of two FEV_{1} values for two doses of formoterol solution aerosol. Under the previously given notation, we have *n*=7, *n*
_{1}=9, *n*
_{2}=8, \(\hat \delta =a\overline {x}_{1}^{(n)}+\left (1-a\right)\overline {x}_{1}^{(n_{1})}-b\overline {x}_{2}^{(n)}-(1-b)\overline {x}_{2}^{(n_{2})}=-0.0840\) (or \(\hat {\delta }=\sum _{j=1}^{n+n_{1}}x_{1j}/(n+n_{1})-\sum _{j=1}^{n+n_{2}}x_{2j}/(n+n_{2})=0.0228\)). Various 95 % CIs for *δ* under *Σ* unknown assumption are presented in Table 14. Examination of Table 14 shows that the actions of two doses of formaterol solutions aerosol are the same because all the derived CIs include zero.

## Discussion

Although testing equivalence of two correlated means with incomplete data has been studied, there is little work done on their interval estimators. To address the issue, this paper proposes various interval estimators of the difference of two correlated means for *Σ* known and unknown cases based on the large sample method, hybrid method and Bootstrap-resampling method. Extensive simulation studies are conducted to evaluate the finite performance of the proposed CIs in terms of the empirical coverage probability, empirical interval width and ratio of the mesial non-coverage probability to the non-coverage probability (RNCP). Empirical results evidence that the Bootstrap-resampling-based CIs *B*
_{1}, *B*
_{2}, *B*
_{4} behave satisfactorily for small to moderate sample sizes in the sense that their coverage probabilities could be well controlled around the pre-specified nominal confidence level and their RNCPs almost lie in the interval [0.4, 0.6]. However, confidence intervals based on the large sample method and hybrid method behave unsatisfactory for small sample sizes because the distributions of statistics *T*
_{1},⋯,*T*
_{5} are asymptotical, and these asymptotical distributions are proper only when *N*
_{
i
}→*∞*. When *Σ* is unknown, using GEE method to estimate variance is less efficient.

It is interesting to investigate confidence interval construction of the difference of two means with incomplete correlated data under missing at random and non-ignorable missing data mechanism assumptions of bivariate variables. We are working on the topics.

## Conclusion

According to the aforementioned findings, we can draw the following conclusions. The Bootstrap-resampling-based CI *B*
_{1} is a desirable interval estimator for the difference of two means with incomplete correlated data.

## References

Senn S. Cross-over Trials in Clinical Research, 2nd Edition. New York: John Wiley & Sons; 2002.

Choi SC, Stablein DM. Practical tests for comparing two proportions with incomplete data. Stat Med. 1982; 7:929–39.

Ekbohm G. On testing the equality of proportions in the paired case with incomplete data. Psychometrika. 1982; 47:115–8.

Campbell G. Testing equality of proportions with incomplete correlated data. J Stat Planning Inference. 1984; 10:311–21.

Bhoj DS, Snijders TAB. Testing equality of correlated proportions with incomplete data on both response. Psychometrika. 1986; 51:579–88.

Thomson PC. A hybrid paired and unpaired analysis for the comparison of proportions. Stat Med. 1995; 14:1463–70.

Pradhan V, Menon S, Das U. Correlated profile likelihood confidence interval for binomial paired incomplete data. Pharmaceutical Stat. 2013; 12:48–58.

Tang ML, Tang NS. Exact tests for comparing two paired proportions with incomplete data. Biometrical Journal. 2004; 46:72–82.

Tang ML, Ling MH, Tian G. Exact and approximate unconditional confidence intervals for proportion difference in the presence of incomplete data. Stat Med. 2009; 28:625–41.

Lin Y, Lipsitz S, Sinha D, Gawande AA, Regenbogen SE, Greenberg CC. Using Bayesian p-values in a 2×2 table of matched pairs with incompletely classified data. J R Stat Soc Series. 2009; C58:237–46.

Li HQ, Chan Ivan SF, Tang ML, Tian GL, Tang NS. Confidence interval construction for rate ratio in matched-paired studies with incomplete data. J Biopharmaceutical Stat. 2014; 24(3):546–68.

Ekbohm G. On comparing means in the paired cased with incomplete data on both responses. Biometrika. 1976; 63:299–304.

Morrison DF. A test for equality of means of correlated variates with missing data on one response. Biometrika. 1973; 60:101–5.

Lin PE, Stivers LE. Testing for equality of means with incomplete data on one variable: A Monte Carlo study. J Am Stat Assoc. 1975; 70:190–3.

Lin PE, Stivers LE. On difference of means with incomplete data. Biometrika. 1974; 61:325–s34.

Shao J, Tu DS. The Jackknife and Bootstrap. New York: Springer-Verlag; 1995.

Newcombe RG. Interval estimation for the difference between independent proportions: comparison of eleven methods. Stat Med. 1998; 17:873–90.

Zou GY, Donner A. Construction of confidence limits about effect measures: a general approach. Stat Med. 2008; 27:1693–702.

Tang NS, Luo XG. Confidence interval construction for sensitivity difference of two continuous-scale diagnostic tests at the fixed level of two specificities. Stat Probab Lett. 2015; 97:32–40.

Rubin DB. Inference and missing data. Biometrika. 1976; 63:581–92.

Seaman S, Galati G, Jackson D, Carlin J. What Is Meant by “Missing at Random”?Statistical Science. 2013; 28:257–268.

Anderson TW. Maximum likelihood estimates for a multivariate normal distribution when some observations are missing. J Am Stat Assoc. 1957; 52:200–4.

Orchard T, Woodbury MA. A missing information principle: theory and applications. Proc 6th Berkeley Symp. 1970; 1:697–715.

Welch BL. The generalization of ’Student’s’ problem when several different population variances are involved. Biometrika. 1947; 34:28–35.

Tsou TS. Exchangeable cluster binary data correlation coefficient estimation with generalized estimating equations. Stat Probab Lett. 2000; 50:179–86.

Liang KY, Zeger SL. Longitudinal data analysis using generalized linear models. Biometrika. 1986; 73:13–22.

Prentice RL. Correlated binary regression with covariates specific to each binary observation. Biometrics. 1988; 44:1033–48.

Pan W, Wall MM. Small-sample adjustments in using the sandwich variance estimator in generalized estimating equations. Stat Med. 2002; 21:1429–41.

Bartlett MS. Approximate confidence interval more than one unknown parameters. Biometrika. 1953; 40:306–17.

Tang ML, Li HQ, Tang NS. Confidence interval construction for proportion ratio in paired studies based on hybrid method. Stat Methods Med Res. 2010; 21:361–78.

## Acknowledgements

The research of Hui-Qiong LI was supported by the Natural Science Foundation of China (11201412, 11561075). The work of the second author was partially supported by the grants from the National Science Foundation of China (11225103).

## Author information

### Authors and Affiliations

### Corresponding author

## Additional information

### Competing interests

The authors declare that they have no competing interests.

### Authors’ contributions

HQL carried out the study, performed the statistical analysis and drafted the manuscript. NST participated in the design of the study, developed methods and revised the manuscript. YJY interpreted results and revised the manuscript. All authors commented on successive drafts, and read and approved the final manuscript.

## Rights and permissions

**Open Access** This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

## About this article

### Cite this article

Li, HQ., Tang, NS. & Yi, JY. Confidence intervals construction for difference of two means with incomplete correlated data.
*BMC Med Res Methodol* **16**, 31 (2016). https://doi.org/10.1186/s12874-016-0125-3

Received:

Accepted:

Published:

DOI: https://doi.org/10.1186/s12874-016-0125-3