Rank-Based Analysis of Unbalanced Repeated Measures Data

In this article, we have developed a rank (intra-subject) based analysis of clinical trials with unbalanced repeated measures data. We assume that the errors within each patient are exchangeable and continuous random variables. This rank-based inference is valid when the unbalanced data are missing either completely at random or by design. A drop in dispersion test is developed for general linear hypotheses. A numerical example is given to illustrate the procedure.


Introduction
In a typical clinical trial with repeated measures, a set of patients is chosen according to selection criteria defining a patient population.These patients are randomized into two or more treatment groups and are then observed at successive time points until the end of the trial.Even though trials might be designed with the intent to collect the same set of measurements from each patient, there are many practical situations where the investigators have to deal with observation (per patient) vectors of unequal length.Some patients might not have all measurements collected at all planned visits or time points.Patients could also be discontinued from the trial for one or more reasons.Over the course of the study, patients may miss visits, stop taking their assigned treatment, or be removed from treatment.Whatever the reasons for dropouts, drug related or not, withdrawals of patients from the trial result in unbalanced data.In this article, we have assumed that these unbalanced data are missing either completely at random or by design.In multivariate settings where missing values occur on more than one variable, the noncompleters are often a substantial portion of the entire data set.If so, ignoring them may be inefficient, causing a large amount of information to be discarded.Moreover, discarding them from the analysis will tend to introduce bias, to the extent that the incompletely observed cases differ systematically from the completely observed cases.The completely observed cases will be unrepresentative of the population for which the inference is usually intended: the population of all cases, rather than the population of cases with no missing data.Therefore, the results of this article have potential application in pharmaceutical, medical, and clinical research.
In this article, using a sum of Jaeckel (1972) type dispersion functions we develop a rank (intra-patient) based analysis of the time effects of the repeated measure design.Our analysis includes estimation of the time effects and their standard error of the estimates.Further, we develop a test based on the drop in dispersion for general linear hypotheses concerning the time effects, including an asymptotically distribution-free test of the interaction between the time effects and treatment groups.As discussed in Sections 4 and 5 our proposed analysis is robust.Furthermore, our analysis is completely invariant to the random effects due to the patients.and it allows each patient to have an observation vector of unequal length.Kloke, McKean, and Rashid (2009) developed a rank-based analysis of a general mixed model which can be used for this repeated measures design.In contrast to the analysis discussed in this paper, it is not invariant to the random effect due to the patients; although, its inference is adjusted for the random effect.Similar to the theory for the method discussed in this paper, the theory assumes that the number of patients become large.We compare this analysis to our proposed analysis in Section 3. Also, Rashid, McKean, and Kloke (2012) proposed an analysis for a multicenter clinical design, which is also based on a sum of Jaeckel dispersion functions.Their design was a mixed model with centers treated as random.Hence, for the situation of this paper, centers correspond to a patient.The asymptotic theory for this analysis assumes that the sample size of a center becomes large while the number of centers remains fixed.In contrast, the methodology developed in this paper assumes that the number of patients become large while the number of repeated measures remains finite.Hence, the theories for the two analyses are quite different.The analysis of Rashid et al. (2012) is briefly discussed in Section 3, also.
In Section 2 we present the repeated measures model discussed in this paper.In Section 3, we discuss the details of our proposed method and in Sections 4 and 5 present its asymptotic theory.In Section 6 we propose a consistent estimate of a necessary scale parameter for our procedure.An example which illustrates our methodology is presented in Section 7.

Repeated Measures Model
Consider an experiment with q treatment groups, j n subjects belong to the j th ( = 1, 2, , ) j q  treatment group and the k th subject in the j th group is supposed to be observed at p time points.A model for this kind of experiment is: = , where ijk Y is the observation corresponding to the ( = 1, 2, , ) ith i p  time treatment, the group treatment and the ( = 1, 2, , ) interaction between the i th time treatment and the j th group treatment, ijk e is the error term corresponding to the ( ) ijk th observation, and jk  is the random subject effect for the kth subject in the jth group.
We assume that ijk e are iid random variables with mean (or median) zero and scale parameter 2 , e  and jk  are iid random variables with mean (or median) 0 and scale parameter 2 .
  Let ( ) g t denote the common probability density function (pdf) of ijk e .We discuss tests of general linear hypotheses, but our main hypothesis of interest is that of interaction between the the group and times effects, i.e., 0 : 0 : = 0, , .
