On efficient robust first order rotatable designs with autocorrelated error

Journal of the Korean Statistical Society 37 (2008) 95–106 www.elsevier.com/locate/jkss

On efficient robust first order rotatable designs with autocorrelated errorI Rabindra Nath Das a , Sung Hyun Park b,∗ a Burdwan University, Burdwan, WB, India b Department of Statistics, Seoul National University, Seoul 151-747, Republic of Korea

Received 1 March 2007; accepted 1 August 2007 Available online 16 April 2008

Abstract Generally, it is very difficult to derive optimal or at least efficient designs for linear models with correlated observations, and for some correlation structure, an exact D-optimal design does not exist. In this paper we have developed the notion of a D-optimal robust first order design (D-ORFOD) for linear model with a general correlated error structure. We have shown that D-optimal robust first order designs are always robust first order rotatable designs (RFORDs) but the converse is not always true. For a first order linear model with autocorrelated error, we have developed a set of efficient RFORDs with efficiency around ninety percent and the developed designs are very close to D-ORFODs. We have also developed a new method of analysis that is the estimation of regression parameters, correlation parameter and error variance, assuming the correlation parameter involved in the correlation structure is unknown. c 2008 The Korean Statistical Society. Published by Elsevier Ltd. All rights reserved.

AMS 2000 Subject Classification: primary 62K20; secondary 05B30 Keywords: Autocorrelated structure; D-optimal; Optimum design; Robustness; Rotatability

1. Introduction In general, it is difficult to derive optimal or at least efficient designs for linear models with correlated observations. Recently, Bischoff (1995) has proved two determinant formulas. These formulas are useful for deriving optimal or efficient designs with respect to the D-criterion. Bischoff (1996) pointed out that a D-optimal design does not exist for a linear model with some correlated errors, and for the tri-diagonal structure he derived maximin designs. Therefore, Kiefer and Wynn (1981) suggested seeking optimal designs only in the class of uncorrelated and homoscedastic errors. For special factorial linear models with correlated observations, exact optimal designs are known, see for example, Kiefer and Wynn (1981, 1984), and the references cited there. Gennings, Chinchilli, and Carter (1989), studied optimality under the non-linear model with correlated observations. Only a little is known about optimal designs of linear regression models when the observations are correlated, see Bischoff (1992, 1995, 1996) and related references cited there. I This work was partially supported by the second stage Brain Korea 21 Project. ∗ Corresponding author.

E-mail addresses: [email protected] (R.N. Das), [email protected] (S.H. Park). c 2008 The Korean Statistical Society. Published by Elsevier Ltd. All rights reserved. 1226-3192/$ - see front matter doi:10.1016/j.jkss.2007.08.003

96

R.N. Das, S.H. Park / Journal of the Korean Statistical Society 37 (2008) 95–106

Herein we have considered the autocorrelated structure, which is mostly used in some econometrics model and is given in Section 3. According to Bischoff (1995), the exact D-optimal does not exist for an autocorrelated error structure. In this paper we have developed robust rotatable designs (in Section 3.2) having efficiency greater than ninety percent. Bischoff has not introduced D-optimal robust designs and he has not established any relation between D-optimal robust designs and robust rotatable designs. Again Bischoff suggested the estimation of parameters by the ordinary least squares (OLS) method, which is not always appropriate under the correlated model. Herein we have developed a new method for estimation of parameters which gives the best linear unbiased estimators of all the parameters except the intercept, which is often unimportant in practice. In general, a linear model with explanatory variables is a response surface model. A short review of response surface methodology (RSM) is given in Box and Hunter (1957), Hader and Park (1978), Myers, Khuri, and And Carter (1989), and Park (1987). Recent books on RSM have been written by Box and Draper (1987), Khuri and Cornell (1996), Myers and Montgomery (2002), Pukelsheim (1993) etc. For correlated observations, rotatability for first and second order designs are given in Das (1997, 2003a,b, 2004), Das and Park (2006) and Panda and Das (1994). In this paper, we have considered the first order linear regression model with correlated observations (in Section 2.1). The main objective of this paper is to develop efficient designs for first order linear models with correlated errors. Our approach is to develop designs which satisfy approximately the conditions of D-optimality, free of the correlation parameter involved in the correlation structure. In this sense our method is completely different from the earlier concept of Kiefer and Wynn (1981). Herein, we have derived the necessary and sufficient conditions for first order D-optimality when errors in observations have a general correlated error structure (in Section 2.3). Then we have established a relation between the first order D-optimal designs and first order rotatable designs, when errors in observations are correlated (in Section 2.3). Then we have specially considered the autocorrelated error structure, which most frequently occurs in the observations under the linear model (in Section 3). In this situation, we have examined the conditions of first order D-optimality (in Section 3.1). We have developed some efficient RFORDs which are nearly D-optimal robust first order designs with autocorrelated error (in Section 3.2). In the process we have developed a new method of analysis that is the estimation of regression parameters, correlation parameter and error variance assuming the correlation parameter involved in the correlation structure is unknown (in Appendix). Concluding remarks is given in Section 4. 2. First order correlated linear model 2.1. Model Suppose there are ‘k’ factors x = (x1 , x2 , . . . , xk )0 which yield a response of yu on the study variable y when ∼

x = x0u = (x1u , x2u , . . . , xku )0 , 1 ≤ u ≤ N . Assuming that the response surface is of first order, we adopt the model

∼

∼

yu = β0 +

k X

βi xiu + eu ; 1 ≤ u ≤ N ,

i=1

or,

Y = X β + e,

∼

∼

(2.1)

∼

