Testing Equality of Treatments under an Incomplete Block Crossover Design with Ordinal Responses

1 Introduction

When studying non-curable chronic diseases, including angina pectoris, epilepsy, hypertension, asthma, etc., we may often consider using a crossover design to reduce the number of patients needed for a parallel group design [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]. When there are more than two treatments under comparison, the trial duration for a crossover design can be much longer than that for a parallel group design if each patient is to receive every treatment in use of the former. The longer the duration of a trial, the higher is the patient risk of being lost to follow-up. Furthermore, a lengthy trial duration can cause the difficulty in recruiting patients into a trial and ensuring patients to closely follow a study protocol. To alleviate these concerns, we may consider using an incomplete block crossover design, in which each patient is to receive only a subset of treatments [1]. For example, consider the double-blind placebo controlled crossover trial comparing 12μg and 24μg of formoterol solution aerosol with a placebo [1]. For practical reasons, it was decided that each patient could receive only two of the three treatments: the placebo, 12μg and 24μg of formoterol solution. Although there were some publications [1, 12, 15, 16, 17, 18] on the incomplete block crossover design, all these focused discussion on either continuous or binary data. The discussion on testing equality of treatments with ordinal responses under an incomplete block crossover design is limited [1, 7, 8].

Because ordinal responses are not on an interval scale, it is generally not appropriate to apply arithmetic operation to ordinal data [19]. In practice, we may commonly assign arbitrary scores to ordinal responses and do hypothesis testing with use of the t-test. Since the relative distances between consecutive categories in ordinal data are not truly comparable, converting ordinal responses into a universally agreeable score scale is difficult. Also, how to provide a meaningful and easily-understood summary measure based on these arbitrarily assigned scores to quantify the relative treatment effect can be challenging. On the other hand, if we dichotomize the ordinal responses into binary outcomes, the test procedures for binary data will probably lose efficiency.

In this paper, we propose use of the generalized odds ratio (GOR) for paired samples [20, 21, 22] to measure the relative treatment effect on patient responses in ordinal data. We focus our discussion on an incomplete block two-period crossover trial comparing three treatments with ordinal responses. We derive asymptotic test procedures based on the weighted-least-squares (WLS) and Mantel-Haenszel (MH) estimators [23] for testing equality between treatments. We further derive the exact test procedures for testing equality of treatments for small-sample cases. We employ Monte Carlo simulation to evaluate the finite-sample performance of these test procedures in a variety of situations. We use the data taken as a part of the crossover trial [24] comparing the low and high doses of an analgesic with a placebo for the relief of pain in primary dysmenorrhea over the first two periods to illustrate the use of these test procedures.

2 Notation, assumption and methods

Consider comparing two experimental treatments A and B with a placebo P under an incomplete block crossover design with two periods. We let X-Y denote the treatment-receipt sequence of receiving treatment X at period 1 and then crossover to receive treatment Y at period 2. Suppose that we randomly assign ng patients to group g(=1,2,3,4,5,6), where g=1 denotes the group with P-A treatment-receipt sequence; g=2 denotes the group with A-P treatment-receipt sequence; g=3 denotes the group with P-B treatment-receipt sequence; g=4 denotes the group with B-P treatment-receipt sequence; g=5 denotes the group with A-B treatment-receipt sequence; and g=6 denotes the group with B-A treatment-receipt sequence. For patient i(=1,2,…,ng) assigned to group g(=1,2,3,4,5,6), we let Yiz(g) denote the ordinal outcome of the patient at period z(=1,2), and take one of possible ordinal values Cj, where C1<C2<C3<⋯<CL. We let Xiz1(g) denote the indicator function of treatment-receipt for treatment A, and Xiz1(g)= 1 for patient i assigned to group g at period z receiving treatment A, and = 0, otherwise. Similarly, we let Xiz2(g) denote the indicator function of treatment-receipt for treatment B, and Xiz2(g)= 1 for the corresponding patient at period z receiving treatment B, and = 0, otherwise. We let 1iz(g) represent the indicator function of period, setting 1iz(g)= 1 for period z = 2, and = 0, otherwise. We further let μi(g)denote the random effect due to the ith subject in group g, and assume μi(g)’s to independently follow an unspecified probability density fg(μ). We assume that one can apply an adequate wash-out period on the basis of our subjective knowledge to nullify the carry-over effect. As noted elsewhere [1, 4, 11, 12, 25, 26, 27], if we cannot ensure this assumption to hold, we may not wish to employ the crossover design. For patient i(=1,2,…,ng) in group g(=1,2,3,4,5,6), we assume that the joint conditional probability of (Yi1(g),Yi2(g)) between periods 1 and 2, given the random effect μi(g)fixed, satisfies

