On the analysis of circular balanced crossover designs

journal of statistical •

Journal of Statistical Planning and Inference 69 (1998) 359-370

ELSEVIER

planning

and inTerence

On the analysis of circular balanced crossover designs J. Kunert* Fachbereich Statistik, Universitiit Dortmund, D44221 Dortmund, Germany Received 17 June 1996; revised 15 August 1997

Abstract

The paper considers Aza'is' (J. Roy. Statist. Soc. B, 49 (1987) 334-345) randomization procedure for circular balanced crossover designs. It is shown that this randomization does not justify the assumption of independent identically distributed errors when the estimates are corrected for carryover effects. This might lead to underestimation of the variance of treatment estimates. Similar to the results of Kunert (Biometrics, 43 (1987) 833-845) and Kunert and Utzig (J. Roy. Statist. Soc. B, 55 (1993) 919-927), we give constants, such that multiplication with this constant makes the usual estimate of variance conservative. © 1998 Elsevier Science B.V. All rights reserved. A M S classification: primary 62J10, 62K10; secondary 62P10 Keywords: Circular balanced crossover design; Direct effect; Carryover effect; Randomization

1. I n t r o d u c t i o n

In crossover designs the effects of treatments on experimental units are compared by applying a sequence of treatments to each unit. With these designs it is hoped that the differences between units can be treated as block effects. However, there are problems. First, there is concern that a treatment applied in period j might have a carryover effect on the response of the unit in periodj + 1. To reduce the impact of such carryover effects and to make the estimation of carryover effects possible, one often designs crossover experiments with a systematic order of the treatments instead of a randomized order. However, this brings the second problem with crossover designs. With a systematic ordering of the treatments, the usual randomization justification of identically distributed independent errors does not hold. This is

* E-mail: [email protected].

0378-3758/98/$19.00 © 1998 Elsevier Science B.V. All rights reserved. PII S 0 3 7 8 - 3 7 5 8 ( 9 7 ) 0 0 1 5 9 - 6

360

J. Kunert /Journal of Statistical Planning and Inference 69 (1998) 359-370

especially the case if one uses an analysis which takes account of the carryover effects; see Kunert (1987) or Kunert and Utzig (1993). Aza'is (1987) noted that the introduction of a preperiod, at which the same treatment is applied as at the last period, allows for some randomization. This randomization leaves the neighbour structure of the designs undisturbed. If carryover effects are neglected in the analysis, then Aza'is' randomization justifies the usual block model, with experimental units as blocks. The present paper considers the case that the estimates are corrected for carryover effects. It is shown that, even with Aza'is' (1987) randomization, the assumption of independent errors cannot then be justified. This is similar to the non-circular case as was considered by Kunert (1987) and Kunert and Utzig (1993). However, we can show that the usual estimate of the variance can be made conservative, if it is multiplied by a conservative constant. This is again similar to the results in Kunert (1987) and Kunert and Utzig (1993). However, in the circular case this constant is generally smaller.

2. The designs and the randomization Let there be n experimental units, each of which can be used in p periods. We have t treatments and assume t ~> p. Assume unit a receives treatment k in period q and that it has received treatment l in the period before. Then we assume that the response Yaq of unit a in period q satisfies Yaq = "rk + Pl + Uaq,

(1)

where Zk is the direct effect of treatment k, Pz is the carryover effect of treatment l and uaq is the theoretical response which we would get if there were no treatment effects. The treatments are allocated to the experimental units in the different periods through an experimental design and a randomization procedure (see e.g. Bailey, 1981). The randomization gives names to the unit/period combinations, the experimental design prescribes which treatment is assigned to each name. Aza'is' (1987) randomization acts as follows: Assume the experimental units are {al, ..., a,} and the periods are {ql . . . . ,qp}. Among the set of all permutations of the numbers 1. . . . ,n, we randomly select one permutation n, say. For each i, 1 ~< i ~< n, we independently select a random number rl E {0, 1, ... ,p - 1}. Then period q~ of unit ai gets the name (n(i), (j + ri)rood p). The design allocates treatment d(i,j) to the unit/period combination with the n a m e (i,j). Eq. (1) and the randomization imply that for the response y~j on the unit/period combination with the name (i,j) we have Yli = Zdli,j~ + Pd,,j- 1) + Vi~,