where Y is the vector of recorded observations on the study variable y, β = (β0 , β1 , β2 , . . . , βk )0 , is the vector of ∼

∼

regression coefficients, X = ( 1 : (xiu ); 1 ≤ i ≤ k, 1 ≤ u ≤ N ) is the design matrix. Further, e is the vector of errors ∼

∼

which are as assumed to be normally distributed with E(e ) = 0 and D(e ) = W with rank (W ) = N . The matrix ∼

∼

∼

W may represent any structure with correlated errors, e.g., tri-diagonal, autocorrelated structure, etc. In general, the matrix W is unknown but for all the calculations here, W is assumed to be known. In practice, W includes a number of unknown parameters, and in the calculations which follow, the expressions for W and W −1 are replaced by those obtained by replacing the unknown parameters by suitable estimates or some assumed values. 2.2. Analysis We will assume X 0 W −1 X is positive definite for a known W . The best linear unbiased estimator (BLUE) of β , for ∼

a known W is βˆ = (X 0 W −1 X )−1 (X 0 W −1 Y ),

∼

∼

(2.2)

97

R.N. Das, S.H. Park / Journal of the Korean Statistical Society 37 (2008) 95–106

with D(βˆ ) = (X 0 W −1 X )−1 = ((vi j ))−1 = ((v i j )), say 0 ≤ i, j ≤ k, where v00 = 1 0 W −1 1 , v0 j = 1 0 W −1 x j ; 1 ≤ ∼

∼

j ≤ k, vi j = xi

0 W −1

∼

x j and xi = (xi1 , xi2 , . . . , xi N ∼

∼

)0 ;

∼

∼

∼

1 ≤ i, j ≤ k.

Originally, the matrix W is unknown but for calculations as in (2.2), the expressions for W and W −1 are replaced by those obtained by replacing the unknown parameter ρ ∈ W by its suitable estimate or some assumed values. An analysis is given in the Appendix, when ρ ∈ W is replaced by some assumed values. 2.3. D-optimal first order regression designs In this section, we are interested in finding regression designs under the correlated regression model (2.1) such that variance of estimated response at x = (x1 , x2 , . . . , xk )0 is optimum in a certain sense, for a wide variety of values for ∼ the dispersion matrix W . The goals of the four optimality criteria are D-criterion goal → minimize |(X 0 W −1 X )−1 |, or equivalently maximize |(X 0 W −1 X )|, A-criterion goal → minimize trace[(X 0 W −1 X )−1 ], G-criterion goal → minimize maxx∈χ [N f 0 (x)(X 0 W −1 X )−1 f (x)], I V -criterion goal → minimize average[N f 0 (x)(X 0 W −1 X )−1 f (x)], over x ∈ χ , where x is any point in the design region χ , N is the design size, and f (x) = ( f 1 (x), . . . , f p (x))0 is a vector of p real valued functions based on the p model terms, W is the known variance covariance matrix of errors. In general, the dispersion matrix W is unknown, it is very difficult to find optimal designs under any one of the above four criteria. In this paper we have established a relation between D-optimal designs and rotatable designs for correlated observations. It is also established herein that D-optimal designs are a sub-class of robust rotatable designs. So our objective is to develop efficient robust rotatable designs which approximately satisfy the conditions of D-optimality, free of the correlation parameter involved in the correlation structure W . Thus our goal is to develop efficient designs through robust rotatable designs. Below we have gone through the first order D-optimality criterion with correlated observations. In this connection we give the following definitions. Definition 2.1 (Rotatable Design). In RSM rotatability is a natural and highly desirable property of regression designs. A design is said to be rotatable if the variance of the estimated response at a point is only a function of distance from the design centre to that point. Definition 2.2 (Robust First Order Rotatable Design). A design D on k factors under the correlated model (2.1) which remains first order rotatable for all the variance–covariance matrices belonging to a well defined class W0 = {W positive definite: W N ×N defined by a particular correlation structure possessing a definite pattern} is called a robust first order rotatable design (RFORD), with reference to the variance–covariance class W0 . D-optimality: We have assumed the normality of errors. Then in the linear context of the model (2.1), a confidence ellipsoid for β of a given confidence coefficient and for a given residual sum of squares, arising from the design matrix ∼

X when W is known, has the form {β : (β −βˆ )0 (X 0 W −1 X )(β −βˆ ) ≤ constant}, ∼

∼

∼

∼

∼

where βˆ is the least squares estimate of β . The content of this ellipsoid is proportional to {det. M(x)}−1/2 , where ∼

∼

M(x) = (X 0 W −1 X ). A natural design criterion, known as D-criterion, is that of making this ellipsoid as small as possible, in other words that of maximizing the determinant M(x). Definition 2.3 (D-optimal Robust First Order Design). A first order regression design D on k factors in the correlated model (2.1) will be D-optimal robust first order design (D-ORFOD) if the determinant (X 0 W −1 X ) is maximum (over the design space χ = |xiu | ≤ 1; 1 ≤ i ≤ k, 1 ≤ u ≤ N ) for all the variance–covariance matrices belonging to a well defined class W0 , as in Definition 2.2, with reference to the variance–covariance class W0 . Note that |M(x)| ≤

k Y i=0

vii ,

(2.3)

98

R.N. Das, S.H. Park / Journal of the Korean Statistical Society 37 (2008) 95–106