whereηAPand ηBPdenote the respective effect of treatments A and B relative to placebo P, as well as γ denotes the effect of period 2 versus period 1. Based on model (1), the GOR of responses [20, 21, 22] on a given patient i in group g when he/she has covariates (Xi21(g),Xi22(g),1i2(g)) at period 2 versus when he/she has covariates (Xi11(g),Xi12(g),1i1(g)) at period 1 is, by definition, equal to

When ηAP = 0, we can see from eq. (2) that the GOR of responses remains unchanged despite of receiving treatment A or placebo P. When ηAP > 0, taking treatment A tends to increase the patient response as compared with taking placebo P, given all the other covariates fixed. When ηAP < 0, taking treatment A tends to decrease the patient response as compared with taking placebo P. Similar interpretations of ηAP are applied to parameters ηBP and γ. We define the GOR of responses for treatment A versus placebo P and that for treatment B versus placebo P as GORAP=exp(ηAP) and GORBP=exp(ηBP), respectively. Also, we define the GOR of responses for treatment B versus treatment A as GORBA=exp(ηBP−ηAP).

On the basis of model (1), for a randomly selected patient i from group g the probability that the patient response Yi1(g)at period 1 is less than his/her response Yi2(g)at period 2 is

Similarly, for a randomly selected patient i from group g the probability that the patient response Yi1(g) at period 1 is larger than his/her response Yi2(g)at period 2 is

Note that because we do not assume any parametric p.d.f. for fg(μ) in the following discussion, our approach is semi-parametric.

For simplicity in notation, we define ΠC(g)=P(Yi1(g)<Yi2(g)) and ΠD(g)=P(Yi1(g)>Yi2(g)). From eqs (3) and (4), we can see that for a randomly selected patient i from group g the GOR of patient responses between periods 2 and 1 is

We denote for a randomly selected patient i from group g the probability P(Yi1(g)=Cr,Yi2(g)=Cs) byπrs(g), where r=1,2,…,L, and s=1,2,…,L. Thus, we have

These represent the probability that a randomly selected patient i from group g has the response at period 2 higher than his/her response at period 1, and the probability that a randomly selected patient has the response at period 1 higher than his/her response at period 2, respectively. When L = 2, GOR(g)reduces toπ12(g)/π21(g), the OR of responses between periods 2 and 1 in binary data with matched-pairs. On the basis of model (5), we can express the GOR of responses for treatment A versus placebo as (Appendix I)

Let nrs(g)denote the number of patients in group g (= 1, 2, 3, 4, 5, 6) with the vector of patient responses (Yi1(g)=Cr,Yi2(g)=Cs)among ng patients. The random cell frequencies {nrs(g)|r = 1, 2, 3,⋯, L, s=1, 2, 3,⋯, L} then follow the multinomial distribution with parameters ng and {πrs(g)|r = 1, 2, 3,⋯, L, s=1, 2, 3,⋯, L}. Note that we can estimate πrs(g)by the unbiased consistent estimator πˆrs(g)=nrs(g)/ng. We define nC(g)=∑r=1L−1∑s=r+1Lnrs(g) and nD(g)=∑r=2L∑s=1r−1nrs(g). When substituting πˆrs(g) for πrs(g) in ΠC(g)and ΠD(g), we obtain ΠˆC(g)=nC(g)/ng and ΠˆD(g)=nD(g)/ng. These lead us to obtain the estimator GORˆ(g)= ΠˆC(g)/ΠˆD(g). Using the delta method [22, 28], we obtain the estimated asymptotic variance Varˆ(log(GORˆ(g)))=(ΠˆC(g)+ΠˆD(g))/(ngΠˆC(g)ΠˆD(g)). When substituting GORˆ(g)for GOR(g) in eq. (7), we obtain the following three consistent estimators for GORAP(=exp(ηAP)) as