(2)

where v~j is a random variable, the distribution of which is a mixture of the distributions of the uaq and is determined by the randomization process. Note that there are two cases in which Eq. (2) does not follow from what we have said so far. The first case is when the name (i,j) is given to the first period ql of a unit. In this case we can make

J. Kunert/Journal o f Statistical Planning and Inference 69 (1998) 359-370

361

Eq. (2) true by introducing a preperiod. Before we start the p r o p e r experiment, we treat this unit with treatment d ( i , j - 1). N o t e that the n a m e (i,j - 1) is given to the last period qp of this unit by the randomization. The second case is period which receives the n a m e (i, 1). The r a n d o m i z a t i o n procedure implies that the preceding period received the n a m e (i, p) and Eq. (1) implies that Yil = zali, ~) + Pdti,pl + Vij"

Hence, we m a k e Eq. (2) true if we define d(i, O) = d(i,p). Let us consider the distribution of vii. The randomization implies that there are n u m b e r s ~o . . . . , ~Ep/21 and fl such that (i) all vii have the same expectation p, say, and the same variance cto. (ii) The covariance between vq and Vkl equals fl if/ # k and ~lJ tl i f / = k and IJ - l] ~< p/2. It equals ~p-lj-ll if IJ - II > p/2. Additionally, Aza'/s' (1987) r a n d o m i z a t i o n has one m o r e step: the labels of the treatments are randomized• We define Y = [ Y l l . . . . , Y l p , Y21 . . . . . Ynp] T, v = [vll . . . . . rip, v21 . . . . . v,p] T, z = [r~ . . . . . zt] T and p = l-p1 . . . . . &]T. F o r every i, 1 ~< i ~< n, let Ti be the treatment design matrix for the observations on unit i. Then T = [ T t x, . . . , T,X] T is the treatment design matrix of the whole design• If F is the carryover effect design matrix, then F can be written as F = [F1x, . . . , F,x] T, where Fi = V T i , 1 <<. i <<. n, and

V=

[i00]0 0 1

"..

• .-

0

"

1

.

0

If we define In as the n-vector of ones and In as the n x n identity matrix and if ' ® ' denotes the K r o n e c k e r product of matrices, then model (2) in vector notation becomes Y=

T z + F p + v,

where E(v) = #lnp and Cov v = In ® Z + fllnplnXp. Here _r is a circulant matrix, such that "~0

~1

~2

~3

• ..

