Analysis of repeated measures incomplete block designs using R-estimators

STATISTICS& PROBABILITY LETTERS ELSEVIER

Statistics & Probability Letters 22 (1995) 213-221

Analysis of repeated measures incomplete block designs using R-estimators M.

Mushfiqur Rashid

Mathematical Sciences Department, Worcester Polytechnic Institute, 100 Institute Road, Worcester, MA 01609, USA

Received June 1993; revised January 1994

Abstract Using a dispersion function, a rank-based inference is developed for repeated measures incomplete block designs (IBD). Asymptotic distributions of rank estimators follow from the approximate linearity of the gradient vector of the dispersion function. Three asymptotically equivalent test statistics are developed for testing the equality of the treatment effects. Estimation and testing are tied together. Keywords: Dispersion function; Approximate linearity; Quadratic approximation; Robust estimates and tests

1. Introduction The most commonly used incomplete block designs (IBD) are proper, binary and equireplicate. Using a dispersion function we develop a robust inference for repeated measures IBD in this paper. We will consider repeated measures IBD which are proper, binary and equireplicate. Suppose there are n subjects and t treatments. Each subject receives k ( < t) treatments and every treatment appears in r(~< n) subjects. It is assumed that the designs are connected. Let Yo be the observation corresponding to the ith treatment and the jth subject. A model for repeated measures incomplete block designs (IBD) is given by Yo=nlj(#+~i+e~j),

i = 1,2 . . . . . t ; j =

1,2 ..... n,

(1)

where nij = 1 if the ith treatment appears in the jth subject and nii = 0 otherwise, /~ = [overall mean], ~ti = [effect of the ith treatment], El= ~~ = 0, ei~ is the random error term. The random error vectors ~j's for different subjects are i.i.d. In parametric inference we assume that ~j,-~M V N I [0, tr2Ak ×k(P)], where tr 2 = var(eo) and Ak ×k(P) is an equicorrelation matrix. The parameter p measures subject effect; if p = 0, there is none. See Berry (1987); and Samuels et al. (1991) for further discussions. F o r inferences concerning repeated measures IBD, one approach is to find a sensible estimator ~ of p and apply the results of the generalized linear model assuming that p =/3. For further discussion see Arnold (1981). Another approach is discussed in papers concerning recovery of inter-block information (Yates, 1940; Rao, 1947a). For excellent development of recovery of inter-block information, see John (1987), and John (1971). 0167-7152/95/$9.50 © 1995 Elsevier Science B.V. All rights reserved SSDI 0167-7 1 5 2 ( 9 4 ) 0 0 0 6 9 - K

214

M.M. Rashid / Statistics & Probability Letters 22 (1995) 213-221

However, there are m a n y instances where the multi-normality assumption of the error vectors may not be valid and in such situations one would prefer to use a distribution-free procedure or an asymptotically distribution-free procedure. Since the purpose of this paper is to develop a rank-based theory for repeated measures IBD, which will be a robust alternative to the parametric analyses described earlier, we assume that (a) the error vectors, e~'s are independent and each error vector has a continuous c.d.f. F ( x l , X 2 . . . . . Xk) = F ( x ' l , X ' 2 . . . . . X t k ) where (x'l, x'2 . . . . . X t k ) is a permutation of(x1, X2 . . . . . Xk) , (b) the bivariate density f ( . , . ) for any two components of ej is continuous in R 2, and (c) ~ ~ f ( x , x) dx < ~ . We develop asymptotic results assuming n ~ oo. In order to develop a robust inference for repeated measures IBD, we will define a rank-based dispersion function for each subject using an idea similar to that of Jaeckel (1972). This dispersion function is a piece-wise linear function of the residuals with coefficients determined by the corresponding ranks where ranking is done separately for each subject. Then we will define an overall dispersion function for the repeated measures I B D based on all the subjects. Like the sum of squared-errors (SSE) of the linear model theory, this dispersion function will play a key role for rank based inference for the repeated measures IBD. In fact, in Section 2, we will obtain rank estimators of the treatment effects and their asymptotic distributions as n --* oo. In Section 3, we develop three asymptotically equivalent test statistics for testing the equality of the treatment effects. The first test statistic is based on the drop in dispersion. The second test statistic is based on the gradient of the dispersion function evaluated at the null hypothesis. This test statistic reduces to Durbin's (1951) statistic. The third test statistic is based on a quadratic in R-estimators. We will see that the asymptotic distributions of R-estimators of ~'s and test statistics will depend on a scale parameter z = 1 / [ x / ~ S ~ _ ~ f ( x , x) dx]. We need a consistent estimator of r to make our tests asymptotically distribution free. We develop a consistent estimator of z in Section 3.2. The parameter z will play a role similar to crz(1 - p) in normal theory repeated measures designs. If the components of each error vector are independent, z becomes 1/[x/~S_~o f 2 ( x ) d x ] and tr2(1 - p) becomes a 2. Therefore, in this paper subject effect is incorporated through z, and if z = 1 / [ x / ~ _ ~ f 2 ( x ) d x ] there is none. We will not consider the case when ej,,,Fj(x~, x2 . . . . . Xk), j = 1, 2 . . . . . n (i.e. when the error vectors for different subjects are not identically distributed) because we will not be able to find a consistent estimator of zj = 1 / [ ~ / ~ fi(x, x)dx], j = 1, 2 . . . . . n, where f j(., . ) i s the bivariate p.d.f, corresponding to any two components of cj, since we are assuming that k is fixed and n ---, ~ . However, in this case one can use Durbin's (1951) statistic for testing equality of the treatment effects.

2. Rank estimation and asymptotic distribution theory In this section we develop R-estimators of ~i's and their asymptotic distributions. The asymptotic distributions of R-estimators will be required to develop test statistics for testing the equality of the treatment effects. The dispersion for thejth subject based on Wilcoxon scores using an idea similar to that of Jaeckel (1972) is

oj(a) =

y. , , J R ( r - " - i= 1 'J [_ k + 1

[ Yi~ - ctl],

(2)

where R ( Y u . - - ~t~) is the intra-subject rank of the residual corresponding to the ith treatment and the jth subject. The combined dispersion function is D(a) = Z~= ~Dj(a). It can be shown that this dispersion function is nonnegative, continuous, location-free and convex in a. Hence, a rank estimate ~i of a can be obtained by minimizing the dispersion function. One can use Newton's method or a nondifferentiable optimization algorithm such as Nelder and Mead (1965) to obtain rank estimate d of a. The domain of D(a) is divided into a finite number of convex polygonal subsets, on each of which D(a) is a linear function of a. The partial

M.M. Rashid / Statistics & Probability Letters 22 (1995) 213-221

215

derivatives exist almost everywhere and the negatives of the partial derivatives are given by the vector S(a) with the ith element

~/-;-g~(R(YiJ--Cti)

st(a) = x / l z 2., nij< j=l (

~}

.... k+l

i = 1, 2, . . . .

, t.

(3)

We need S(a) to construct a quadratic a p p r o x i m a t i o n to D(a) as well as in determining the a s y m p t o t i c distributions of R-estimators of ~[s. Let a ° be the true value of a. In the following subsections we will develop the limiting distribution for x / ~ ( d - a ° ) . The strategy for determining the limiting distribution is as follows: (1) from a linear a p p r o x i m a tion to S(a) we will construct a quadratic a p p r o x i m a t i o n to D(a); (2) for the value ~ which minimizes the quadratic a p p r o x i m a t i o n , it will be shown that x/~(~i - a °) has a limiting multi-normal distribution; and then (3) we will show that x/~(¢i - ~i) ~ 0 as n ~ oo. 2.1. Asymptotic distribution o f S(a °) under the true value a ° In this subsection we will develop the a s y m p t o t i c distribution of the negative of the gradient vector S(a °) under the true value a °. W i t h o u t loss of generality, let us assume that a ° = 0. U n d e r the true value a ° = 0 and the a s s u m p t i o n (a), the rank vectors for different subjects will be independently and uniformly distributed over the p e r m u t a t i o n of the integers 1 through k. Therefore, it is easy to show that Eo[S(0)] = 0,

Varo[si(0)]

-

r(k - 1). k+l

'

and

Covo[s/(0),sr(0)] -

- 2ii,

k+l'

(4)

where r is the replications of each treatment and 2ii, = Y.3= 1 nijni,g is the n u m b e r of times both the ith treatment and the i'th t r e a t m e n t a p p e a r in the same subject (i # i' = 1, 2 . . . . . t ) . Let m = l i m . ~ r/n, and Z = [tr,,] where tr,, = (1/(k + 1)) lim.~o~ - 2ii,/r(i # i') and aii = (k - 1)/(k + 1). Further, the rank of Z is t - 1 since we are considering connected designs. Lemma 1. Suppose the assumption (1/x/~)S(0) ~ M V N ( O , m r ) as n ~ oo.

(a) holds for

the

model (1), and

true

value a ° = 0 .

Then

L e m m a 1 is a special case of a t h e o r e m in Skillings and M a c k (1981, p. 172). Therefore its p r o o f is omitted. 2.2. Linear approximation to S(a) and quadratic approximation to D(a) A linear a p p r o x i m a t i o n to the gradient vector is i m p o r t a n t to the development of the a s y m p t o t i c distribution theory of r a n k estimators and tests for repeated measures IBD. First we would like to find the expected value of the gradient vector under the true value a ° = O. The result is given in the following lemma. Lemma

2. Suppose the assumptions (a)--(c) hold for the model (1), and the true value a ° = O. Then

eoES(a)-I =

r -

-~ T

a + o(11 a If),

where r is the replications of each treatment and o( [Ia II)/ll a II ~ 0 as IIa II ~ 0. Proof. We have, using (3) j=l

'J~

k + 1

i =

1, 2 .....

t.

(5)

M.M. Rashid / Statistics & Probability Letters 22 (1995) 213-221

216

But

Eo[R(YIj - 0tl)] = 1 +

~,

ni'jPo(Yi,j -- o~i, < Yij -- O~i)

i'~:i=l

= 1+

~

ni'jg(o4, ~i').

i'¢:i=1

(6)

E x p a n d i n g e ( ~ , 0ti,) as a function of cq and oq, a r o u n d (0,0) and retaining the first two terms, we have t

g(ai, ~i') =

o0

f(Yi'j, Y i i ) d y r j dyij + (cti, - ~i)

f ( Y i j , Yij)dyij + 0 [or2 + ~2] 1/2 --o0

1 = ~ + (~i' -- cq)

f ( y , y)dy + 0 [ct2 + ~2] t/2, --ct3

w h e r e f ( x , y) is a bivariate density function corresponding to the c.d.f.F. Hence, using (5) and (6), we have =

f(y, y)dy + o(llal[)

~ n ~ n v j ( ~ . - ~i)

1 1 (k + 1) z [ - '~" -- )-i2 . . . . .

--

r(k - 1). . . . . - 4 d o + o(llall).

(7)

Therefore, we have

Eo [s(a)] where z =

:

-

r z "c

a

+

o(llall),

1/[x/~I~_oof(y, y)dy].

El

N o w we state and prove some useful theorems. T h e o r e m 1 gives a key result for linearization of the gradient vector. T h e o r e m 1.

Suppose the assumptions (a)-(c) hold for the model (1) and the true value of A is d ° = 0. For every

e >0,

limPo[ x/~{!S(-~A)-~S(O,}+mszt n-. oo LII

>~e]=0.

(8,

z

Proof. Using L e m m a 2, we can write for large n,

E°x//-n[~S(~A)-~S(O']

- - m'~A ' r'

(9,

Also, for i = 1, 2 . . . . . t,

Varo

12r of

2

M.M. Rashid / Statistics & Probability Letters 22 (1995) 213-221

217

assuming without loss of generality the ith treatment appears on the first subject and r/n ~ m as n ~ 0o. But, for i = 1, 2 ..... t,

{R(Yil--~)--R(Yil))

2~k2 ,

and lim R

Yil

--

Ai

--

R(Y,)

nel,(I rt,_.7~
= lim

n~oo

n---~ oo

= 0.

d,=l

Hence by Lebesgue's dominated convergence theorem,

Ai

lim Eo R Y/1 -

n~ oo

-- R(Yil)

Y

= O.

Hence lim n---~ 0t3

Varo( E s,( )-!s,(0,]) =0,

i = 1,2 . . . . . t.

(10)

Using (9), (10) and Markov's inequality,

S

limPo

1 d --1S(0

n~oo

n

Hence the result.

+--Xd

)e

=0.

"c

[]

The following theorem enables us to linearize the gradient vector of the dispersion function of repeated measures IBD.

Theorem 2. Suppose the assumptions (a)-(c) hold for the model (1) and the true value old is .40 = 0. For every e>Oandc>O, limPo { Supll~ll.
t>e } = 0 .

(11)

"~

The proof of Theorem 2 rests essentially on monotone components is S(d) and Theorem 1. Jureckova (1969, 1971) proved a similar result like Theorem 2 for the linear model with i.i.d, errors. Theorem 2 suggests how to construct a quadratic approximation to D(.). This result is given in Theorem 3.

Theorem 3. Suppose the assumptions (a)-(c) hold for the model (1) and the true value of d is d o = 0. For every e>Oandc>O lim Po~ Sup I[D.(d)-Qn(d)I[ ~ > e t = O,

n~oo

~.llJII ~< c

(12)

M.M.Rashid/ Statistics&ProbabilityLetters22 (1995)213-221

218

where

On(A)=~/r'~j~_li~=lnijl ~ k - ~ } - ~

fYiJ-~t

and

Qn(a) = On(0) +

d ' Z d - ~ n d'S(0).

(13)

The Qn(d) is a quadratic in d and has the property that Qn(O) = Dn(O) and the gradient of Qn(A) is a linear approximation to that of D . ( d ) .

Jaeckel (1972) proved Theorem 3 for the linear model with i.i.d, errors. Under the assumptions (a)-(c), Theorem 3 can be proved assuming (11) and then following the approach in Lemma 1 of Jaeckel (1972). By exploiting convexity of Dn(3) and without using the concordance conditions (3a-3c) of Jureckova (1971) on the design matrix Heiler and Willers (1982) proved that for the linear model with i.i.d, errors, the results given in (8), (11) and (12) are equivalent. Now, we state a similar result for the repeated measures IBD. 4. Suppose the assumptions (a)-(c) hold for the model (1), and true value of A is d ° = O. Then (8), (11) and (12) are equivalent.

Theorem

Proof.

Both D.(A) and Q.(A) are convex functions of A with D,(0) = Q.(0) and

O[D.(A)-On(d)]ad =__~-1 ( ~ n ) S d

+__~nnlS ( 0 ) - m z d .

Therefore, assuming (8) and following the approach in Lemma 3.1 of Heiler and Willers (1982) (12) follows. By Lemma 4.3 of Heiler and WiUers (1982), it is straightforward to show that (12) implies (11). Further, (11) implies (8), this implication is trivial. Hence (8), (11) and (12) are equivalent. [] 2.3. Asymptotic distribution of the rank estimators

Lemma 3. Let ~, - be a generalized inverse o r s satisfying S - S,Z- = E - , and suppose the assumptions (a)-(c) hold for the model (1). Then x / ~ ( a - a °) L M V N

O, -m- , ~ ' -

a s n---* zx3,

where a ° is the true value of a. P r o o f . L e t a ° = O.

Let

mn

Q(a) = O(O) + 2~z a'Za - a'S(O).

(14)

Suppose ~ and a minimize Q(a) and D(a) respectively. Then ~ = ( z / r ) Z - S ( O ) . Using Lemma 1, x / ~ t ~ MVN(O, (z2/m)2~-) and n ---, ~ . As a conclusion of Theorem 3, by letting A = ~v/-na,the minima of

M.M. Rashid / Statistics & Probability Letters 22 (1995) 213-221

219

Q(a) and O(a) coincide asymptotically (similar to T h e o r e m 3 of Jaeckel, 1972), i.e. x/~(fi - d) ~ 0 as n ~ oo. Hence x/~(a - 0) ~~ MVN(O, (zE/m)Z - ) and n ~ oo . []

3. R a n k tests and e s t i m a t i o n o f

In this section we develop three different statistics for testing hypothesis a = 0. We also discuss a simple consistent estimator of z.

3.1. Rank tests First we consider a test based on D(0) - D(a), the d r o p in the dispersion function. Suppose that 27- is a generalized inverse of 27 satisfying 27-27Z- = 27-. To test the hypothesis a --- 0 versus a # 0, we c o m p a r e D(fi) with D(0). The result is given in the following theorem. T h e o r e m 5. Suppose the assumptions (a)-(c) hold for the model (1) and the null hypothesis Ho: a = 0 is true.

Then, D* = 2 [D(0) - D(d)]/~ ~ X2_ 1

as n ~ oo

Proof. Since D(0) = Q(0), we have, O(O) - O(d) = [Q(0) - Q0i)] + I-Q(~) - Q(d)] + [Q(d)

- O(,i)'],

(15)

where Q(a) is defined in (14). The second difference in the right of (15), Q(~i) - Q(a), can be written as x//-n[d - a"J'[(1/x/~)S(0) - (z/2)m27x/~(d - 4)]. N o w using L e m m a 1 and x/~(d - ~i) ~ 0 as n ~ oo, we have Q(~) - Q(a) = %(1). Also from L e m m a 3, x/~d is b o u n d e d in probability. Hence by T h e o r e m 3, we assert that the Q(a) - D(a) = %(1). Therefore, D(0) - D(a) = Q(0) - Q(fi) + o9(1). Further it is easy to show that Q ( 0 ) - Q ( ~ ) = (1/2mn) zS'(0) 27-S(0). Therefore, 2 J - D ( 0 ) - D(~)]/z = (1/mn) S'(0) Z - S ( O ) + op(1). Using L e m m a 1, (1/mn) S'(O) Z-S(O) ~ Z2_1 as n ~ oo. Therefore,

2[D(0) -- D(d)]/z

~

•2_ 1

as n ~ oo

(16)

If we replace z by a consistent estimator ~ in (16), we still have an asymptotic chi-square distribution for large n, because the rank of 2: = t - 1. M c K e a n and H e t t m a n s p e r g e r (1976) developed D* for the linear model with i.i.d, errors. A second a p p r o a c h to testing Ho: a = 0 is based directly on the gradient vector evaluated at Ho. Therefore, using L e m m a 1, under the null hypothesis a = 0, S* =-1 S'(0)27-S(0) ~ Z 2 _ 1

as n ~ oo

(17)

r

because the rank of 27- is t - 1. The test statistic S* reduces to Durbin's (1951) statistic. A third a p p r o a c h to testing Ho: a = 0 is based directly on the R-estimator ~i. Using L e m m a 3, under the Ho, the quadratic

W = r~t'27d/'c2 ~ Z2_I

as

n ~ oo.

(18)

If we replace z by a consistent estimator ~ in (18), by Slutsky's theorem it will again asymptotically follow a chi-square distribution with t - 1 degrees of freedom as n ~ oo. []

220

M.M. Rashid / Statistics & Probability Letters 22 (1995) 213-221

3.2. Estimation o f z In Sections 2 and 3.1 we have seen that the p a r a m e t e r 272 1/[12(S~o~f (x, x ) d x ) 2] plays a p r o m i n e n t role in rank based inference for repeated measures I B D with exchangeable errors within subjects. Assuming without loss of generality that both treatment i and treatment i' are applied on the j t h subject ( j = 1,2 . . . . . 2,,), it can be shown that for repeated measures IBD, 1 / [ x / ~ z ] is the density of Z~ i'i') = Yi~ - Y~,j at ~i - 0~i,, for i -~ i' = 1, 2 . . . . , t. Since treatment i goes with treatment i'2ir times, Z] ~'e), Z~'i') ..... Z~;i, '~ are i.i.d, continuous and symmetric r a n d o m variables with median ~ - ~,. Therefore estimating z involves estimating the inverse of the density function at the median. W i t h o u t loss of generality let us assume that Z~ "~'~, Zti2'i') ..... Z~;~ ') are ordered observations. Then following Bloch and Gastwirth (1968), a consistent estimator of z, when 2,, is large (i.e. when n is large), is given by =

2,, =

-

F

+,-

Z[~]_~ ,]

i4:i'=

1,2 . . . . . t ,

(19)

I = O(2ai,) and 0 < d < 1. See Bloch and Gastwirth (1968) for a choice of l and d. We use 2/t(t - 1) S,~<~, s]i~ii! as a consistent estimator of z. Also 1 / x / ~ z , is the density ofe~ - e~,j (i :~ i' = 1, 2 . . . . . t; j = 1, 2 . . . . . n) at the median = 0. Therefore one where

can use pairwise differences of the residuals across the n subjects as a r a n d o m sample from a continuous p o p u l a t i o n with median = 0 and construct an estimator of z similar to (19) because &~- &~, converges to ~ - cq, as n --* oo. In that case an overall estimator of z will be based on the average of k(k - 1)/2 estimators.

4. Conclusion and further research In this p a p e r we have developed a robust alternative to n o r m a l theory inference for repeated measures I B D using Wilcoxon scores. O n e could use general score functions to define the dispersion function. Multiple c o m p a r i s o n s of treatment effects can be carried out by R-estimators &is.

Acknowledgements I would like to t h a n k the referee whose c o m m e n t s and suggestions were very helpful for improving quality of the paper. I a m also thankful to m y colleague Dr. Balgobin N a n d r a m for reading the revised version with patience.

References Arnold, S.F. (1981), The Theory of Linear Models and Multivariate analysis (Wiley, New York), p. 232. Bloch, D.A. and J.S. Gastwirth (1968), On a simple estimate of the reciprocal of the density function, Ann. Math. Statist. 39 (3), 1083-1085. Berry, D.W. (1987), The consultants forum: logarithmic transformation in ANOVA, Biometrics 43, 439456. Durbin, J. (1951), Incomplete blocks in ranking experiments, British J. Psychology 4, 85 90. Heiler, S. and R. Willers (1982), Asymptotic normality of R-estimators in the linear model, Selecta Statist. Canad. 6, 151-175. Jaeckel, L.A. (1972), Estimating regression coefficients by minimizing the dispersion of residuals, Ann. Math. Statist. 43, 1449-1458. John, J.A. (1987), Cyclic Designs (Chapman & Hall, London, ch. 8). John, P.W.M. (1971), Statistical Design and Analysis of Experiments (Macmillan, New York) p. 233, 244. Jureckova, J. (1969), Asymptotic linearity of a rank statistic in regression parameter, Ann. Math. Statist. 40, 1889-1900. Jureckova, J. (1971), Nonparametric estimate of regression coefficients, Ann. Math. Statist. 42, 1328-1338.

M.M. Rashid / Statistics & Probability Letters 22 (1995) 213-221

221

McKean, J.W. and T.P. Hettmansperger (1976), Tests of hypotheses based on ranks in the general linear model, Comm. Statist. Theory Methods 15 (8), 693-709. Nelder, J.A. and R. Mead (1965), A simplex method for function minimization, Comput. J. 7, 308-313. Rao, C.R. (1947a), General methods of analysis for incomplete block designs, J. Amer. Statist. Assoc. 42, 541-561. Samuels, M.L., G. Casella and G.P. McCabe (1991), Interpreting blocks and random factors, J. Amer. Statist. Assoc. 86, 798-821. Skillings, J.H. and A.G. Mack (1981), On the use of a Friedman-type statistic in balanced and unbalanced block designs, Technometrics 23, 171-177. Yates, F. (1940), The recovery of interbiock information in balanced incomplete block-designs, Ann. Eugenics 10, 317-325.