For convenience, we define three 2 × 2 tables consisting of cell frequencies (f11k,f12k,f21k,f22k)(for k = 1, 2, 3) corresponding to eq. (8) as

and

When testing H0:GORAP=1 versus Ha:GORAP≠1, we first consider use of the WLS summary test procedure [23] based on eq. (9). We will reject H0:GORAP=1at the α-level if

where

for k = 2, 3, LGORˆAP(1)=log((f111f221)/(f121f211))1/2, LGORˆAP(k)=log((f11kf22k)/(f12kf21k))for k = 2, 3, and Zα is the upper 100(α)th percentile of the standard normal distribution. Note that if fijk= 0 for some observed frequencies in a 2 × 2 table k, we cannot employ the test procedure (10). We may apply the commonly-used ad hoc arbitrary adjustment for sparse data by adding 0.50 to each observed frequency fijk in this particular table k.

When the observed marginal frequencies nC(g) and nD(g) are not large, the WLS test procedure may lose accuracy because the weights Wkin eq. (10) can be subject to a large variation. This may lead us to consider use of the MH summary test procedure [28, 29]. When comparing treatment A with placebo, we will reject the null hypothesis H0:GORAP=1 at the α-level if the test statistic

When both the numbers of patients nC(g) and nD(g) are small, the asymptotic WLS and MH test procedures may lose accuracy. Thus, we may consider use of the following exact test procedure. Define ndis(g)=nC(g)+nD(g). Given ndis(g) fixed, we can show that nC(g) follows the binomial distribution with parameters ndis(g) and ΠC(g)/(ΠC(g)+ΠD(g))(= GOR(g)/(1+GOR(g))). UnderGORAP=1, the conditional probability distribution of f11k, given f+1k,f+2k,f1+kand f2+k fixed, is given by the hypergeometric distribution:

where ak≤f11k≤bk, ak=max0,f1+k−f+2kand bk=minf+1k,f1+k for k = 1, 2, 3.

Thus, the joint conditional probability distribution off111,f112andf113 is simply [30]

wheref_11=(f111,f112,f113)′, f_u+=(fu+1,fu+2,fu+3)′for u=1, 2, and f_+v=(f+v1,f+v2,f+v3)′ for v=1, 2. Given an observed value f_11o= (f111o,f112o,f113o)′, if the following p-value, calculated as [23, 30]

where C = f_11P(f_11f_+1,f_+2,f_1+,f_2+)≤Pf_11o|f_+1,f_+2,f_1+,f_2+, is less than a small given α-level, we will reject H0:GORAP=1. Note that the exact test (14) is actually a direct extension of Fisher’s exact test to a series of 2 × 2 tables [23, 28, 30] with modifications to accommodate ordinal responses.

As shown in Appendix I, we can easily modify asymptotic test procedures (10) and (11) and the exact test procedure (14) to account for testingH0:GORBP=1 (or H0:GORBA=1) with replacing fijk by fijk∗ (or fijk∗∗), where fijk∗ and fijk∗∗ are defined in (26) and (28) (Appendix I), respectively.

3 Monte Carlo simulation