Model (2.1) can be rewritten as a repeated measures model as follows: = , is an identity matrix of order p and 1 p is a 1 p  vector of unity.
If the normality assumption is not reasonable, the experimenter would prefer to use a distribution-free procedure or an asymptotically distribution free procedure.For the theory of such procedures in this paper, we make the following assumptions, B1: The error vectors jk  ( = 1, 2, , ; = 1, 2, , j j q k n   ) for individual subjects are independent and continuous random vectors with the c.d.f. for incorporating missing values in the context of incomplete repeated measures data

Rank-Based Invariant Estimation
= ( , , , ) , treatment combination on the k th patient.To incorporate possible missing data, let * ( ) = ( ) Kloke, McKean, and Rashid (2009) extended the usual rank-based estimates of Jaeckel (1972) to general mixed models with covariates, including Model (2.1), for general scores functions.In this paper, we consider Wilcoxon scores only.For the Wilcoxon scores, using the missing data indicator ijk I , Jaeckel's objective function is given by where the subscript JR denotes the joint ranking of all of the residuals and n denotes the total sample size; i.e., convex function of  and the rank-based estimator of  is the minimizing value which we denote by ˆJR  .Kloke et al. (2009) developed the asymptotic theory of these estimates and their corresponding tests of general linear hypotheses for the general mixed model.Note that unlike a rank transform, the involved rankings are ranks of residuals not ranks of responses; see McKean and Vidmar (1992) for discussion.The estimates of the scale parameters involved in the associated JR inference are adjusted for the random effects, so the JR inference is asymptotically distribution-free.The JR estimator, though, is not invariant to the random effects.
For multicenter clinical trials, Rashid, McKean, and Kloke (2012) proposed an R estimator which is invariant to the random center effect.This procedure forms a Jaeckeldispersion function for each center by ranking the residuals within that center.Then the overall dispersion function is the sum of these center-specific dispersion functions.The resulting minimizer is called the MR estimator, where M denotes multiple.The asymptotic theory for the MR estimator is based on a fixed number of centers with the center sample sizes converging to infinity at relatively the same rate.In the situation of this paper, however, the number of repeated measures is fixed at p and for many practical cases p is small.Hence, asymptotic theory should cover the case where the number of subjects within each group gets large.For example, these are the assumptions underlying the theory for the JR estimator.Hence, the MR estimator is not appropriate in this case.As we show next, though, a similar type estimator is invariant to the random subject effects.

Rank-Based Invariant Estimates
Instead of doing a joint ranking of residuals, as in the JR estimator, we consider ranking the residuals within each subject.For the jk th subject, rank the residuals * that is, the dispersion functions jk D are invariant to the random effects for all combinations , j k .Further, for each patient this dispersion function is a nonnegative, continuous, and convex function of residuals * ijk W .An overall dispersion function based on all the patients in the q groups is defined as the sum of all ( ) ( = 1, 2, , ; = 1, 2, , )

