ExactCIdiff: An R Package for Computing Exact Confidence Intervals for the Difference of Two Proportions

Comparing two proportions through the difference is a basic problem in statistics and has applications in many fields. More than twenty confidence intervals (Newcombe, 1998a,b) have been proposed. Most of them are approximate intervals with an asymptotic infimum coverage probability much less than the nominal level. In addition, large sample may be costly in practice. So exact optimal confidence intervals become critical for drawing valid statistical inference with accuracy and precision. Recently, Wang (2010, 2012) derived the exact smallest (optimal) one-sided 1− α confidence intervals for the difference of two paired or independent proportions. His intervals, however, are computer-intensive by nature. In this article, we provide an R package ExactCIdiff to implement the intervals when the sample size is not large. This would be the first available package in R to calculate the exact confidence intervals for the difference of proportions. Exact two-sided 1− α interval can be easily obtained by taking the intersection of the lower and upper one-sided 1− α/2 intervals. Readers may jump to Examples 1 and 2 to obtain these intervals.


Introduction
The comparison of two proportions through the difference is one of the basic statistical problems. Onesided confidence intervals are of interest if the goal of a study is to show superiority (or inferiority), e.g., that a treatment is better than the control. If both limits are of interest, then two-sided intervals are needed.
In practice, most available intervals, see Newcombe (1998a,b), are approximate ones, i.e., the probability that the interval includes the difference of two proportions, the so-called coverage probability, is not always at least the nominal level although the interval aims at it. Also, even with a large sample size, the infimum coverage probability may still be much less than the nominal level and does not converge to this quantity. In fact, the Wald type interval has an infimum coverage probability zero for any sample sizes and any nominal level 1 − α even though it is based on asymptotic normality, as pointed out by Agresti and Coull (1998) and Brown et al. (2001). See Wang and Zhang (in press) for more examples. Therefore, people may question of using large samples when such approximate intervals are employed since they cannot guarantee a correct coverage.
Exact intervals which assure an infimum coverage probability of at least 1 − α do not have this problem. But they are typically computer-intensive by nature. In this paper, a new R package ExactCIdiff (Shan and Wang, 2013) is presented which implements the computation of such intervals as proposed in Wang (2010Wang ( , 2012. The package is available from CRAN at http://CRAN.R-project. org/package=ExactCIdiff/. This package contains two main functions: PairedCI() and BinomCI(), where PairedCI() is for calculating lower one-sided, upper one-sided and two-sided confidence intervals for the difference of two paired proportions and BinomCI() is for the difference of two independent proportions when the sample size is small to medium. Results from ExactCIdiff are compared with those from the function ci.pd() in the R package Epi (Carstensen et al., 2013), and the PROC FREQ procedure in the software SAS (SAS Institute Inc., 2011).
Depending on how the data are collected, one group of three intervals is needed for the difference of two paired proportions and another group for the difference of two independent proportions. Pointed out by Mehrotra et al. (2003), an exact inference procedure may result in poor powerful analysis if an impropriate statistic is employed. Wang's one-sided intervals (Wang, 2010(Wang, , 2012, obtained through a carefully inductive construction on an order, are optimal in the sense that they are a subset of any other one-sided 1 − α intervals that preserve the same order, and are called the smallest intervals. See more details in the paragraph following (6). From the mathematical point of view, his intervals are not nested, see Lloyd and Kabaila (2010); on the other hand, for three commonly used confidence levels, 0.99, 0.95, 0.9, the intervals are nested based on our numerical study.
Although R provides exact confidence intervals for one proportion, e.g., the function exactci() in the package PropCIs (Scherer, 2013), the function binom.exact() in the package exactci (Fay, 2010(Fay, , 2012 and the function binom.test() in the package stats (Version 2.15.2), there is no exact confidence interval available in R, to the best of our knowledge, for the difference of two proportions, which is widely used in practice. ExactCIdiff is the first available R package to serve this purpose. The R package ExactNumCI (Sun and Park, 2013)   confidence interval for the difference of two independent proportions, however, pointed out by a referee, the coverage probability of a 95% confidence interval, when the numbers of trials in two independent binomial experiments are 3 and 4, respectively, is equal to 0.8734 when the two true proportions are equal to 0.3 and 0.5, respectively.
In the rest of the article, we discuss how to compute intervals for the difference of two paired proportions θ P defined in (1), then describe the results for the difference of two independent proportions θ I given in (7).

Intervals for the difference of two paired proportions
Suppose there are n independent and identical trials in an experiment, and each trial is inspected by two criteria 1 and 2. By criterion i, each trial is classified as S i (success) or F i (failure) for i = 1, 2. The numbers of trials with outcomes (S 1 , S 2 ), (S 1 , F 2 ), (F 1 , S 2 ) and (F 1 , F 2 ) are the observations, and are denoted by N 11 , N 12 , N 21 and N 22 , respectively. Thus X = (N 11 , N 12 , N 21 ) follows a multinomial distribution with probabilities p 11 , p 12 , p 21 , respectively. Let p i = P(S i ) be the two paired proportions. The involved quantities are displayed in Table 1. The parameter of interest is the difference of p 1 and p 2 : (1) To make interval construction simpler, let T = N 11 + N 22 and p T = p 11 + p 22 . We consider intervals for θ P of form [L(N 12 , T), U(N 12 , T)], where (N 12 , T) also follows a multinomial distribution with probabilities p 12 and p T . The simplified sample space is S P = {(n 12 , t) : 0 ≤ n 12 + t ≤ n} with a reduced parameter space H P = {(θ P , p T ) : p T ∈ D(θ P ), −1 ≤ θ P ≤ 1}, where D(θ P ) = {p T : 0 ≤ p T ≤ 1 − |θ P |}. The probability mass function of (N 12 , T) in terms of θ P and p T is p P (n 12 , t; θ P , p T ) = n! n 12 !t!n 21 ! p n 12 12 p t T p n 21 21 .
Suppose a lower one-sided 1 − α confidence interval [L(N 12 , T), 1] for θ P is available. It can be shown that [−1, U(N 12 , T)] is an upper one-sided 1 − α confidence interval for θ P if and [L(N 12 , T), U(N 12 , T)] is a two-sided 1 − 2α interval for θ P . Therefore, we focus on the construction of L(N 12 , T) only in this section. The R code will provide two (lower and upper) one-sided intervals and a two-sided interval, all are of level 1 − α. The first two are the smallest. The third is the intersection of the two smallest one-sided 1 − α/2 intervals. It may be conservative since the infimum coverage probability may be greater than 1 − α due to discreteness.

An inductive order on S P
Following Wang (2012), the construction of the smallest 1 − α interval [L(N 12 , T), 1] requires a predetermined order on the sample space S P . An order is equivalent to assigning a rank to each sample point, and this rank provides an order on the confidence limits L(n 12 , t)'s. Here we define that a sample point with a small rank has a large value of L(n 12 , t), i.e., a large point has a small rank. Let R(n 12 , t) denote the rank of (n 12 , t). Intuitively, there are three natural requirements for R: 2) R(n 12 , t) ≤ R(n 12 , t − 1), 3) R(n 12 , t) ≤ R(n 12 − 1, t + 1), as shown in the diagram below: Therefore, R(n − 1, 1) = 2 and a numerical determination is needed for the rest of R(n 12 , t)'s. Wang (2012) proposed an inductive method to determine all R(n 12 , t)'s, which is outlined below.
Step k: For k > 1, suppose the ranks, 1, . . . , k, have been assigned to a set of sample points, denoted by S k = ∪ k i=1 R i , where R i contains the ith largest point(s) with a rank of i. Thus, S k contains the largest through kth largest points in S P . The order construction is complete if S k 0 = S P for some positive integer k 0 , and R assumes values of 1, ..., k 0 .
Step k + 1: Now we determine R k+1 that contains the (k + 1)th largest point(s) in S P .
Part a) For each point (n 12 , t), let N (n 12 ,t) be the neighbor set of (n 12 , t) that contains the two points next to but smaller than (n 12 , t), see the diagram above. Let N k be the neighbor set of S k that contains all sets N (n 12 ,t) for (n 12 , t) in S k .
Part b) To simplify the construction on R, consider a subset of N k , called the candidate set from which R k+1 is going to be selected.
Then define R k+1 to be a subset of C k that contains point(s) with the largest value of L * P . We assign a rank of k + 1 to point(s) in R k+1 and let S k+1 be the union of R 1 up to R k+1 .
Since S P is a finite set and S k is strictly increasing in k, eventually, S k 0 = S P for some positive integer k 0 (≤ (n + 1)(n + 2)/2) and the order construction is complete.

The computation of the rank function R(N 12 , T) in R
There are three issues to compute the rank function R(n 12 , t): i) compute the infimum in f * (n 12 ,t) (θ P ); ii) determine the smallest solution of equation (4); iii) repeat this process on all points in S P .
These have to be done numerically.
Regarding i), use a two-step approach to search for the infimum when p T belongs to interval D(θ P ), i.e., in the first step, partition D(θ P ) into, for example, 30 subintervals, find the grid, say a, where the minimum is achieved; then in the second step, partition a neighborhood of a into, for example, 20 subintervals and search the minimal grid again. In total we compute about 50 (= 30 + 20) function values. On the other hand, if one uses the traditional one-step approach, one has to compute 600 (= 30 × 20) function values to obtain a similar result.
Regarding ii), the smallest solution is found by the bisection method with different initial upper search points and a fixed initial lower search point −1. The initial upper search point is the lower confidence limit of the previous larger point in the inductive search algorithm.
Regarding iii), use unique() to eliminate the repeated points in N k and use which() to search for R k+1 from C k (smaller) rather than N k .