and equality holds if and only if v0i = 0; 1 ≤ i ≤ k and vi j = 0; 1 ≤ i 6= j ≤ k, and vi j ’s are as in (2.2). Note that the inequality in (2.3) is known as Hadamard’s inequality (Anderson, 1984, p. 54), which can be easily proved as −1 |M(x)| = |M11 − M12 M22 M21 |.|M22 | ≤ v00 |M22 |, where M22 is k × k and applying induction on |M22 |. Hence we have the following theorem. Theorem 2.1. The necessary and sufficient conditions for a D-optimal robust first order design in the model (2.1) are (i) v0 j = 1 0 W −1 x j = 0; 1 ≤ j ≤ k, for all ρ, ∼

∼

(ii) vi j = xi 0 W −1 x j = 0; 1 ≤ i 6= j ≤ k, for all ρ, ∼

∼

and (iii) vii = xi 0 W −1 xi = maximum possible value (constant) = µ, say; ∼

∼

1 ≤ i ≤ k, for all ρ, over the design space χ , as in Definition 2.3,

(2.4)

where ρ is the correlation parameter in W ∈ W0 , as in Definition 2.2. Proof directly follows from (2.3) and Definition 2.3. The variance of the estimated response yˆx at x = (x1 , x2 , . . . , xk )0 of a first order D-optimal design with correlated ∼ error is given by k 1X 1 + x 2 = f (r 2 ), say, v00 µ i=1 i

V ( yˆ x ) =

(2.5)

Pk

= 1 0 W −1 1 and µ as in (2.4). ∼ ∼ Following (2.5), we have the following theorem.

where r 2 =

2 i=1 x i , v00

Theorem 2.2. A D-optimal robust first order regression design is always a robust first order rotatable design but the converse is not always true. Proof directly follows from (2.4) and (2.5). If vii is constant = λ 6= µ, 1 ≤ i ≤ k, where µ is as in (2.4), then the design is a robust first order rotatable but not a D-optimal robust first order design. Let V ( yˆ x )0 be the variance of the estimated response yˆx at x of a D-ORFOD ‘d0 ’ as defined in Definition 2.3 and ∼

V ( yˆ x ) be the variance of the estimated response yˆx at x of any other RFORD ‘d’ (possibly near D-optimal), under ∼ Pk the model (2.1). Note that V ( yˆ x )0 and V ( yˆ x ) are both functions of r 2 = i=1 xi2 and the correlation parameter or parameters dictated by the correlation structure class W0 in Definition 2.2. The efficiency ratio of the given RFORD ‘d’ is denoted by “ER” and is naturally given by Max. {V ( yˆ x )0 }

ER% =

0≤r 2 ≤k

Max. {V ( yˆ x )}

× 100,

(2.6)

0≤r 2 ≤k

and “ER” is computed for the class of all covariance matrices in Definition 2.2, i.e., in the well defined range of correlation parameter or parameters as dictated by the class. Here V ( yˆ x )0 may denote the variance of yˆ x for the hypothetical D-optimal design which may not be possible to arrive at explicitly. Note that an exact D-optimal design does not exist for some correlation structures (see Bischoff (1995)), but we can derive the V ( yˆ x )0 from (2.5). Below we will study the necessary and sufficient conditions for D-optimality of first order regression designs, the variance function and the corresponding robust D-optimal and nearly D-optimal robust first order regression designs under the autocorrelated correlation structure of errors commonly encountered in practice. 3. Autocorrelated structure In this section we will study D-optimal robust first order regression designs when errors have autocorrelated structure given by W0 = {D(e ) = [σ 2 {ρ |i− j| }1≤i, j≤N ] = W N ×N (ρ)}. ∼

(3.1)

R.N. Das, S.H. Park / Journal of the Korean Statistical Society 37 (2008) 95–106

99

This structure is very well known and commonly used in Econometrics (see Kiefer and Wynn (1981, 1984)). First order autoregressive AR(1) process in errors gives rise to an autocorrelated structure as described below. The AR(1) process is defined as et = ρet−1 + t where, E( ) = 0 and E(  0 ) = σ2 I, i.e., t ’s are uncorrelated and have ∼

∼∼

∼

equal variances. The {et } forms an autoregressive process with the autocorrelated structure for the variance covariance matrix. Use of the autoregressive AR(1) process for errors in regression problems in the usual normal set up is very widely known. Observe that, W N−1×N (ρ) = {σ 2 (1 − ρ 2 )}−1 [(1 + ρ 2 )I N − ρ 2 A0 − ρ B0 ],

(3.2)

where I N is the N × N identity matrix, A0 is the N × N matrix with elements a11 = a N N = 1 and all other elements are 0 (zeros), and B0 is the N × N matrix with bi j = 1 for |i − j| = 1 and all other elements 0 (zeros). Following (2.4), the necessary and sufficient conditions for a D-optimal robust first order regression design under the correlation structure (3.1) simplify to (1) v0 j =

N X

x ju − ρ

u=1

(2) vi j =

N X

N −1 X

x ju = 0; N −1 X

xiu x ju + ρ 2

u=1

(3) vii

1 ≤ j ≤ k, for all ρ,

u=2

( xiu x ju − ρ

u=2

N −1 X

xiu x j (u+1) +

u=1

N −1 X

) xi(u+1) x ju

= 0;

u=1

1 ≤ i, j ≤ k, i 6= j, for all ρ, " # N N −1 N −1 X X X 2 2 2 2 2 −1 = {σ (1 − ρ )} xiu + ρ xiu − 2ρ xiu xi(u+1) = constant = µ;

u=1

1 ≤ i ≤ k, for all ρ,

u=2

u=1

(3.3)

where ρ ∈ W N ×N (ρ), as in (3.1). 3.1. D-ORFODs with autocorrelated structure In this subsection we find the necessary and sufficient conditions for D-optimal robust first order regression designs under the structure (3.1). The following theorem gives the necessary and sufficient conditions for D-ORFOD, as given in Definition 2.3 under the structure (3.1). Theorem 3.1. The necessary and sufficient conditions for a D-optimal robust first order regression design under the correlation structure (3.1), over the design space χ as in Definition 2.3, are (i) v0 j = 0;

