A Comparison of Exact Tests for Trend with Binary Endpoints Using Bartholomew’s Statistic

1 Introduction

Suppose realizations of a binary variable are recorded in a k-group setting where the groups are defined by ordered categories. The data associated with such experiments can be organized as a 2×k table; see Table 1. Most commonly, xi, the number of responses in group i, is assumed to follow a binomial distribution with parameters πi and ni, i=1,2,…,k. Furthermore, we restrict interest to scenarios where xi is assumed independent of xj for all i≠j. The parameter πi is often referred to as the response probability in situations where the quantity represents the likelihood of positive response when the subject is exposed to the stimuli or intervention associated with the group. Under these situations, it is often of interest to test the equality of the response probabilities, and the null hypothesis takes the form

Consider an experiment with four dose groups, each fixed by the design to have 10 subjects. There is expected to be an increasing probability of response with increasing dose level, and doses administered are 0, 2.5, 25, and 250 units. In such scenarios, use of tests of the null in favor of the general alternative may result in an inefficient procedure in that they ignore the ordinal information which defines the groups. Rather, a test for the ordered alternative (a.k.a. “simple order hypothesis” [1]) is appropriate, and the alternative may be written as

Although the standard test due to Cochran and Armitage [2, 3] may be applied, the inherent weakness of this approach is the dependence on a choice of scores, which can be somewhat arbitrary and present the researcher with a dilemma in that the operating characteristics are dependent on the choice. Given that the procedure is based upon an assumed linear relationship between the chosen score and the probability of response, an incorrect choice in scores can result in losses in efficiency and vastly different conclusions for a given data set. In the aforementioned example, say responders totaled 0, 1, 4, and 3 in the four groups receiving the lowest to highest doses, respectively; see Table 2. Using dosage amount as the score, the Cochran–Armitage procedure does not detect any trend in response probabilities, with p-value 0.2615. However, if the transformation log{dose+0.01} is taken as a score, the resulting p-value is 0.0332.

Bartholomew [4] proposed a test utilizing isotonic regression that addressed trend alternatives via a likelihood-ratio approach, the test statistic being a variant of the standard Pearson chi-square statistic. In that the null distribution depends on an unknown nuisance parameter in finite samples, Bartholomew published tables of asymptotic critical values for his statistic in the balanced sample size case [5] for use in practice. Poon [6] conducted a Monte Carlo study indicating that the trend test based on these critical values sometimes violated the nominal type I error rate; Collings et al. [7] corroborated these simulation results. Leuraud and Benichou [8, 9] also used simulation to assess the statistical properties of Bartholomew’s test, finding the procedure to have good power properties under a number of alternatives.

Given the aforementioned type I error control problems, it is natural to alternatively consider use of the exact probability distribution. Complications arise in that, conditional on n1,n2,…,nk, the distribution functions of x1,x2,…,xk are dependent upon the unknown nuisance parameter π1=π2=⋯=πk=π. A number of classic solutions exist for the special case of k=2. For testing equality of two binomial probabilities, Fisher’s famous approach [10] conditions on the marginal totals of an observed 2×2 table. Eliminating dependence on the nuisance parameter π by restricting attention to 2×2 tables with the same marginal totals as the observed table, the p-value for Fisher’s exact test is the sum of hypergeometric probabilities. Several exact tests for trend have been developed, many listed by Agresti and Coull [11], all of which used the exact conditional approach. The term conditional in this context refers to the fact that the corresponding null distributions are those conditioning on the total number of successes, and these procedures can be associated with a high degree of discreteness producing conservative tests.

For the 2×2 case, Barnard [12] proposed a more powerful test of H0:π1=π2=π that called for conditioning only on n1 and n2. An expression containing πx1+x2(1−π)n−x1−x2 is then maximized over the binomial parameter π to obtain the p-value. Berger and Boos [13] augment this procedure by maximizing only over a 100(1−β)% confidence interval for π, then penalizing the resulting p-value by adding β. In the present context, we consider the more general problem of testing equality of k≥2 binomial parameters against a trend (simple order) alternative.