The smallest one-sided interval under the inductive order
For any given order on a sample space the smallest one-sided 1 − α confidence interval for a parameter of interest can be constructed following the work by Buehler (1957); Chen (1993); Lloyd and Kabaila (2003) and Wang (2010). This interval construction is valid for any parametric model. In particular, for the rank function R(n 12 , t) just derived, the corresponding smallest one-sided 1 − α confidence interval, denoted by L P (n 12 , t), has a form L P (n 12 , t) = −1, if no solution for (6); the smallest solution of (6), otherwise, where f (n 12 ,t) (θ P ) = 1 − α and f (n 12 ,t) (θ P ) = 1 − sup p T ∈D(θ P ) ∑ {(n 12 ,t )∈S P :R(n 12 ,t )≤R(n 12 ,t)} p P (n 12 , t ; θ P , p T ), that are similar to (4) and (5). Two facts are worth mentioning. a) Among all one-sided 1 − α confidence intervals of form [L(N 12 , T), 1] that are nondecreasing regarding the order by the rank function R, L ≤ L P . So [L P , 1] is the best. b) Among all one-sided 1 − α confidence intervals of form [L(N 12 , T), 1], [L P , 1] is admissible by the set inclusion criterion (Wang, 2006). So [L P , 1] cannot be uniformly improved. These properties make [L P , 1] attractive for practice. The computation of L P is similar to that of the rank function R.