1 ≤ j ≤ k, for all ρ,

(ii) vi j = 0;

1 ≤ i, j ≤ k, i 6= j, for all ρ,

(iii)

N X

2 xiu = N;

u=1

(iv)

N X u=1

2 xiu = N;

N −1 X

2 xiu = N − 2;

N −1 X

u=2

u=1

N −1 X

N −1 X

2 xiu = N − 2;

u=2

xiu xi(u+1) = −(N − 1);

xiu xi(u+1) = (N − 1);

1 ≤ i ≤ k, if ρ > 0, 1 ≤ i ≤ k, if ρ < 0,

(3.4)

u=1

where ρ ∈ W N ×N (ρ) and v0 j , vi j are as in (3.3). Note that, in most the practical situations, ρ is positive for the autocorrelated structure (3.1). Here we study DORFOD for positive and negative values of ρ. Following (2.5) and using (3.4), the variance function of a D-ORFOD under the structure (3.1), over the factor space as in Definition 2.3, is   1+ρ (1 − ρ 2 )r 2 V ( yˆ x )01 = σ 2 + , if ρ > 0, (3.5) N − (N − 2)ρ N + (N − 2)ρ 2 + 2ρ(N − 1)   (1 − ρ 2 )r 2 1+ρ 2 + and V ( yˆ x )02 = σ , if ρ < 0, (3.6) N − (N − 2)ρ N + (N − 2)ρ 2 − 2ρ(N − 1)

100

where r 2 =

R.N. Das, S.H. Park / Journal of the Korean Statistical Society 37 (2008) 95–106

Pk

2 i=1 x i .

3.2. Efficient RFORDs with autocorrelated error From (3.4) it is clear that exact hypothetical D-optimal robust first order regression designs with autocorrelated error in observations, may not be possible to arrive at explicitly. Practically it is impossible in this case. Following (3.4), it is clear that one should construct two sets of RFORD or nearly D-ORFOD under the structure (3.1). One set is related to ρ > 0 and the other set is related to ρ < 0. In this subsection we have developed a set of efficient RFORDs under autocorrelated structure (3.1), using augmentation of columns of Hadamard matrices, and the designs thus obtained are not D-optimal but nearly so, having high efficiency around above ninety percent. General methods have not been developed. For the sake of simplicity, here as examples we describe the explicit constructions of nearly D-optimal designs or efficient RFORDs using augmentation of columns of Hadamard matrices H4 , H8 and H16 only. Obviously, similar methods can be applied to other Hadamard matrices to get useful designs, but here other designs are not shown as it would increase the length of the paper. (i) Let h j be the jth column of a Hadamard matrix of order i = 4, 8, 16 where the Hadamard matrices are formed by the Kronecker product method as given below:   1 1 H2 = , H4 = H2 ⊗ H2 , H8 = H2 ⊗ H4 , H16 = H2 ⊗ H8 , etc., (3.7) 1 −1 where ⊗ denotes a Kronecker product. Below we have developed some efficient RFORDs which are very close to D-optimal robust first order designs with design points 11, 15, 19, 27 and 35 from H4 , H8 and H16 only. Other possible designs with different design points have not been considered here as they are not so efficient. From our construction, it is clear that the minimum number of design points is 11 for k = 2 factors. Note that the minimum number of design points of a RFORD with 2 or 3 factors under autocorrelated error is 9 (see Panda and Das (1994)), but those RFORDs are not efficient. These designs satisfy the conditions of robust first order rotatability (i.e., (3.3) with vii = constant 6= µ as in (3) of (3.3)) but are far from satisfying (3.4). 3.2.1. Efficient RFORDs or nearly D-ORFODs for ρ ≥ 0 and ρ ≤ 0 Some efficient RFORDs or nearly D-optimum robust first order regression designs for k = 2 and k = 3, 4 factors with different number of design points for ρ ≥ 0 (ρ ≤ 0) are given in Tables 3.1 and 3.3 (Tables 3.5 and 3.6), respectively. More explicitly the designs D1 and D2 (in Table 3.1) are given in Tables 3.2 and 3.4, respectively. The quantity in () (in 3.1, 3.3, 3.5 and 3.6) is the actual quantity needed for satisfying (3.4). The differences of PN P NTables −1 2 2 and x x always 3 and 1 respectively, from the actual value in (), regardless of the ρ and k as there u=1 iu u=2 iu areP N −1 are 3 zeros in each xi . Also u=1 xiu xi(u+1) is always at least reduced by 4 from (N − 1) in (), due to 3 zeros in each xi but the differences are different (for ρ > 0, ρ < 0 and k) due to the different ordered composition of −1, 0 and +1 in each xi . Note that the designs D1 to D5 (each with 2 factors) which are developed here are efficient RFORDs with design points, respectively, 11, 15, 19, 27 and 35. But one can construct RFORDs with design points 12, 13 and 14 from D1 by adding to it one, two and three central points, respectively. These developed designs (from D1 ) are RFORDs but they lose their efficiency, as the addition of central points makes satisfying (3.4) less likely. Similarly one can construct RFORDs with design points 16, 17 and 18 from D2 by adding to it one, two and three central points, respectively, and so on. The expressions for “efficiency ratios” (ER) have not been shown here for the designs in Tables 3.3 and 3.6 as it would increase the length of the Table. They can be easily derived similarly as for the designs in Table 3.1. The computation of ER of the designs D1 through D9 and the designs D10 through D18 are given in Tables 3.7 and 3.8, and also their respective graphical presentations are give in Figs. 3.1 and 3.2, respectively. For positive and negative values of ρ, we have developed here some of the most efficient RFORDs or nearly D-optimal robust first order designs up to 4 factors. We have also developed such designs for more than 4 factors and also from H32 which are not shown here, because it would increase the length of this paper.