To evaluate and compare the performance of the WLS, MH and exact procedures for testing equality between treatments, we employ Monte Carlo simulation. By use of the conditional arguments, we do not need to estimate the nuisance period effect γ in use of these test procedures. We arbitrarily set γ equal to 0.10 in the simulation. Furthermore, our approach is valid for any assumed distribution for the random effects μi(g). We consider the cases in which random effects μi(g)are independent and identically distributed (i.i.d.) as the normal distribution with mean 0 and standard deviation σ=0.5,1; as well as μi(g)are i.i.d from a gamma distribution with shape parameter α=½ and scale parameter β=1,2. We cover the situations in which the relative effect of treatment A versus the placebo, ηAP=0.0,0.50; the relative effect of treatment B versus treatment A, ηBA=0.0,1.0 (such that the relative effect of treatment B versus the placebo, ηBP=ηAP+ηBA); and the number of patients n (=n1=n2=…=n6) per group n=10,15,25. Note that when n=10, these include the cases in which the expected number of patients with discordant responses ndis(g)(=nC(g)+nD(g)) is as small as approximate 3.71–4.84 patients. Note also that all test procedures proposed here depend on only the marginal totals nC(g) and nD(g)instead of individual cell frequency nrs(g). Thus, there is no need to consider the number of ordinal levels L, and so is nrs(g). For each configuration determined by a combination of the above parameter values, we write programs in SAS [31] and generate 10,000 simulated samples of n patients per group, each having the bivariate responses (Yi1(g), Yi2(g))^’ with probability P(Yi1(g)<Yi2(g)|μi(g)) and P(Yi1(g)>Yi2(g)|μi(g)) given by model (1), to calculate the simulated Type I error and power at the 0.05-level for a given test procedure. Recall that the power function of a test procedure with rejection region is, by definition, the probability that the sample points fall into the rejection region. The power function will give Type I error when H0 is true, and will give power when H0 is false. Therefore, the simulated Type I error for a given test procedure can be calculated as the proportion of 10,000 simulated samples for which we reject H0 when H0 is true. Similarly, the simulated power for a given test procedure can be calculated as the proportion of 10,000 simulated samples for which we reject H0 when H0 is false. For readers’ information, the SAS program for our simulation can be accessible at http://edoras.sdsu.edu/~kjl/exactso.htm.

4 Results

We summarize in Table 1 the estimated Type I error (in boldface) and power of using the WLS, MH and exact procedures for testing H0:ηAP=0 and testingH0:ηBP=0 at the 0.05-level when μi(g)are i.i.d. as the normal distribution with mean 0 and standard deviation σ=0.5,1. For example, when ηAP=0 and ηBP=1.0, the entries corresponding to procedures for testing H0:ηAP=0 are Type I errors, while those corresponding to procedures for testing H0:ηBP=0 are powers (Table 1). We can see that both the MH and exact tests can perform well, while the WLS test can be conservative especially when n is small (say, 10). We note that the MH test can be consistently of more power than the WLS and exact tests in almost all the situations considered in Table 1. For example, when σ=0.50, ηAP=0.50, ηBP=1.50 and n = 15, the powers for testing H0:ηAP=0 and H0:ηBP=0 are 0.175 and 0.825 for the MH test, while these powers are 0.124 and 0.781 for the WLS test, and are 0.132 and 0.731 for the exact test. We also note that the power for all test procedures increases as the number of subjects n increases, but decreases as the variation σ of responses between patients increases. Since all the findings on the performance of the WLS, MH and exact procedures with respect to Type I error, as well as the relative order of powers between these test procedures hold when μi(g)are i.i.d. from the gamma distribution, we do not present these results for brevity. These results are, however, to readers upon request.

5 An example

Consider the data (Table 2) taken as a part of a crossover trial comparing an analgesic at low (L) and high (H) doses with a placebo (P) for the relief of pain in primary dysmenorrhea patients over the first two-periods [24]. Here, we refer the low and high doses as treatments A and B. There were 86 patients randomly assigned to the six groups: P-L (g = 1); L-P (g = 2); P-H (g = 3); H-P (g = 4); L-H (g = 5); and H-L (g = 6). At the end of each treatment period, each patient was assessed the extent of relief on the ordinal scale: none (coded as 1), moderate (coded as 2) and complete (coded as 3). When applying the WLS and MH procedures, as well as the exact procedure to test H0:GORAP=1, we obtain p-values 0.0119, 0.0013 and 0.0002. Similarly, when applying these procedures with replacing fijk by fijk∗ to test H0:GORBP=1, we obtain all p-values as 0.0017, 0.0000 and 0.0000. Thus, there is evidence that taking either low or high dose of the analgesic can help the relief of pain in primary dysmenorrhea as compared with the placebo. Furthermore, when applying the corresponding test procedures with replacing fijk by fijk∗∗ to test H0:GORBA=1, we obtain p-values 0.7994, 0.9986, and 1.0000. There is no evidence that taking the high dose can help, as compared with the low dose, the relief of pain among these patients with primary dysmenorrhea.