An inductive order and the corresponding smallest interval
Following Wang (2010), a rank function R I (X, Y) is to be introduced on S I . This function provides an order of the smallest one-sided interval L I (x, y). In particular, a point (x, y) with a small R I (x, y) is considered a large point and has a large value of L I (x, y). Similar to the rank function R in the previous section, R I should satisfy three rules: a) R I (n 1 , 0) = 1, b) R I (x, y) ≤ R I (x, y + 1), as shown in the diagram below: Repeating the process in the previous section, we can derive this new rank function R I on S I and the corresponding smallest one-sided 1 − α confidence interval [L I (X, Y), 1] for θ I by replacing (n 12 , t) by (x, y), D(θ P ) by D I (θ I ) and p P (N 12 , T; θ P , p T ) by p I (x, y; θ I , p 2 ). The only thing different is that for the case of n 1 = n 2 = n, R I generates ties. For example, R I (x, y) = R I (n − y, n − x) for any (x, y). However, the procedure developed is still valid for this case. Technical details were given in the Sections 2 and 3 of Wang (2010).

Examples Example 1: Exact intervals for the difference of two paired proportions θ P
We illustrate the usage of the PairedCI() function to calculate the exact smallest lower one-sided confidence interval [L P , 1] for θ P in (1) with the data from Karacan et al. (1976). In this study, 32 marijuana users are compared with 32 matched controls with respect to their sleeping difficulties, with n 11 = 16, n 12 = 9, n 21 = 3, and n 22 = 4. The second argument in the function is t = n 11 + n 22 = 20.
Function PairedCI() has the following arguments: PairedCI(n12, t, n21, conf.level = .95, CItype = "Lower", precision = . 1, grid.one = 3 , grid.two = 2 ) The arguments n12, t, and n21 are the observations from the experiment. The value of conf.level is the confidence coefficient of the interval, 1 − α, which is equal to the infimum coverage probability here. One may change the value of CItype to obtain either an upper one-sided or a two-sided interval. The precision of the confidence interval with a default value 0.00001 is rounded to 5 decimals. The values of grid.one and grid.two are the number of grid points in the two-step approach to search the infimum. The higher the values of grid.one and grid.two, the more accurate is the solution but the longer is also the computing time. Based on our extensive numerical study, we find that grid.one = 3 and grid.two = 2 are sufficient enough for the problem.
The use of marijuana helps sleeping because the interval [0.00613, 1] for θ P is positive.
The upper one-sided 95% interval and the two-sided 95% interval for θ P are given below for illustration purpose.

