Testing Equality in Ordinal Data with Repeated Measurements: A Model-Free Approach

1 Introduction

In clinical trials or health-related studies, we often encounter the patient response on an ordinal scale, for example, worse, same or better. Arbitrarily assigning scores (such as –1, 0, 1) to these ordinal categories for arithmetic operation can be inappropriate due to the fact that the relative distance between any two successive ordinal categories is not really equal or even comparable. One may also have difficulty in interpreting the mean of these arbitrary scores in terms of practical meaning. Grouping multiple categories into a single category to reduce ordinal outcomes into dichotomous responses can cause the loss of efficiency.

To reduce the number of patients in clinical trials, taking more than one measurement on each patient frequently arises in practice. Research on repeated ordinal responses has been intensive [1, 2, 3–5, 6, 7, 8, 9]. Francom et al. [2] applied a family of structural log-linear models and required investigators to assign the scores to represent the relative distance between ordinal categories. Agresti [1] addressed use of generalized linear models with different link functions and linear predictors. Ware et al. [9] as well as Kenward and Jones [4] discussed various approaches to analyze repeated categorical measurements. Parsons et al. [7] considered the proportional odds logistic regression with a range of working correlation models. All these publications focused discussions on model-based methods. By contrast, the methods proposed here is model-free and does not assume any parametric form for the data structure. Furthermore, there is no need to assume or specify any particular dependence structure between repeated measurements. Also, rather than concentrating attentions on testing equality between treatments in ordinal data with repeated measurements, there were publications [10, 11] discussing and deriving procedures for testing positive (quadrant) dependence under various parameter constraints and marginal modeling in two-way tables with ordinal categories. Agresti and Coull [12] focused discussion on model-based approach as well and addressed testing hypothesis against order-restricted alternatives. The purposes of these papers are different from what we focus here is to derive model-free procedures for testing equality of treatments, while the dependence between measurements taken within patients is nuisance effect.

Using the generalized odds ratio (GOR) [13], we develop in this paper model-free procedures for testing equality of treatments and derive interval estimators for the relative treatment effect in ordinal data with repeated measurements. We further develop a simple procedure for testing the treatment-by-period interaction. We discuss the usefulness and limitations of test procedures developed here. To illustrate the use of these procedures, we consider two real-life data sets, one taken from a trial studying the gender effect on pain scores on an ordinal scale after hip joint resurfacing surgeries [7] and the other taken from a double-blind randomized trial comparing an active hypnotic drug with a placebo in insomnia patients with respect to the ordinal category of time to falling asleep [2].

2 Notation and Methods

Consider comparing two treatments in a randomized clinical trial, in which we randomly assign ng patients to group g (=1, 2). For a randomly selected patient i(=1, 2,…, ng) from group g, we let Yizg denote the patient response at period z (=1, 2), where Yizg takes one of L ordinal labels Cj (j=1, 2,…, L) with C1<C2<C3<⋯<CL. Let nrsg denote the number of patients with (Yi1g=Cr,Yi2g=Cs) among ng patients in group g. The random vector (n11g,n12g,⋯n1Lg,n21g,n22g,...,n2Lg,⋯,nL1g,nL2g,...,nLLg)′ then follows the multinomial distribution with parameters ng and (π11g,π12g,⋯π1Lg,π21g,π22g,...,π2Lg,⋯,πL1g,πL2g,...,πLLg)′, where πrsg denotes the cell probability that a randomly selected patient from group g has the vector of responses (Yi1g=Cr,Yi2g=Cs). We let “+” denote the summation over that particular subscript. For example, πr+g(=∑sπrsg) denotes the probability P(Yi1g=Cr) that a randomly selected patient from group g (=1, 2) has the response Yi1g=Crat period 1. Following Agresti [13], we define ΠC1=∑r=1L−1∑r′=r+1Lπr+1πr′+2, denoting for period 1 the probability that the response of a randomly selected patient from group 2 is larger than that from group 1. We further define ΠD1=∑r=2L∑r′=1r−1πr+1πr′+2 denoting for period 1 the probability that the response of a randomly selected patient from group 1 is larger than that from group 2. The GOR of responses between groups 2 and 1 at period 1 is simply defined as G1=ΠC1/ΠD1. When there is no difference in treatment effects at period 1 between the two groups, G1=1. If the treatment at period 1 in group 2 tends to increase the response as compared with that in group 1, G1>1. If the former tends to decrease the response as compared with the latter, G1<1. Similarly, we define the GOR of responses between groups 2 and 1 at period 2 as G2=ΠC2/ΠD2, where ΠC2=∑s=1L−1∑s′=s+1Lπ+s1π+s′2 and ΠD2=∑s=2L∑s′=1s−1π+s1π+s′2. Note that when L=2, Gz(z=1, 2) reduces to the regular odds ratio (OR) of responses in dichotomous data.

Note that we can estimate πrsg by the unbiased consistent sample proportion estimator πˆrsg=nrsg/ng, and thereby, we can estimate G1 at period z=1 between groups 2 and 1 by