Note that ( )
D  is also a non-negative, location-invariant, continuous, piece-wise linear and convex function of the ijk W , because the sum of several convex functions is again a convex function.Furthermore, ( ) Although, we have written ( ) D  as a function of  , note that just as the dispersion function is invariant to the random effect it is also invariant to the group effects j  .But it is not invariant to the time effects i  or the interaction effects ij  .Let  be the subvector of  which includes these time and interaction effects.Then our rank-based- Asymptotic theory for   is given in Section 4. This includes standard errors which can be used to determine asymptotically distribution free confidence intervals for the time and interaction effects.Consider, next, the hypotheses of interaction (2.2).Denote the minimizing value of the full model, ( is the upper  critical point of a 2  -distribution with ( 1)( 1) q p   degrees of freedom.This test is robust; see the discussion of Section 5.
For computation, the R algorithm discussed in Rashid et al. (2012) for the MR estimates can be used to obtain the RBIR estimates and the minimizing values of the dispersion function.

Special Case
In this subsection, we will examine how the R -estimator is related to the actual observations for a special case when = 1 q (only one grouping factor) and = 2 p (two levels of the repeated measures factor).( , ) = ( , ) = (0, ) = 12 is the intra-patient standardized ranks corresponding to is the intra-patient standardized ranks corresponding to
Thus, the R -estimate of  can be obtained by solving the following equation The above equation is equivalent to where I is an indicator function which takes a value 1 if In general (when ), no closed form R -estimate of the contrasts ( = 1, 2, , ) A competing R estimate, in this simple case, is the Hodges-Lehmann estimator of  , which is the median of the Walsh averages of the differences 2 theory for the Hodges-Lehmann estimate requires symmetry of the distribution of these differences, which follows from exchangeability of the errors ijk e .We illustrate this discussion with an example. be the median corresponding to the baseline and 2  be the median corresponding to the post-drug.Let 1 2


Our estimate of  is the median of the differences which is 7.The LS estimate is the mean of the differences which is 11.5.This disparity in the estimates is due to the large outlier, the postdrug PVC response for the 12 th patient.If this patient's values are omitted, the median of the differences remains the same, but the mean of the differences is 7.91.The LS estimator is highly influenced by outliers, while the rank-based estimate is robust.

Asymptotic Distribution Theory of the RBRI Estimators
As we outline in this section, with more details in the appendix, the asymptotic theory for the RBRI estimator follows from the asymptotic theory of the gradient of the dispersion function, much like the theory for rank-based estimators in general which is detailed, for example, in Chapters 3 and 5 of Hettmansperger and McKean (2011).So our discussion is brief.
As discussed at the end of Section 3, for the asymptotic theory, we assume that j n   , where j n is the number of patients in the j th group, = 1, , j q  .We assume that these j n s grow at similar rates; that is, for some postive constants j  > 0, = 1, , , Let 0  be the true value of  .Without loss of generality, we assume that 0 = 0  .

Theorem 4.1.2 Under exchangeablity of errors within each patient and assuming
as , N   where MVN stands for a multi-variate normal distribution, We assume that the rank of ( = 1, 2, , ) j j q   is 1. p  In other words, the levels of the timing treatment (repeated measures factor) within the j th grouping treatment are connected.Without this assumption, we are unable to perform the test for either the interaction between the grouping treatment and the timing treatment or the equality of the levels of the timing treatment.
Based on this, as discussed in Chapter 3 of Hettmansperger and McKean (2011), it follows that the influence function of  is bounded in response space.Hence, since the design matrix only contains dummy variables, the estimator   is bounded in both response and factor space.Therefore, the estimator is robust. By using arguments similar to those given in Rashid (1996), it can be shown that for large N ˆˆ

Drop in dispersion or reduction in dispersion tests
Theorem 5.1.3Under exchangeablity of errors within each patient and assuming as N   and where 2 r  stands for a chi-square distribution with r degrees of freedom.
The drop-in-dispersion test for general linear hypotheses depends on the R -estimates in the same way that the classical normal theory F test depends on the generalized least squares estimates.It is also analogous to 2log   in maximum likelihood procedure and has a similar interpretation.Note that the * D test is a mixed procedure involving both the ranks and the residuals.As discussed on page 197 of Hettmansperger and McKean (2011), however, the influence function of the drop in dispersion test statistic is bounded in response space and for this repeated measures design also in factor space too.Hence, the drop in dispersion test statistic is robust.On the other hand, the maximum likelihood procedure is not robust.It is quite sensitive to outliers as the example in Section 7 illustrates.

Estimation of 
In this section, we briefly sketch another consistent estimator of the scale parameter  .
For our discussion, we can assume without loss of generality that the true regression coefficients ik  are 0.
Denote the residuals for the k th subjects of the j th group by = , = 1, , .
Note that we can write the differences as these residuals as ˆ= ( ) Each bracketed term is invariant to the random effects, so the differences of these residuals are invariant to the random effects.Furthermore, the pdf of difference of the random errors in the first bracket is symmetric about 0.
Let M denote the number of such residuals and, for convenience, renumber them as 1 to M .Next form the Walsh averages of these residuals, i.e., Denote the the number of Walsh averages by   Then as a final estimator take the average of these estimators.Since each of the estimators is consistent and there are only a finite number of them, this average is consistent also.This was the estimator used in the example of Section 7. In practice, as long as the missing is within practical bounds, this estimator is consistent.

Example
One of us consulted on a data set that would serve as an example but due to confidentiality matters it cannot be used.Instead, we have simulated a data set which is similar to real data set in most aspects.The design for the real data set was a repeated measures design with two groups (Placebo and Treated) and the primary outcome (a lipid) measured at 4 time periods on a total of 81 patients (22 in the Placebo arm and 59 in the Treated arm).There were many outliers in the real data set.The interaction between the groups and times was the effect of interest.
For our simulated data set, we used the same design, including sample sizes and the number of repeated measures.We set the cell medians (means) at the values in We generated data with a high intraclass correlation coefficient.The random errors were simulated from the Laplace distribution with variance set at 1.30, while the random effects are generated normal variates with variance set at 16 also.So the intraclass correlation coefficient is 0.924.In the cell with Group at level 2 and Time at level 1 two points were contaminated.
Using the algorithm discussed in Kloke et al. (2009), we computed the JR-fit for the vector of parameters  .Figure 7.1 show the q q  and Interaction Plots based on the JR Fit.The q q  clearly shows the two contaminated data points.Actually, except for these two points, the rest of the data are good.The profile plots clearly indicate interaction between the group and time.The top profile is that of the treated group, while the bottom profile is based on the fit for placebo group.The interaction plot detected the drop in response between times 2 and 3 for the placebo group.Table 7.3 displays these time effects estimates and their standard errors.For comparison, we also included the JR estimates and their standard errors.The estimates and standard errors for the two methods are similar.They both detect the patterns in the table of the true contrasts.The RBRI seems to be a little more precise.Both robust tests (RBRI and JR) found interaction to be highly significant.They easily detected the pattern in the table of population medians.On the other hand, the mle analysis was impaired by the two outliers.Practical interpretations would differ between the robust and traditional (mle) analyses.In most applications, with this p -value, the traditional analysis would accept no interaction and concentrate on the main effects.On the other hand, the robust analyses would turn to the estimation of contrasts similar to those in Table 7.3.

Conclusion
When the response variable is continuous, the assumption of multivariate normal is not always reasonable.In addition, in other situations involving a continuous variable, the actual distribution may be unknown.Thus, the use of standard parametric procedures is patient to criticism regarding both optimality and validity.
In this article, we have developed rank based inferences for unbalanced repeated measures models with incomplete observations.These inferences are robust alternatives to parametric analyses which are invariant to the random effect due to patients.Note that nonparametric methods may also be of interest as a means of confirming the results of a parametric analysis or in assessing the sensitivity of the results to the assumptions of the selected parametric model.

Proof of Theorem 5.1
Suppose   minimizes ( ) Q  and H   minimizes ( ) Q  under the restriction = 0. A By using arguments as those given in Rashid(1996), it can be shown that the asymptotic behavior of ˆ( ) ( ) H D D    can be determined by that of ( ) ( ) It can be shown that,


Then the dispersion function for the k th subject is given by


We obtain the R estimate of  by minimizing the above dispersion function with respect to . Note the partial derivative of ( ) D  is given by * 2 =1

D
 exist almost everywhere and are given by the coefficients of .
are the rank-based analogues of the traditional log likelihood tests; see McKean and Hettmansperger (1976) or Chapter 3 of Hettmansperger and McKean (2011) for a review.In this subsection, we outline these tests based on the RBRI estimators for general linear hypotheses of the form

1 
of ordered Walsh averages.McKean and Hettmansperger (1976) proposed estimating  by standardized length of a confidence interval for the intercept given by * .The generally recommended value of  is 0.05.For a fixed effects linear model with symmetric errors and residuals based on a Wilcoxon fit, ˆMH  is a consistent estimator of  .In our situation, the errors are symmetrically distributed.but due to the random effects there are dependencies of the random errors within a subject.One way of avoiding this is to first consider the case of no missing data.Then consider the sets of residuals, one for each pairing; that is, for each pair , = 1, , ; <

Figure 7 . 1 :
Figure 7.1: q q  and Interaction Plots of the JR Fit.The top and bottom profiles are the profiles for the treated and placebo groups, respectively.Using the MR algorithm, as discussed in the last paragraph of Section 3.1, we obtained the RBRI fit of the vector of parameters  .Based on this fit, the McKean- Hettmansperger estimate of  is 1.479, as discussed in the last paragraph of Section 6.


be the true value of .j Without loss of generality, we assume that 0

m
For fixed j and k , under exchangeablity of errors within each patient, ijk R are uniformly distributed over the integers from 1 through .jkBy using arguments similar to those given inRashid (2000), it can be shown that

Following
Rashid (1996), it can be shown that ( ) Q  where   is a generalized inverse of  (defined in (Error!Reference source not found.)satisfying = .4.12) follows from (8.25).

Table 3 .1: Before and After EKG Readings of Patients
16)where  stands for asymptotically equivalent.It can be seen that the reduced model RBRI estimate, ˆ, H  is asymptotically equivalent to a linear function of the full model

Table 7 .2: True Cell Medians for Times and Their Contrasts
For each of these sets obtain the McKean-Hettmansperger consistent estimator of  .

Table 7 .
2, which are close to the estimated cell medians of the real data set.This table also shows the true values of the contrasts

Table 7 .4: Estimates of Contrasts and Standard Errors for the time Effects for the RBRI and JR Procedures.
For the test of the hypotheses of interaction (2.2), Table7.4displays the test statistics and their p -values in the next table.Results are given for the RBRI and JR-analyses and also the maximum likelihood analysis (assuming normal error distributions).The mle computation was performed by the R function lme.