Example 2: Exact intervals for the difference of two independent proportions θ I
The second data set is from a two-arm randomized clinical trial for testing the effect of tobacco smoking on mice (Essenberg, 1952). In the treatment (smoking) group, the number of mice is n 1 = 23, and the number of mice which developed tumor is x = 21; in the control group, n 2 = 32 and y = 19. The function BinomCI() computes exact confidence intervals for θ I in (7), the difference of proportions between two groups.
In this clinical trial, the maximum likelihood estimate for the difference between two tumor rates θ I is calculated asθ # the lower one-sided 95% interval > lci <-lciall$ExactCI # extracting the lower one-sided 95% interval > lci [1] .133 1.
The lower one-sided 95% confidence interval for θ I is [0.133, 1]. Therefore, the tumor rate in the smoking group is higher than that of the control group.
The following code is for the upper one-sided and two-sided 95% confidence intervals.

Comparison of results with existing methods
Our smallest exact one-sided confidence interval [−1, U I ] for θ I is first compared to an existing asymptotic interval (Newcombe, 1998b) using the coverage probability. The coverage of an upper confidence interval [−1, U(X, Y)] as a function of θ I is defined as: Ideally, a 1 − α interval requires that Coverage(θ I ) is always greater than or equal to 1 − α for all the possible values of θ I .
The coverage for the exact upper 95% confidence interval [−1, U I ] and the asymptotic upper confidence interval based on the tenth method of Newcombe (1998b), which is the winner of his eleven discussed intervals, is shown in Figure 1. The two intervals are calculated by BinomCI() and the function ci.pd() in the package Epi. The left plot of Figure 1 shows the coverage against θ I ∈ [−1, 1] based on our exact method. As expected, it is always at least 95%. However, the coverage for the asymptotic interval may be much less than 95% as seen in the right plot of Figure 1. The coverage of the majority of θ I values is below 95% and the infimum is as low as 78.8% for a nominal level of 95%. The similar results are observed for the asymptotic confidence intervals based on other methods, including the one proposed by Agresti and Caffo (2000).
In light of the unsatisfied coverage for the asymptotic approaches, we next compare our exact intervals to the exact intervals by the PROC FREQ procedure in the software SAS. First revisit Example 2, where SAS provides a wider exact two-sided 95% interval [0.0503, 0.5530] for θ I using the EXACT RISKDIFF statement within PROC FREQ. This is the SAS default. The other exact 95% interval in SAS using METHOD = FMSCORE is [0.0627, 0.5292], which is narrower than the default but is wider than our two-sided interval. Also SAS does not compute exact intervals for θ P at all.
Two exact upper intervals produced by BinomCI() in the R package ExactCIdiff and the PROC FREQ procedure in SAS are shown in Figure 2. The smaller upper confidence interval is preferred due to the higher precision. Almost all the points in the figure are below the diagonal line, which confirms a better performance of the interval by BinomCI(). The average lengths of the two-sided interval for n 1 = n 2 = 10 and α = 0.1 are 0.636 and 0.712, respectively, for our method and the SAS default procedure. The newly developed exact confidence intervals have better performance than other asymptotic or exact intervals due to their guarantee on the coverage or because of their shorter length.

Summary
A group of three exact confidence intervals (lower one-sided, upper one-sided, and two-sided) are computed efficiently with the R package ExactCIdiff for each of the differences of two proportions: θ P and θ I . Each one-sided interval is admissible under the set inclusion criterion and is the smallest in a certain class of intervals that preserve the same order of the computed interval. Unlike asymptotic intervals, these intervals assure that the coverage probability is always not smaller than the nominal level.
], with a sample size of 320, takes less than one minute. Our exact interval is constructed by an inductive method. By nature, when there are many sample points, i.e., the sample size is large, deriving an order on all sample points is very time consuming. Thus the confidence limit on a sample point, which is located at the beginning (ending) part of the order, needs a short (long) time to calculate. Roughly speaking, when the sample size is more than 100, one would expect a long computation time for a two-sided interval. More details may be found at: http://www.wright.edu/~weizhen.wang/software/ExactTwoProp/examples.pdf.