An approximate unconditional approach to eliminating the nuisance parameter is the notion of simply replacing it with an estimate, for instance one obtained via maximum likelihood [14]. This procedure provides an estimated (E) p-value. Lloyd [15] proposed using this E p-value as an ordering criterion and then maximizing the null likelihood, resulting in the E + M (Estimation + Maximization) p-value. Shan et al. [16] investigated this approach in the context of detecting monotone trend using the Cochran–Armitage test statistic and subsequently using the Baumgartner-Weiss-Schindler statistic [17].

Due to the conservative nature of conditional tests, and the violations of type I error that accompany tests based on asymptotic critical values as observed by Poon [6], it is desirable to develop exact unconditional tests for trend for use with qualitatively ordered samples using the Bartholomew statistic.

Section 2 gives background information on Bartholomew’s test for trend. Section 3 details the proposed exact testing procedures. Section 4 reports type I error and power properties of the examined procedures. Section 5 presents an example of the tests’ application, and Section 6 contains concluding remarks.

2 Background

Bartholomew’s [4] test of equality against the simple order alternative as applied to binomial counts is based on the statistic

where x=(x1,x2,…,xk), πˉ is the maximum likelihood estimator (MLE) of the nuisance parameter (πˉ=∑i=1kxi/∑i=1kni), and the πˆi∗s are the isotonic regression estimates of the πi s. This specific application of Bartholomew’s test, as described in Barlow et al. [1], begins with determination of the vector πˆ∗=(πˆ1∗,πˆ2∗,…,πˆk∗), the MLE of π=(π1,π2,…,πk) subject to the order restriction [18],

If the usual MLEs πˆi=xi/ni are observed to be isotonic (i.e. πˆ1≤πˆ2≤⋯≤πˆk), then πˆ∗=πˆ, for πˆ=(πˆ1,πˆ2,…,πˆk). If the πˆis are not in the order corresponding to Ha, then the order-restricted MLE πˆ∗ may be obtained using the pool-adjacent-violators algorithm (PAVA) [19]. For more on the estimates computed by PAVA for independent binomial samples and the MLE of π under the simple order restriction, see Silvapulle and Sen [20]. Recall the hypothetical example presented earlier, with ni=10,i=1,2,3,4, and πˆ=(0,0.1,0.4,0.3). These πˆis do not concur with the order restriction shared by the parameter space, since 0.3/≥0.4. Following the steps dictated by PAVA [4], the violating quantities, namely 0.4 and 0.3, are pooled together, giving (4+3)/(10+10)=0.35. Replacing both 0.4 and 0.3 by the pooled value 0.35 gives an array of nondecreasing estimates, and hence πˆ∗=(0,0.1,0.35,0.35).

Although the statistic TB is essentially Pearson’s chi-square statistic applied to the isotonic regression estimates πˆ∗[1, 7, 11], the distribution of TB is not as convenient as that of Pearson’s chi-square. In the context of binary data, the sample proportions πˆi would have only an asymptotic normal distribution, and the isotonic modification is an additional complicating factor. The asymptotic null distribution of TB may be expressed as ∑a=1kpaχa−12 (referred to as the chi-square-bar distribution [4]) where χa−12 is a chi-squared random variable with a−1 degrees of freedom (χ02≡0), and pa is the probability that the vector πˆ∗ has a distinct ordered values. While Barlow et al. [1, pp. 134–149] present formulas and recursions to calculate pa, the expression is a function of the number of groups k, and

even for moderate k this may require several numerical integrations.

[21]

Determining values for pa can be highly prohibitive, necessitating the use of simulation to obtain approximate values [22]. Critical values are presented by Bartholomew [5] for the equal sample size case after a precise determination of the required pa values.

Suppose the number of observed responders in the k groups of a particular data set is x0=(x01,x02,…,x0k), and we define s=∑i=1kx0i, n=∑i=1kni. By conditioning on s, which is a sufficient statistic, we obtain the exact conditional p-value as,