0{2

~1

~0

~1

~2

"'"

~3

~2

~2

~1

~0

~1

O~0

O{1

O~2

0~1

~0

Z =

- I lv IT

(X2

~X3

~1

~2

362

J. Kunert / Journal of Statistical Planning and Inference 69 (1998) 359-370

We call the design a balanced circular design if (B1) the design is a balanced incomplete block design with experimental units as blocks, (B2) each treatment is preceded by every other treatment equally often, but never by itself We also define the matrix B = I, ® lp. With this notation (B1) implies that T T T = F T F -=- np It t

and TXBBTT

= TTBBTF

= FTBB~F

np - t ( t : 1)((p - 1)lt 1~ + (t - p)I,).

Condition (B2) implies that np T X F - t ( t ~ 1)(1, l,v - I,).

Circular balanced designs and their construction are frequently treated in the literature; see e.g. Lawless (1971). Aza'is, et al. (1993) present some examples and constructions of circular balanced designs.

3. The simple block model Azais (1987) used his randomization to validate the analysis of balanced circular designs in the simple block model Yij = Zd(i,j) ~- Yi ~- % ,

(3)

where the %, 1 ~< i ~ n, 1 ~
J. Kunert / Journal o f Statistical Planning and Inference 69 (1998) 359 370

363

For a balanced block design we have

Defining Y = pr±(B)Y and for a design satisfying (B1) we hence get t-1 ~A(simple) -- • ( p

- - 1) T T y "

Unfortunately, the

estimate lT~(simple)is not unbiassed under model (2). We have

t-1 TA El "['(simple)-n(p -- 1) ITT x E ~

t--1 - n ( p - 1) Ix{Txpr±(B)Tz + TTpr±(B)Fp} t-1 x = lxz -~ n[p-__-'l)l TTpr±(B)Fp. Hence,

b(l,p) = n~p--ll)1TTTprI(B)Fp is the bias

of/T~lsimple)in the presence of p. Define

b* = max] b(/, p)1, l,p

where the maximum is taken over all l and p for which ITl = pTp = 1. Aza'is and Druilhet (1997) have shown that TTpr±(B)F is completely symmetric and, among all designs satisfying (B1), has minimal trace if the design satisfies (B2), that is if it is circular balanced. They point out that this implies that the circular balanced design minimizes ~o(TTpr±(B)F) for all q~ satisfying (a) q~ is invariant to simultaneous permutation of rows and columns, (b) ~o(aTTpr±(B)F) is increasing in a > 0, (c) q) is convex. It is obvious, that b* satisfies (a) and (b). To see that it is also convex, note that

max liT(am + (1 - a)N)p] = max ltT(am)p + IT((1 -- a)N)p] <~amaxllT Mp] + (1 -- a)maxllT Np[ for any two matrices M and N. Hence, it follows that the circular balanced design minimizes the maximum bias b*. Note that b* is a different concept from the criterion M B used by Aza'is and Druilhet (1997), which was also considered by Kunert (1985). However, it can be shown that both are multiples of the largest singular value of TTpr±(B)F, so b* and M B are equivalent to each other.

364

J. Kunert/Journal of Statistical Planning and Inference 69 (1998) 359-370

4. The analysis allowing for carryover effects A simple analysis which allows for carryover effects would extend model (3) to Ylj = Zd(~,j)+ Pd(i,j- 1) 3V 7i ~- eij,

(4)

where Zd(i,~, 7i and eij are as in (3) and Pd(i,i- lj is the carryover effect as in (2). It was shown by Kunert (1984) in the special instance where p = t, that circular balanced designs are universally optimal for the estimation of direct effects in model (4) over all designs with the same number of treatments, periods and experimental units. Druilhet (1997) has shown that they are also universally optimal ifp < t. Unfortunately, as we will see now, the randomization model (2) cannot validate the variance estimation from model (4). It turns out, however, that the approach of Kunert and Utzig (1993) can be used for the circular model as well. 4.1. The variance of the estimate which is corrected for carryover The best linear unbiased estimate of a contrast lTz in model (4) is lT'~

......

~c. . . . . ted = (T T pr±([B, F]) T)+ T T pr±([B, F]) Y. We have TTprl([B, F]) = T x pr±(B) - T x prl(B)F(F T pr±(B)F)+ F x prl(B). Since for a circular balanced design TTpr±(B)F_

t - n 1 ( l t - t 11 tiT )

and

it follows that T T pr±([B, F]) T = ( t - 1 ) ( p - 1 )

(p-2)

(1)

It-tltl~

Consequently, (t--1)(p--1)

I

TT)7+

i

FT17

t

.

ted where

.L Kunert / Journal of Statistical Planning and Inference 69 (1998) 359 370

365

It is easy to show that l T^ z . . . . . . ted is unbiassed in model (2). The covariance matrix of . . . . . . ted in model (2) equals Cov(? ...... t~a) =

(t -- 1)2(p -- 1)2(

~£pT~Z-2)Y

x

=

TT +

1

)

FT Pr±(B)C°v(Y)Pr±(B)

p-1

( ') T+

F

p-1

( t - 1 ) 2 ( p - 1)2( ~ 2 p 2 ~ - - ~)~

TT +

1

FT

p-1

) (') (in@S)

T+

p-1

F

where

N o t e that S has row and column sums zero but has the same circulant structure as Z. We denote the (l,1)-element of S by s and the (1,p)-element of S by a. N o t e that all diagonal elements of S are equal to s, while all elements of the first sub-diagonal are equal to a. The randomization of treatment labels implies that we can write C o v ~ ...... tea = t r ( C o v t . . . . . .

ted)

(1) It -

t-1

1, 1T

7

"

We have

t r T T ( I . ® S ) T = ~ tr(T'X, STi). i=l

Each row of Ti contains exactly one 1, all other elements are zero. Hence, we can write Ti = [Ip,O]IIi, where Hi is a t × t permutation matrix. This implies that tr T~STi = t r S = ps and tr(TT(I, ® S)T) = rips. Similarly, tr(FT(I, ® S)F) =nps and tr(TT(I. ® S)F) = n tr(SV) = npa. In all, tr(Cov'f~ . . . . . t~d)

(t -- 1)2(P -- 1)2 (

= n2p2(p _ 2) 2

2

(t--1)2(p--1)2( l + ( p - 1 ) 2 ~-p~ (t-

2) 2

1)2

1

nps+ p-- 1 npa +--npS(p_1)2

(p _ 1)2

2 s + p_l

-- n ~ p ---~-'2 ((1 + (p -- 1)2)s + 2(p -- 1)a). z)

t a

)

366

J. Kunert / Journal o f Statistical Planning and Inference 69 (1998) 359-370

It follows that the estimate I T"~corrected A has variance var IT'~c..... ted = tr(Cov~ ...... ted) IT I t--1 l(t-

1)((1 + ( p -

4.2.

1)2)s + 2 ( p -

1)a)lT1"

p ( p -- 2) 2

n

The expectation of the estimated variance

In model (4) we estimate a 2 through ~2

1

yT prJ_([B, F, T ] ) Y,

=

np -- 2t - n + 2

and the variance of/T'~c . . . . .

ted

through

vfir (lT'~c..... ted) = ~2 (t -- 1)(p -- 1)1Tl. n p ( p - 2) If model (4) were valid, then the expectation of vfir(lTic ..... t~a) in model (2) should be equal to var(lV¢c ..... tea). Note that in model (2) we have Epr-:([B,r,

T ] ) Y = 0,

hence E Y T p r ' ( [ B , F, T ] ) Y

= tr {pr±([B, F, T])Cov Y} = tr { p r l ( B ) C o v Y - p r ± ( B ) F ( F T pr±(B)F) + F T p r ± ( B ) C o v Y -

pr(pr'([BF])T)Cov Y }

=tr(I.®S)

/ = nps~l

t -- 1 T(I.®S)F n~-p---1)trF

t--1 "~ n-(p----1))

np(p (t-

2)

1 ) ( p - 1)

tr{Cov(¢......

n p ( p - - 2) . . . . . (t ~ 1)---(p--1) t r t " ° v t z . . . . . . ted'.

It follows that V~(/Tf¢ ..... ted) has the expectation EV~(/T¢c ..... ted) -- (t -- 1)(p -- 1)ITI 1 n p ( p -- 2) np -- 2 t -- n + 2 E Y T p r ± ' ' B[[' F ' (t -- 1)(np -- n -- t + 1) s . lT1 = (p - - 2 ) n ( n p -- n - - 2 t + 2)

t--1 np - - 2 t -- n + 2

var(IT¢¢..... ,~d).

T])Y

ted)}

J. Kunert / Journal of Statistical Planning and Inference 69 (1998) 359-370

367

Define var(/T'~ . . . . .

ted)

Q = E va'r(IT'~ c. . . . .

ted)"

Then model (4) is valid in the r a n d o m i z a t i o n model (2), if for any X, i.e. for any s and a, the quotient Q is 1. We calculate 1 ( (np--n--t + 1)p(p--2) ) 1/Q = (np _ n _ 2t + 2) \(pZ _ 2p + 2) + 2(p _ l)a/s - t + 1_. It follows that Q = 1 if and only if a = - s/(p - 1). This shows that model (4) in general is not validated by the randomization. N o t e that a is the covariance of directly adjacent entries of Y. Hence, it is likely in practice that a > - s/(p - 1).

4.3. Upper bounds for the bias o f the estimated variance N o t e Q is increasing in a/s. So in practice Q is liable to be greater than 1. However, if we can determine an upper b o u n d for a/s, we get an upper b o u n d for Q. As in Kunert and Utzig (1993), we hence can give upper bounds for Q, depending on n, t and p. Multiplication of vS,r ( l r 0 by this b o u n d Q*(t,n,p), say, will m a k e the analysis conservative. F o r details, see K u n e r t and Utzig (1993). Case 1: p = 3. Then we have

S=

s

.

a

D u e to the fact that S has row and column sums zero, a must be equal to - s/2. Hence, Q*(t,n, 3) = 1 for all n and all t >/3. So the analysis in model (4) is valid for p=3. Case 2: p = 4. In this case

S=

s

a

a

s

b

a

Since S has row and column sums zero, it holds that b = - (s + 2a). Hence, a cannot be greater than 0; otherwise b would be less than - s and S would not be a nonnegative-definite matrix. It follows that

1 1 ( 8 ( 3 n Q>~3n - 2 t + 2

- t +1) ]-6

with equality holding if a = 0, b = - s.

+ l ) 1 =2 n

- 9t + 15n-

1 0 t + 10"

J. Kunert / Journal of Statistical Planning and Inference 69 (1998) 359 370

368

H e n c e , Q can be g r e a t e r t h a n one, b u t there is a n u p p e r b o u n d 15n-10t+10 12n-9t+9

Q <<.Q * ( t , n , 4 ) =

Case 3: p = 5. In t h a t case

bS S=

a b b

ab s a a s b a

.

w h e r e b = - ½ (s + 2a). T h e n o n z e r o e i g e n v a l u e s of S are (5s + x/-5(s + 4a))/4

(twice)

(5s - x f 5 ( s + 4a))/4

(twice).

and

It follows t h a t S is n o n n e g a t i v e definite o n l y if

a/s <<.( x f 5 - 1)/4 = a*

say.

Hence, 30n - (t - 1)(16 + 4a*)

1

34n + 16Ha* - ( t -

1)(17 + 8a*)

and 34n + 4(,,/-5 - 1)n - (t - 1)(15 + 2x/5)

Q <~ Q*(t,n, 5 ) =

-

(t -

I)(15

+

Case 4: p = 6. T h e n S is of the f o r m

S=

[iacbil a b c b

s a b c

a s a b

b a s a

w h e r e c = - (s + 2a + 2b). T h e n o n z e r o e i g e n v a l u e s of S are 2s+3a+b -- 3 a -

(twice) 3b

(twice)

J. Kunert/Journal of Statistical Planning and Inference 69 (1998) 359-370

369

and 2s + 4b

(simple).

It follows t h a t S can only be n o n n e g a t i v e definite if b >>. - s/2 and, consequently, if a/s <~1. This implies t h a t 1

120n

-

1)55

(t-

/> 155n - ( t -

1)62

and

Q <~ Q*(t,n, 6) =

155n - 62(t - 1) 120n - 55(t - 1)"

Case 5: p >/7. Here we use an u p p e r b o u n d which in general is n o t sharp. N o t e that S can only be n o n n e g a t i v e definite if a ~< s. This implies that 1

n(p -

>~

1)(p

-

2) -

(t -

1)2(p

-

1)

n p ( p - 1) - 2(t - 1)p

and Q ~< Q* (t, n, p) =

np(p-

1 ) - 2(t - 1)p

n(p - 1)(p - 2) - 2(t - 1)(p - 1)"

4.4. Discussion o f the bounds F o r every t, n, a n d p we have f o u n d a b o u n d Q*(t, n,p) for Q. T o see h o w large these b o u n d s are, we list some values in T a b l e 1. N o t e t h a t a circular b a l a n c e d design can Table 1 Some values of the conservative multiplier p

2

t=4

t=5

t=6

t=7

t=8

t=9

t=10

4

1 2 3 oo

1.67 1.33 1.30 1.25

1.46 1.31 1.28 1.25

1.39 1.30 1.28 1.25

1.35 1.29 1.27 1.25

1.33 1.28 1.27 1.25

1.32 1.28 1.27 1.25

1.31 1.27 1.27 1.25

5

1 2 3 3c

1.53 1.37 1.34 1.30

1.45 1.35 1.33 1.30

1.42 1.34 1.32 1.30

1.39 1.33 1.32 1.30

1.38 1.33 1.32 1.30

1.37 1.33 1.32 1.30

6

1 2 3 oo

1.43 1.34 1.32 1.29

1.40 1.33 1.32 1.29

1.38 1.33 1.31 1.29

1.36 1.32 1.31 1.29

1.35 1.32 1.31 1.29

7

1 2 3 oc

1.56 1.46 1.44 1.40

1.53 1.45 1.43 1.40

1.51 1.44 1.43 1.40

1.49 1.44 1.42 1.40

J. Kunert / Journal of Statistical Planning and Inference 69 (1998) 359-370

370 only exist if

2-

np t(t

1)

-

is a n integer. Hence, we can restrict a t t e n t i o n to such n for which n =

2 t(t -- 1)

,

2e[~.

P Note that not all such n are integers. O f course, a circular b a l a n c e d design can only exist, if n is an integer. So we do not get designs for every entry of Table 1. However, Table 1 shows that the b o u n d is decreasing in 2. If a given design has t, n a n d p such that 2 > 3, then Q*(t, n,p) is between the entries of Table 1 for 2 = 3 a n d for 2 ~ 0o. F o r p >/7, if we leave t a n d 2 fixed, then the b o u n d Q * ( t , 2 ( t ( t - 1)/p),p) is decreasing in p. Therefore, we only list the cases p = 4, 5, 6 a n d 7. If we c o m p a r e the entries of T a b l e 1 to the c o r r e s p o n d i n g values in the table of K u n e r t a n d Utzig (1993), we find that they are generally smaller. This is due to the different r a n d o m i z a t i o n . Aza'is' (1987) r a n d o m i z a t i o n procedure forces more structure on the covariance matrix S. This allows better b o u n d s for Q. However, it does n o t make the b o u n d s equal to 1. M o n o d (1996) gives examples of designs where a r a n d o m i z a t i o n similar to Azais' (1987) procedure justifies model (4). However, these designs need pairs of consecutive identical treatments a n d p > t. The present paper shows that his results c a n n o t be extended to the case of circular balanced designs.

References Azais, J.-M., 1987. Design of experimentsfor studying intergenotypic competition. J. Roy. Statist. Soc. Ser. B 49, 334~345. Azais, J.-M., Bailey, R.A., Monod, H., 1993. A catalogue of efficient neighbour-designswith border plots. Biometrics 49, 1252-1261. Azais, J.-M., Druilhet, P., 1997. Optimality of neighbour balanced designs when neighbour effects are neglected. J. Statist. Plann. Inference, to appear. Bailey, R.A., 1981. A unified approach to design of experiments.J. Roy. Statist. Soc. Ser. A 144, 214~233. Druilhet, P., 1997. Optimality of circular neighbor balanced designs., to appear. Kunert, J., 1984. Designs balanced for circular residual effects. Commun. Statist. Theory Methods 13, 2665 2671. Kunert, J., 1985. Optimal experimental design when the errors are assumed to be correlated. Statist. Decisions 2 (Suppl.), 287 298. Kunert, J., 1987. On variance estimation in crossover designs. Biometrics 43, 833 845. Kunert, J., Utzig, P., 1993. Estimation of variance in crossover designs. J. Roy. Statist. Soc. Set. B 55, 919 927. Lawless, J.E., 1971. A note on certain types of BIBD's balanced for residual effects. Ann. Math. Statist. 42, 1439 1441. Monod, H., 1996. Valid randomization for models with neighbour effects from treatments. Statistics 27, 311 324.