101

R.N. Das, S.H. Park / Journal of the Korean Statistical Society 37 (2008) 95–106 Table 3.1 Nearly D-optimal designs for k = 2 factors when ρ ≥ 0 Design

No. of runs

PN

P N −1 2 u=2 xiu

P N −1

11

8

8

−4

(6= 11)

(6= 9)

(6= −10)

2 u=1 xiu

u=1 xiu xi(u+1)

ER%

D1 : (4)0

(4)0

(4)0

(4)0

x1 = (0, h 2 , 0, h 4 , 0)0 , x2 = (0, h 4 , 0, h 2 , 0)0

2(1−ρ 2 ) ] 11+9ρ 2 +20ρ 2(1−ρ 2 ) 1+ρ ] [ 11−9ρ + 8+8ρ 2 +8ρ

1+ρ + [ 11−9ρ

D2 : (4)0 (8)0 x1 = (0, h 4 , 0, h 2 , 0)0 , (4)0

15

(8)0

x2 = (0, h 2 , 0, h 6 , 0)0

2(1−ρ 2 ) ] 15+13ρ 2 +28ρ 2(1−ρ 2 ) 1+ρ [ 15−13ρ + ] 12+12ρ 2 +16ρ 1+ρ + [ 15−13ρ

12

12

−8

(6= 15)

(6= 13)

(6= −14)

D3 : (8)0 (8)0 x1 = (0, h 2 , 0, h 6 , 0)0 (8)0

19

(8)0

x2 = (0, h 6 , 0, h 2 , 0)0

2(1−ρ 2 ) ] 19+17ρ 2 +36ρ 1+ρ 2(1−ρ 2 ) [ 19−17ρ + ] 16+16ρ 2 +24ρ 1+ρ [ 19−17ρ +

16

16

−12

(6= 19)

(6= 17)

(6= −18)

D4 : (16)0 (8)0 x1 = (0, h 2 , 0, h 10 , 0)0 , (8)0

(16)0

x2 = (0, h 6 , 0, h 2

27

, 0)0

2(1−ρ 2 ) ] 27+25ρ 2 +52ρ 1+ρ 2(1−ρ 2 ) ] [ 27−25ρ + 24+24ρ 2 +40ρ 1+ρ [ 27−25ρ +

24

24

−20

(6= 27)

(6= 25)

(6= −26)

D5 : (16)0 (16)0 x1 = (0, h 2 , 0, h 10 , 0)0 , (16)0

(16)0

x2 = (0, h 10 , 0, h 2

35

, 0)0

2(1−ρ 2 ) ] 35+33ρ 2 +68ρ 2 1+ρ 2(1−ρ ) ] [ 35−33ρ + 32+32ρ 2 +56ρ 1+ρ [ 35−33ρ +

32

32

−28

(6= 35)

(6= 33)

(6= −34)

Table 3.2 Nearly D-optimal design D10 for k = 2 factors and N = 11 Run no.

1

2

3

4

5

6

7

8

9

10

11

x1 x2

0 0

1 1

−1 −1

1 −1

−1 1

0 0

1 1

−1 −1

−1 1

1 −1

0 0

Fig. 3.1. The efficiency of the designs D1 through D9 for ρ ≥ 0.

102

R.N. Das, S.H. Park / Journal of the Korean Statistical Society 37 (2008) 95–106

Table 3.3 Nearly D-optimal designs for k = 3 and 4 factors when ρ ≥ 0 Design

No. of factors

No. of runs

PN

P N −1 2 u=2 xiu

P N −1

D6 : (8)0 (4)0 x1 = (0, h 2 , 0, h 4 , 0)0 , (4)0 (8)0 x2 = (0, h 3 , 0, h 6 , 0)0 ,

3

15

12 (6= 15)

12 (6= 13)

−4 (6= −14)

4

19

16 (6= 19)

16 (6= 17)

−8 (6= −18)

D8 : (16)0 (8)0 x1 = (0, h 2 , 0, h 6 , 0)0 ,

4

27

24

24

−16

(6= 27)

(6= 25)

(6= −26)

D9 : (16)0 (16)0 x1 = (0, h 2 , 0, h 6 , 0)0 ,

4

32

32

−24

(6= 35)

(6= 33)

(6= −34)

(4)0

2 u=1 xiu

(8)0

u=1 xiu xi(u+1)

x3 = (0, h 4 , 0, h 8 , 0)0

D7 : (8)0 (8)0 x1 = (0, h 2 , 0, h 4 , 0)0 , (8)0 (8)0 x2 = (0, h 4 , 0, h 2 , 0)0 , (8)0 (8)0 x3 = (0, h 6 , 0, h 8 , 0)0 , (8)0

(8)0

x4 = (0, h 8 , 0, h 6 , 0)0

(16)0

(8)0

x2 = (0, h 4 , 0, h 2 , 0)0 , (16)0 (8)0 x3 = (0, h 6 , 0, h 14 , 0)0 , (16)0 (8)0 x4 = (0, h 8 , 0, h 10 , 0)0

(16)0

x2 = (0, h 6

(16)0

(16)0

, 0, h 2

(16)0

35

, 0)0 ,

x3 = (0, h 10 , 0, h 14 , 0)0 , (16)0

(16)0

x4 = (0, h 14 , 0, h 10 , 0)0

Table 3.4 Nearly D-optimal design D20 for k = 2 factors and N = 15 Run no.

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

x1 x2

0 0

1 1

−1 −1

−1 1

1 −1

0 0

1 1

−1 −1

1 1

−1 −1

1 −1

−1 1

1 −1

−1 1

0 0

Fig. 3.2. The efficiency of the designs D10 through D18 for ρ ≤ 0.

4. Concluding remarks In this article, we have considered the first order linear regression model with correlated observations. We have considered the most common correlation structure (namely the autocorrelated structure) of the observations or errors.

103

R.N. Das, S.H. Park / Journal of the Korean Statistical Society 37 (2008) 95–106 Table 3.5 Nearly D-optimal designs for k = 2 factors when ρ ≤ 0 Design

No. of runs

PN

P N −1 2 u=2 xiu

P N −1

11

8

8

0

(6= 11)

(6= 9)

(6= 10)

12

12

4

(6= 15)

(6= 13)

(6= 14)

16

16

8

(6= 19)

(6= 17)

(6= 18)

2 u=1 xiu

u=1 xiu xi(u+1)

D10 : (4)0

(4)0

(4)0

(4)0

(4)0

(8)0

(4)0

(8)0

(8)0

(8)0

(8)0

(8)0

x1 = (0, h 3 , 0, h 4 , 0)0 , x2 = (0, h 4 , 0, h 3 , 0)0 D11 : x1 = (0, h 3 , 0, h 7 , 0)0 ,

15

x2 = (0, h 4 , 0, h 5 , 0)0 D12 : x1 = (0, h 5 , 0, h 8 , 0)0 ,

19

x2 = (0, h 7 , 0, h 5 , 0)0 D13 : (16)0 (8)0 x1 = (0, h 5 , 0, h 13 , 0)0 , (8)0

(16)0

x2 = (0, h 7 , 0, h 9

27

, 0)0

24

24

16

(6= 27)

(6= 25)

(6= 26)

D14 : (16)0 (16)0 x1 = (0, h 9 , 0, h 13 , 0)0 , (16)0

(16)0

x2 = (0, h 13 , 0, h 9

, 0)0

35

32

32

24

(6= 35)

(6= 33)

(6= 34)

ER 2(1−ρ 2 ) ] 11+9ρ 2 −20ρ 2(1−ρ 2 ) 1+ρ ] [ 11−9ρ + 8+8ρ 2

1+ρ [ 11−9ρ +

2(1−ρ 2 ) ] 15+13ρ 2 −28ρ 2(1−ρ 2 ) 1+ρ [ 15−13ρ + ] 12+12ρ 2 −8ρ

1+ρ [ 15−13ρ +

2(1−ρ 2 ) ] 19+17ρ 2 −36ρ 1+ρ 2(1−ρ 2 ) [ 19−17ρ + ] 16+16ρ 2 −16ρ 1+ρ [ 19−17ρ +

2(1−ρ 2 ) ] 27+25ρ 2 −52ρ 1+ρ 2(1−ρ 2 ) [ 27−25ρ + ] 24+24ρ 2 −32ρ 1+ρ [ 27−25ρ +

2(1−ρ 2 ) ] 35+33ρ 2 −68ρ 1+ρ 2(1−ρ 2 ) [ 35−33ρ + ] 32+32ρ 2 −48ρ 1+ρ [ 35−33ρ +

Several causes of autocorrelation that may occur in the observations are given in Section 3. Assuming a general correlated error structure we have derived the necessary and sufficient conditions for D-optimality of first order regression designs. We have developed a class of designs which are not exact optimum robust designs but are nearly so. The designs developed here have higher efficiency for positive values of the correlation coefficient ρ (which is mostly seen in practical situations) than the negative values of ρ. Also efficiency decreases with the increase in number of factors for fixed number of design points, and for the same factor the efficiency increases with the increase of number of design points, and for k = 2 factors efficiency for some designs are always greater than ninety percent. The designs developed here are robust rotatable. Exact optimum robust designs have not been developed here. It is very difficult to construct exact optimum robust designs, in fact it is practically impossible. Using this method one can construct designs with efficiency greater than ninety five percent (like D5 , in Table 3.7, for ρ ≥ 0.3) and these designs are very close to D-ORFODs. Finally we have developed a method of analysis of the designs which gives the BLUE of all the regression parameters except the intercept, which is often unimportant. When we are remote from the optimum operating conditions for the system (e.g., quality-improvement, clinical trail experiments, etc.), we usually assume that a first order model is an adequate approximation to the true surface in a small region of the explanatory variables x’s. The designs developed here are more appropriate for searching a new region in which the process or product is improved. Acknowledgements The authors are very much indebted to two referees who contributed valuable suggestions and comments that helped improve this paper very much.

104

R.N. Das, S.H. Park / Journal of the Korean Statistical Society 37 (2008) 95–106

Table 3.6 Nearly D-optimal designs for k = 3 and 4 factors when ρ ≤ 0 Design

No. of factors

No. of runs

PN

P N −1 2 u=2 xiu

P N −1

D15 : (8)0 (4)0 x1 = (0, h 2 , 0, h 7 , 0)0 , (4)0 (8)0 x2 = (0, h 3 , 0, h 4 , 0)0 ,

3

15

12 (6= 15)

12 (6= 13)

0 (6= 14)

4

19

16

16

4

(6= 19)

(6= 17)

(6= 18)

24

24

12

(6= 27)

(6= 25)

(6= 26)

32

32

20

(6= 35)

(6= 33)

(6= 34)

(4)0

2 u=1 xiu

(8)0

u=1 xiu xi(u+1)

x3 = (0, h 4 , 0, h 3 , 0)0

D16 : (8)0 (8)0 x1 = (0, h 3 , 0, h 7 , 0)0 , (8)0

(8)0

(8)0

(8)0

(8)0

(8)0

x2 = (0, h 7 , 0, h 3 , 0)0 ,

x3 = (0, h 4 , 0, h 5 , 0)0 ,

x4 = (0, h 5 , 0, h 4 , 0)0

D17 : (16)0 (8)0 x1 = (0, h 3 , 0, h 13 , 0)0 , (8)0

(16)0

(8)0

(16)0

(8)0

(16)0

x2 = (0, h 4 , 0, h 9

x3 = (0, h 5 , 0, h 7

x4 = (0, h 7 , 0, h 5

4

, 0)0 , , 0)0

D18 : (16)0 (16)0 x1 = (0, h 5 , 0, h 13 , 0)0 , (16)0

(16)0

x2 = (0, h 13 , 0, h 5 (16)0

x3 = (0, h 7

(16)0

x4 = (0, h 9

(16)0

, 0, h 9

(16)0

, 0, h 7

27

, 0)0 ,

4

35

, 0)0 ,

, 0)0 , , 0)0

Table 3.7 Computation of efficiency ratio of the designs from D1 to D9 ρ 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Efficiency of D1

D2

D3

D4

D5

D6

D7

D8

D9

80.000 78.392 78.842 80.805 83.810 87.385 91.056 94.394 97.080 98.951

85.714 84.913 85.612 87.339 89.681 92.260 94.736 96.850 98.447 99.483

88.889 88.415 89.113 90.577 92.562 94.460 96.317 97.852 98.976 99.677

92.308 92.091 92.684 93.780 95.120 96.490 97.722 98.710 99.409 99.826

94.118 93.996 94.493 95.361 96.398 97.440 98.360 99.085 99.591 99.885

84.21 79.29 77.37 77.76 79.97 83.50 87.78 92.13 92.13 98.58

86.95 82.87 81.26 81.54 83.32 86.23 89.81 93.49 96.67 98.90

90.91 88.22 87.32 87.75 89.21 91.35 93.82 96.21 98.16 99.44

93.02 91.02 90.42 90.84 92.02 93.72 95.59 97.35 98.75 99.64

Table 3.8 Computation of efficiency ratio of the designs from D10 to D18 ρ 0.0 −0.1 −0.2 −0.3 −0.4 −0.5 −0.6 −0.7 −0.8 −0.9

Efficiency of D10

D11

D12

D13

D14

D15

D16

D17

D18

80.00 70.75 63.93 58.99 55.50 53.10 51.52 50.55 50.05 49.90

85.71 79.46 75.03 71.87 69.71 68.26 67.35 66.82 66.58 66.54

88.89 84.17 80.87 78.59 77.04 76.02 75.38 75.03 74.88 74.78

92.31 89.14 86.98 85.51 84.53 83.90 83.51 83.30 83.23 83.24

94.12 91.74 90.13 89.05 88.33 87.87 87.60 87.45 87.41 87.43

84.21 73.55 65.82 60.28 56.38 53.69 51.91 50.80 50.19 49.96

86.96 78.49 72.52 68.34 65.45 63.49 62.22 61.45 61.04 60.91

90.91 85.14 81.15 78.40 76.52 75.25 74.47 73.99 73.76 73.70

93.02 88.65 85.66 83.61 82.22 81.30 80.72 80.38 80.21 80.18

R.N. Das, S.H. Park / Journal of the Korean Statistical Society 37 (2008) 95–106

105

Appendix. Estimation of regression parameters and error variance From the nature of experimentation, we can guess the pattern of the correlation structure of errors in observations but the correlation parameter or parameters are not known. The situations described in Section 3 may give rise to an autocorrelation structure of observations as given in (3.1). As the form of the correlation structure is known, we can generate observations using the designs developed in Section 3.2. Below we have developed a method for estimating regression parameters in the case of a first order model as in (2.1) under the autocorrelation structure as in (3.1) and error variance σ 2 . The model, as explained in (2.1), is yu = β0 +

k X

βi xiu + eu ;

1 ≤ u ≤ N,

or,

i=1

Y = X β + e,

∼

∼

(1)

∼

where Y is the vector of recorded observations on the study variable y, β = (β0 , β1 , β2 , . . . , βk , )0 , is the vector of ∼

∼

regression coefficients of order k + 1 and X is the design matrix as developed in Section 3.2. W0 = {D(e ) = [σ 2 {ρ |i− j| } 1≤i, j≤N ] = W N ×N (ρ)}, ∼

is known as the autocorrelation error structure class and σ 2 may be known or unknown. Therefore, the generated observations on the study variable y by a well planned design matrix X , as developed in Section 3.2, are y1 , y2 , . . . , y N . For N observations of y’s, let us define z u = yu − ρyu−1 ; 2 ≤ u ≤ N , or, z u = β0 (1 − ρ) +

k X

(2)

βi (xiu − ρxi(u−1) ) + (eu − ρe(u−1) ),

(3)

i=1

or, z u = α0 +

k X

βi siu + u ,

say, 2 ≤ u ≤ N ,

(4)

i=1

where α0 = β0 (1 − ρ), siu = (xiu − ρxi(u−1) ); 1 ≤ i ≤ k and u = (eu − ρe(u−1) ), 2 ≤ u ≤ N , or,

Z = S η + , ∼

∼

(5)

∼

where Z = (z 2 , z 3 , . . . ., z N )0 , η = (α0 , β1 , β2 , . . . , βk )0 , S = ( 1 : (siu ); 1 ≤ i ≤ k; 2 ≤ u ≤ N ), ∼

∼

∼

E(u ) = 0, 2 ≤ u ≤ N ; Var (u ) = σ 2 (1 − ρ 2 ), 2 ≤ u ≤ N ; Cov. (u , u 0 ) = 0, 2 ≤ u 6= u 0 ≤ N ; and E( ) = 0 , D( ) = σ 2 (1 − ρ 2 )I N −1 . ∼

∼

∼

Here z u ’s, siu ’s and ρ are unknown but yu ’s and xiu ’s are known. The scheme for the calculations which can be followed is: (i) Assume some value of ρ ∈ (−1, 1) and then calculate siu = (xiu − ρxi(u−1) ); 1 ≤ i ≤ k, and z u = yu − ρyu−1 ; 2 ≤ u ≤ N , using the assumed value of ρ. (ii) As S-matrix and Z in (4) are known from (i), using the ordinary least squares (OLS) method, the best linear ∼ unbiased estimate (BLUE) of η is ∼

ηˆ = (S 0 S)−1 (S 0 Z ).

∼

(6)

∼

αˆ 0 Then αˆ 0 = βˆ0 (1 − ρ), where ρ is the assumed value of ρ as in (i) and βˆ0 = (1−ρ) . PN Pk (iii) Using the estimate of η as ηˆ as in (ii), calculate S(ηˆ , ρ) = u=2 (z u − αˆ0 − i=1 βˆi siu )2 . ∼

∼

∼

The same sequence of calculations (i.e., (i)–(iii)) may be followed for different permissible values of ρ ∈ (−1, 1), by using a computer program. Final estimates are obtained by searching through the results of the computer program for a value of ρ such that a combination of ρ and η is found that minimizes S(ηˆ , ρ). ∼

∼

106

R.N. Das, S.H. Park / Journal of the Korean Statistical Society 37 (2008) 95–106

An unbiased estimate of [σ 2 (1 − ρ 2 )] = σ12 is ( Z 0 Z ) − ηˆ (S 0 S)ηˆ 0

σˆ 12

=

∼ ∼

∼

∼

(N − 1) − (k + 1)

.

(7)

Therefore, an estimate of σ 2 is σˆ2 =

0 ( Z 0 Z ) − ηˆ (S 0 S)ηˆ 1 ∼ ∼ ∼ ∼ , . (1 − ρˆ 2 ) (N − 1) − (k + 1)

(8)

where ρˆ and ηˆ are the final estimates of ρ and η and Z is computed on the basis of the final estimate of ρ. It is clear ∼ ∼ ∼ from (5) that the estimates of all regression parameters except β0 are the best linear unbiased estimates. References Anderson, T. W. (1984). Introduction to multivariate statistical analysis (2nd ed.). New York: Wiley and Sons. Bischoff, W. (1992). On exact D-optimal designs for regression models with correlated observations. Annals of the Institute of Statistical Mathematics, 44, 229–238. Bischoff, W. (1995). Determinant formulas with applications to designing when the observations are correlated. Annals of the Institute of Statistical Mathematics, 47, 385–399. Bischoff, W. (1996). On maximin designs for correlated observations. Statistics and Probability Letters, 26, 357–363. Box, G. E. P., & Draper, R. N. (1987). Emperical model building and response surfaces. New York: Wiley and Sons. Box, G. E. P., & Hunter, J. S. (1957). Multifactor experimental designs for exploring response surfaces. Annals of Mathematical Statistics, 28, 195–241. Das, R. N. (1997). Robust second order rotatable designs: Part-I. Calcutta Statistical Association Bulletin, 47, 199–214. Das, R. N. (2003a). Slope rotatability with correlated errors. Calcutta Statistical Association Bulletin, 54, 57–70. Das, R. N. (2003b). Robust second order rotatable designs: Part-III. Journal of the Indian Society of Agricultural Statistics, 56, 117–130. Das, R. N. (2004). Construction and analysis of robust second order rotatable designs. Journal of Statistical theory and Applications, 3, 325–343. Das, R. N., & Park, S. H. (2006). Slope rotatability over all directions with correlated errors. Applied Stochastic Models in Bussiness and Industry, 22, 445–457. Gennings, C., Chinchilli, V. M., & Carter, W. H., Jr. (1989). Response surface analysis with correlated data: A non-linear model approach. Journal of the American Statistical Association, 84, 805–809. Hader, R. J., & Park, S. H. (1978). Slope-rotatable central composite designs. Technometrics, 20, 413–417. Kiefer, J., & Wynn, H. P. (1981). Optimum balanced block and latin square designs for correlated observations. The Annals of Statistics, 9, 737–757. Kiefer, J., & Wynn, H. P. (1984). Optimum and minimax exact treatment designs for one-dimensional autoregressive error processes. The Annals of Statistics, 12, 431–449. Khuri, A. I., & Cornell, J. A. (1996). Response surfaces. New York: Marcel Dekker. Myers, R. H., Khuri, A. I., & Carter, W. H., Jr. (1989). Response surface methodology: 1966–1988. Technometrics, 31, 137–157. Myers, R. H., & Montgomery, D. C. (2002). Response surface methodology. New York: Wiley and Sons. Park, S. H. (1987). A class of multifactor designs for estimating the slope of response surfaces. Technometrics, 29, 449–453. Panda, R. N., & Das, R. N. (1994). First order rotatable designs with correlated errors. Calcutta Statistical Association Bulletin, 44, 83–101. Pukelsheim, F. (1993). Optimal experimental design. New York: Wiley and Sons.