where Ωcond(x0)={x:TB(x)≥TB(x0) and ∑x=s}.

Poon [6] generated tables of binomial counts from the conditional distribution and approximated the type I error and power of the exact test using the Bartholomew statistic under a variety of parameter settings, each with 500 simulations. Conditioning on the number of responses provides a p-value expression that is free of the nuisance parameter π. The exact type I error rate and power for this approach are computed according to the prescription of Mehta et al., who evaluated exact unconditional power of the conditional Cochran–Armitage test [23].

3 Methods of unconditional exact testing

Three unconditional exact testing procedures are now described in detail and distinguished from the conditional version.

4 Test properties

The exact methods detailed in Section 3 and the exact conditional approach from Section 2 are compared with the test using the asymptotic critical values presented in Bartholomew [5]. For additional comparison, the Cochran–Armitage test with equally spaced scores based on the asymptotic and exact conditional null distributions are included. Type I and type II error rates were calculated through complete enumeration, i.e. presented rates are exact and not obtained through Monte Carlo simulation. In the case of the conditional procedure, the unconditional error rates are presented for comparability which are the rates across possible values of s. Complications arise in computing TB when the observed π¯ is equal to 0 or 1. In these instances, since such tables show no evidence of increasing trend, we set TB≡0.

Figure 1 displays exact type I error plotted over the range of the nuisance parameter for the four exact tests and the test based on asymptotic critical values. For scenarios defined by the number of groups (3 or 4) and the sample size within each group (10 or 20), nuisance parameter values between 0 and 1 by an increment of 0.01 were used in implementation of the exact tests. An alternate approach is a two-stage procedure involving an initial pass to identify possible regions for the global maximum, followed by a more refined search of the identified regions. The critical values of TB used to evaluate the asymptotic test were 3.820 and 4.528, respectively, for k=3 and k=4 [5]. The profiles corresponding to the exact conditional tests using Bartholomew’s statistic and the Cochran–Armitage statistic are labeled, respectively, as C(TB) and C(TCA); likewise, the tests based on asymptotic critical values of these same statistics are labeled A(TB) and A(TCA). Profiles based on the unconditional maximization approach, confidence interval approach, and the application of Lloyd’s approach are labeled M, CI, and E + M, respectively.

Figure 1

Exact type I error rates of the Bartholomew test based on the unconditional maximization (M), confidence interval (CI), expectation and maximization (E + M), conditional (C(TB)), and asymptotic-based (A(TB)) approaches, and the Cochran–Armitage test based on conditional (C(TCA)) and asymptotic-based (A(TCA)) approaches

As can be seen in Figure 1, the unconditional tests preserve the nominal α-level in all four scenarios. In the first panel of Figure 1, the M, E + M, and A(TB) profiles coincide. While adequate in this situation (k=3, ni=10), the test using the asymptotic critical value of TB does not control type I error in the remaining panels, as the A(TB) profile is seen to exceed 0.05. The Cochran–Armitage test based on an asymptotic critical value violates the nominal level of type I error in all scenarios. Type I error of the unconditional tests is uniformly less conservative than that of the conditional tests. For extreme values of the nuisance parameter, the unconditional tests are seen to have significant advantage over the conditional test.

Power is first summarized in Figure 2 by choosing a parameter setting similar to those listed in Neuhauser [31], of the form πa=[0.1, 0.1+γ2, 0.1+γ] for the case of k=3, and πa=[0.1, 0.1+γ3, 0.1+2γ3, 0.1+γ] for the case of k=4, with γ∈(0,09). These settings are linear in nature, describing the “center” of the alternative space. It is difficult to distinguish the power curves in the first panel of Figure 2, except to note that in general, power is lowest for C(TCA) and A(TB), both tests utilizing the Cochran–Armitage statistic. This holds for the remaining three panels of Figure 2. With some effort, power for C(TB) can be identified as having slightly worse power than A(TB), M, CI, and E + M, which are all nearly identical. The third panel of Figure 2 indicates that for k=4 and ni=10, A(TB) and CI enjoy a minor power advantage over the other tests, the top performer being dependent on the value of γ. We remind the reader of the type I error violations associated with the test based on TB in conjunction with asymptotic critical values, thereby detracting from the power advantages.