where ΠˆC1=∑r=1L−1∑r′=r+1Lπˆr+1πˆr′+2 and ΠˆD1=∑r=2L∑r′=1r−1πˆr+1πˆr′+2. Using the delta method [14], we can show that an estimated asymptotic variance of Gˆ1(1) with the logarithmic transformation is given by [13, 15]

Note that we define ∑r=L+1Lπˆr+1=∑r′=L+1Lπˆr′+2=0. Similarly, we define ∑r=10πˆr+1=∑r′=10πˆr′+2=0.

Following the same arguments as for deriving Gˆ1 (1) and Var^(log(G^1))(2), we obtain a consistent estimator for G2 between groups 2 and 1 at period 2 as

where ΠˆC2=∑s=1L−1∑s′=s+1Lπˆ+s1πˆ+s′2 and ΠˆD1=∑s=2L∑s′=1s−1πˆ+s1πˆ+s′2. Again, using the

delta method, we can show that an estimated asymptotic variance of Gˆ2(3) with the logarithmic transformation is given by

Also, we define ∑s=L+1Lπˆ+s1=∑s′=L+1Lπˆ+s′2=0 and ∑s=10πˆ+s1=∑s′=10πˆ+s′2=0 in (4).

Note that Gˆ1(1) and Gˆ2(3) are actually correlated. For clarity, we include the tedious details in derivation of Cov^(log(G^1),log(G^2))in Appendix I.

3 Examples

To illustrate the use of point estimators Gˆ1 and Gˆ2, test procedures and interval estimators developed here, we first consider the trial studying hip joint resurfacing surgeries [7]. To alleviate pain and debilitation caused by osteoarthritis, rheumatoid arthritis, fractures and other hip related problems, hip replacement surgery is a widely-used procedure. To assess the failure rate and prognosis after hip joint resurfacing, we may use Harris scores [16], of which one important component is the pain scores coded on an ordinal scale: none (1), slight (2), mild (3) and moderate or marked (4) pain in hip joint. We summarize in Table 1 the data regarding the frequency distribution according to pain scores taken at two and five years after hip joint resurfacing for 58 patients by gender published elsewhere [7]. Because there were very few patients with moderate or marked pain (4), we grouped these patients and patients with mild pain (3) into one category without loss of much information. It is of interest to study whether there is a difference in pain scores between genders over time after surgeries. In terms of our notation, we define “females” and “males” as group (g =) 1 and 2, as well as “two years” and “five years” after the surgery as period (z =) 1 and 2. From Table 1, the numbers of patients are n1=21 and n2=37 for the two comparison groups. We can calculate, for example, the estimated cell proportion πˆ111=n111/n1=7/21. The marginal percentage for each row (or column) is simply equal to the corresponding marginal total divided by the total number of patients ng(g=1, 2). For example, the marginal percentage πˆ1+1 in the first row for females is 0.476 (=10/21). Using marginal percentages for the rows between the two sub-tables (Table 1), we can then calculate the GOR estimate Gˆ1 (1) of pain scores at two years post surgeries as

[0.476×(0.135+0.108)+0.333×0.108]/[0.333×0.757+0.190×(0.757+0.135)]=0.360. Similarly, using marginal percentages for columns between the two sub-tables, we obtain Gˆ2=0.637 at five years post surgeries. As an example, we can interpret Gˆ1=0.360 as that the odds of a randomly selected male with the pain score higher than a randomly selected female at two years post surgeries is 0.360, given no ties in pain scores between genders. Furthermore, when employing interval estimators (9), we obtain 95 % confidence intervals [0.133, 0.974] and [0.250, 1.623] for G1 and G2, respectively. Since the upper limit of the 95 % confidence interval for G1 falls below 1, there is significant evidence at the 5 % level that males tend to have pain scores lower than females at two years after surgeries. However, since Gˆ2(=0.637) is closer to 1 than Gˆ1(=0.360), this difference in pain scores between genders seems to decrease as the time increases from two years to five years after surgeries. Note that the above resulting 95 % confidence interval for G2 covers 1. Thus, the difference in pain scores between genders is no longer significant at the 5 % level after five years post surgeries, despite that males is still likely to have pain scores lower than females (because Gˆ2 < 1). When using test procedures (5) and (6) for testing the overall equality of pain scores across two periods between genders, we obtain p-values as 0.086 and 0.131. On the basis of the above results, we may conclude that although males tend to fall, as compared with females, in categories with low pain scores especially at two years post surgeries, the difference in pain scores between genders seems to become smaller in the long-term results.

When assuming a normal random effects proportional odds model [17], we may employ Proc Glimmix in SAS [18] to study the difference in pain scores between genders on the basis of the model-based approach. Using the data in Table 1, we have obtained the parameter estimate (and its estimated standard error (SE)) for the relative gender effect of males versus females to be 0.961 (SE=0.6148). This leads the p-value for testing the equality of pain scores across the two periods between genders is 0.124, which is similar to those obtained by use of test procedures (5) and (6). Also, note that the parameter estimate 0.961 is larger than 0. Thus, males tend to fall in categories with lower pain scores than females. This inference is identical to that obtained on the basis of the GOR focused here.