Note that if we assume the proportional odds model using cumulative logit [28, 31] with normal random effects [19] due to patients, we can apply Proc GLIMMIX [31] to analyze the data in Table 2. When employing this SAS procedure (http://edoras.sdsu.edu/~kjl/exaordex1.htm), we obtain the parameter estimates (and their estimated standard error (SD)) of the relative effect for the low dose versus the placebo, the high dose versus the placebo, and the high dose versus the low dose are: –1.7656 (SD = 0.3898), –2.3404 (SD = 0.4011), and –0.5748 (SD = 0.3537). On the basis of these estimates (of which the signs are all < 0), we may conclude that both the low dose and high dose can significantly improve, as compared with the placebo, the relief of pain at the 5 % level; both p-values are < 0.0001. Furthermore, though the high dose can improve the outcome of patients as compared with the low dose, this improvement is not significant at the 5 %-level. All the above results are essentially similar to those reported previously for the complete block crossover design over a three-period trial [32].

6 Discussion

We do not recommend, as noted previously [1, 4, 11, 12, 25, 26, 27], use of the crossover design if we cannot ensure ourselves to nullify the carry-over effects with an adequate wash-out period. On the other hand, if there are carry-over effects due to earlier treatments, we note that the test procedures proposed here can still be valid for use under the simple carry-over model (Appendix II). Also, we note that although one may apply the estimator as given in Appendix II for the difference in carry-over effects to test whether there are differential carry-over effects, we do not recommend using this test to determine whether the assumption of no carry-over effects holds. This is because the concerns raised by Freeman [33] and Senn [34] for using the two-stage test procedure suggested by Grizzle [2].

When employing the test procedures developed here, we do not need to assume any parametric distribution for the random effects due to patients. Thus, our procedures is, as noted before, semi-parametric. Also, the number of patients for a crossover trial is often small. The exact test procedure (14) can be of use in practice. By contrast, one needs to assume the random effects due to patients independently follow the normal distribution in use of Proc GLIMMIX [35]. Also, note that the proportional odds model can be badly violated by many bivariate distributions [20, 36, 37]. Furthermore, note that Wald’s test can be invalid for use if the number of patients in a trial is small and the data are sparse.

We note that when the number of subjects per group n is small (say, 10), use of the WLS procedure can be conservative, while the MH and exact test procedure can perform well (Table 1). We further note that the MH test procedure is generally of more power than the other two procedures in the situations considered here. Because use of the MH procedure does not involve any sophisticated numerical procedure and can be calculated by a hand calculator, we may recommend the MH test procedure for general use when n is not large. We may use of the exact test procedure if one has the concern of normal approximation for a small n. When the number of subjects n per group is large (say, 40), however, we want to note that the WLS procedure can be of more power than the MH and exact procedures on the basis of Monte Carlo simulation (not shown here).

If we wish to study the relative period effect, we may apply similar ideas as above to derive the corresponding procedures for testing H0:γ=0. For example, we can easily see from eqs (16)–(21) (in Appendix I) that the GOR of responses for period 2 versus period 1 is

Following similar arguments as for comparing the treatment effect, we can derive on the basis of (15) the WLS, MH and exact procedures for testingH0:GOR21=1. Finally, note that when the patient response is dichotomous, the GOR reduces to the common OR for paired sample data. Thus, the MH and exact test procedures include those for testing equality of treatments in binary data [18] under the random effects logistic regression model as special cases.

In summary, we have derived the WLS, MH and exact procedures for testing equality between treatments under an incomplete block crossover design with ordinal responses. We have evaluated and compare their performance in a variety of situations based on Monte Carlo simulation. We have noted that the proposed test procedures are valid for use in the presence of simple carry-over effects. We have compared the proposed test procedures with use of Wald’s test procedures assuming the normal random effects proportional odds model. We have noted that the proposed test procedures include those for testing equality of treatments in binary data as special cases. The results, findings and discussions should have use for biostatisticians and clinicians when they employ a two-period crossover design to compare three treatments in ordinal data.