To investigate power further, in Figure 3 we consider a parameter setting closer to the edge of the alternative space, choosing πa=[0.1, 0.1, 0.1, 0.1+γ] when k=4, and πa=[0.1, 0.1, 0.1, 0.1+γ] when k=4, again with γ∈(0,09). The power properties are largely unchanged from the earlier parameter setting, with the tests based on TCA suffering a power disadvantage, and power curves for the remaining tests grouping tightly together. The CI test shows slightly superior power properties in the third panel of Figure 3. Power for C(TB) can be seen below the main cluster of curves in all four panels.

Figure 2

Exact power for center of alternative space of the Bartholomew test based on the unconditional maximization (M), confidence interval (CI), expectation and maximization (E + M), conditional (C(TB)), and asymptotic-based (A(TB)) approaches, and the Cochran–Armitage test based on conditional (C(TCA)) and asymptotic-based (A(TCA)) approaches

Figure 3

Exact power for edge of alternative space of the Bartholomew test based on the unconditional maximization (M), confidence interval (CI), expectation and maximization (E + M), conditional (C(TB)), and asymptotic-based (A(TB)) approaches, and the Cochran–Armitage test based on conditional (C(TCA)) and asymptotic-based (A(TCA)) approaches

5 Example

We consider the hypothetical example given in Fleiss et al. [32], tabulating one-month release rates as a function of initial severity measured on an ordinal scale. The observed counts appear in Table 3, and testing results are found in Table 4. For the unconditional exact approaches, all 2×4 tables with the identical column totals are enumerated, and the Bartholomew statistic corresponding to each table is computed. p-Values for the M and E + M tests are calculated by maximizing the expressions in Section 3.3 over π∈(0,1) by increments of 0.01. The p-values for the conditional tests are calculated by enumerating all 2×4 tables with the same row and column totals as the original data.

Figure 4 depicts the p-value profiles of the M and E + M tests. As the height of these curves depends on π, the M and E + M p-values are the respective maxima of these profiles. The limits of the 99.9% Wald confidence interval for π obtained from the observed data are (0.493, 0.807), represented in the figure by vertical lines. The CI p-value is obtained by adding β=0.001 to the maximum of the M test profile between these confidence limits. The horizontal dashed lines in the figure represent p-values for the conditional exact tests based on statistics TB and TCA. The values of these statistics obtained from the example data are 28.81 and 23.91, respectively. Significance can be assigned to these values using standard asymptotic methods by appealing to the limiting chi-square and chi-square-bar distributions. As the group sample sizes in this example are not balanced, the p-value from the asymptotic chi-square-bar distribution can be identified as <0.005 using the Table A.2 of Barlow et al. [1]. All tests detect a highly significant decreasing trend in the probability of release as initial severity increases at a 0.05 nominal significance level.

Figure 4

p-Value profiles for example in Fleiss et al.

6 Conclusion and recommendation

The test using the asymptotic critical values of Bartholomew does not have acceptable type I error control. The exact conditional test is guaranteed to control type I error, but is conservative in all instances, sometimes severely so. The maximization, confidence interval, and Estimation + Maximization tests all maintain type I error control and show markedly less conservatism than the conditional test. Furthermore, each unconditional test exhibits a modest gain in power over conditional tests.

Based on the type I error and power properties, no single unconditional test is uniformly preferable to the other two. The E + M test requires far more computational effort than the M and CI tests. The CI test involves a penalty term that can affect the p-value by orders of magnitude, in addition to the minor task of requiring calculation of a confidence interval for the nuisance parameter. By way of simplicity and computational convenience, we recommend Maximization as the unconditional method of choice.