We may sometimes encounter a trial in which the patient response taken at period (z=) 1 actually represents the baseline response. In this case, we may apply the test procedure (7), developed for testing group-by-period interaction, to study whether there is a relative treatment effect. To illustrate this point, we consider the data (Table 2) taken from a double-blind randomized clinical trial comparing an active hypnotic drug (g=1) with a placebo (g=2) in patients with insomnia [1, 14, 2]. The outcome of interest is to respond the question “How quickly did you fall asleep after going to bed?” and is recorded on a four-point ordinal scale (< 20, 20–30, 30–60, > 60 in minutes). Each participated subject was asked this question twice, one after a one-week placebo washout period for both groups and the other at the conclusion of a two-week treatment period. Using the data in Table 2, we obtain Gˆ1=1.031 and Gˆ2=1.883. Because subjects received placebo at the first one-week washout period in both groups, the estimate for the GOR of responses is expected to be, as shown Gˆ1=1.031, around 1 due to randomization. At the end of two-week treatment period, the estimate for the GOR of responses becomes Gˆ2=1.883 with 95 % confidence interval using interval (9) given by [1.270, 2.792]. Because this lower confidence limit falls above 1, one may conclude that the hypnotic drug changes the GOR over periods and tends to reduce time to falling asleep. When using test procedures (6) and (7), we obtain p-value as 0.002 and 0.004; these small p-values also suggest that the hypnotic drug can significantly reduce the time to falling asleep. We note that all these test results are essentially similar to those based on the proportional odds model with the cumulative logits found elsewhere [1]. When using the summary test procedure (5) using data in Table 2, we obtain the p-value 0.055. Note that since the GOR at the initial period is around 1, there is non-negligible probability that the time to falling asleep for patients in the placebo group is smaller than those in the other group at the initial washout period. As noted previously, use of a summary test procedure (5) over two periods in this case may not only be senseless (due to different treatments received at two periods), but also lack power.

4 Discussion

We can employ, as demonstrated here, the GOR to measure the relative treatment effect without the need to assume any specific parametric model. Since the GOR has a simple interpretation and is easily understood, the GOR is of use in ordinal data. In fact, the GOR is closely related to the gamma correlation [19], a commonly-used measure of the strength of association between two ordinal variables. We refer readers to some publications on estimation of the GOR and its applications under other situations [20, 21, 22–24]. When repeated measurements are taken at the same time (or there are no period effects), one may employ the Dirichlet-multinomial distribution to model the intraclass correlation between repeated measurements [25]. As considered in the above two examples, however, repeated measurements on patients are often taken at different time intervals in clinical trials. The period effect is likely to exist and is required to be incorporated to avoid bias in data analysis [26, 27]. The methods based on the Dirichlet-multinomial model for cluster sampling without accounting for the period effect would not be appropriate for use in situations focused here.

The proportional odds model is probably one of the most commonly-used models to analyze ordinal data. Just like all model-based approaches, we can easily extend the proportional odds model to account for confounders (if there were) or accommodate other general situations. However, the implicit assumption of the proportional odds model can be badly violated by many bivariate distributions [13, 28, 29]. Furthermore, when applying Proc Glimmix in SAS [18] based on the random effects proportional odds model to ordinal data with repeated measurements, we need to assume that the random effects (accounting for the intraclass correlation between repeated measurements) due to patients follow a normal distribution. This normal assumption for random effects can be difficult to be justified. By contrast, the proposed method is model-free. It does not require the random effects due to patients to follow the normal distribution, not does assume any parametric models for the data structure. Thus, our methods are applicable despite of various parametric models for the underlying data structure and distribution assumptions for the patient random effects. Furthermore, the point estimators, test procedures and interval estimators developed here can all be expressed in closed forms. Readers may employ these test procedures and estimators by use of a pocket calculator even without knowledge of any statistical software. The interpretation of the GOR is, as illustrated in examples, easily understood. When there are confounders in a trial of a large size, we may extend the methods proposed here by use of stratified analysis with strata determined by the combined levels of confounders [23]. But we want to note that the model-based approach can be preferable to the model-free approach proposed here if there are many covariates to adjust for a trial of a small or moderate size.

Finally, we note that using similar arguments as above, it is straightforward to extend the results to accommodate the cases with three or more periods. We outlines the extension of results presented here to accommodate three periods in Appendix II.

In summary, we have developed model-free test procedures for testing equality of treatments in ordinal data with repeated measurements. We have further derived interval estimators for the relative treatment effect measured by the GOR. We recommend use of the summary test procedure to improve power when treatments received at two periods in a group are the same and there is no treatment-by-period interaction. However, we should cautiously employ this summary test procedure when there is a treatment-by-period interaction. The bivariate test procedure can be of use in this case. We further outline the extension of results to accommodate three periods. The results, findings and discussions should have use for biostatisticians and clinicians when they encounter ordinal responses with repeated measurements.