KATHOLIEKE UNIVERSITEIT LEUVEN DEPARTEMENT TOEGEPASTE ECONOMISCHE WETENSCHAPPEN RESEARCH REPORT 0126 AN ADJUSTMENT ALGORITHM FOR THE CONSTRUCTION OF OPTIMAL RUN ORDERS by l. TACK M. VANDEBROEK D/2001/2376/26 An adjustment algorithm for the construction of optim::ll run orners L. Tack M. Vandebroek Kaiholieke Universiieii Leuven, Belgium Using a systematic run order can be the proper way to conduct an experiment when a temporal trend is present. The construction of run orders that are optimally balanced for time trend effects is based on the maximization of the information on the important parameters whereas the parameters of the postulated time trend are treated as nuisance parameters. In this paper, an adjustment algorithm is presented to improve the efficiency of the run orders obtained from a search over a predefined set of candidate points. This is done by repeatedly moving the design points or the time points of the candidate list a small amount along their axes as long as an improvement in the efficiency is obtained. It is illustrated that the adjustment algorithm involves substantial increases in the efficiency of the run orders. The use of the adjustment algorithm in addition to a search over a coarse grid of candidate points is especially recommended in situations where the computation time has to be kept within reasonable limits. Keywords: 'Droptimality; adjustment algorithm; run order; time trend 1 1 Introduction Performing a series of measurements in a time sequence possibly creates time order dependence in the results. For example, when a batch of material is created at the beginning of an experiment and the treatments are to be applied to the experimental units formed from the material over time, there could be an unknown effect due to aging of the material. Unwanted time trend effects can also be caused by poisoning of a catalyst, steady build-up of deposits in a test engine, equipment wear-out, analyst fatigue, warm-up in laboratories, etc. For instance, Freeny and Lai (1997) mention an experiment taken from the electronics industry in which a photolithographic polisher shows the tendency to drift lower through time. Another example comes from Joiner and Campbell (1976) who describe an experiment in which the measurements drift with time due to the build-up of carbon in a spectrophotometer. The next section shortly reviews the literature on the existence and the construction of run orders that are optimally balanced for time trend effects. Section 3 considers the problem in light of optimal design theory. Section 4 presents the adjustment algorithm and in section 5 some examples are given to illustrate the use of the algorithm. 2 Trend-free run orders When the experimenter has knowledge about the nature of the time trend it is recommended to construct a run order in which the estimates of the treatment differences or the factorial effects are little disturbed by the presence of the time trend. This time dependence is usually approximated by a polynomial function of order q. Let g(t) be the q x 1 vector representing the polynomial expansion for the time trend, expressed as a function of time t E [-1,1]. With f3 the q x 1 vector of parameters of the polynomial time trend, f(x) the p x 1 vector representing the polynomial expansion of x = [Xl,'" ,xf]' for the response and a the p x 1 vector of important parameters, let the statistical model be y = f'(x)a + g'(t)(3 + E. (1) The error terms are assumed to be independent and identically distributed with mean zero and constant variance Contrary to simply observing the system twice for successive replicate points, it is understood that the transition from one run to the next one involves some intervention. As a result, two observations on the same design point are clearly subject to different errors from two runs and the assumption of uncorrelated error terms is justified. The fact that there is no interaction between the controllable variables Xi, with i E {I, ... ,j}, and time t usually holds in practice. For n observations, it is convenient to rewrite (1) as (J;. y = Fa + Gf3 + c, 2 (2) where F and G represent the n x p and the n x q design matrices respectively. The least-squares estimators Ii are given by (3) where In is the n-dimensional unity matrix. The factorial effects are now said to be q-trend-free if the least-squares estimator Ii is the same as when the time trend of qth order is not present. It is easy to verify that the least-squares estimator (3) is equal to the estimator Ii = (F/F)-lF/y in the absence of trend effects if and only if F'G = O. Consequently, a run order is called q-trend-free if the columns of F are orthogonal to the columns of G. Cox (1951) was the first to study the construction of run orders for the estimation of treatment effects in the presence of a polynomial time trend. He presents trend-free run orders for a number of simple design problems for categorical variables. The most important contributions to trend-resistant experimental design for quantitative variables come from Cheng (1985) and Daniel and Wilcoxon (1966). Cheng (1985) formulates the Generalized Foldover Scheme with which generator sets can be created to construct trendfree run orders of full or fractional factorial designs. Based on the trend-free properties of the standard run order of factorial designs, Daniel and Wilcoxon (1966) develop another method for the construction of trend-free run orders but every run order derived with their method can also be obtained by applying the Generalized Foldover Scheme. A more extensive survey of the literature on trend-resistant design of experiments is given in Tack and Vandebroek (2001). 3 Vt-optimal run orders Although the above references on the existence and the construction of trend-free run orders are of great use to the practitioner, they all have a number of important shortcomings. Firstly, none of the references deals with the numerous situations in which the required orthogonality cannot be attained. Secondly, there is no consideration of the connection between trend-resistance and the generalized variance or the precision of the parameter estimates. Thirdly, all approaches assume that the treatments or factor level combinations have already been chosen. Finally, the approaches rely on assumptions that simplify reality too much. For instance, it is assumed that the number of time points is equal to the number of observations and that the time points are equally spaced. Besides, only 2- and 3-level factorial designs are considered, interactions of higher order are assumed to be negligible and constrained experimental regions cannot be dealt with. To fill these gaps, Atkinson and Donev (1996) present a generic exchange algorithm to compute run orders that are optimally balanced for time trends. Tack and Vandebroek (2001) extend the exchange algorithm to allow for cost considerations and they develop a methodology to construct run orders that maximize the experimental efficiency defined as the amount of information obtained per unit cost, rather than focusing too much on 3 a purely statistical efficiency. Contrary to offering a catalogue of specific solutions to design problems with a special structure, they provide the experimenter with a broadly applicable method to tackle a wide range of practical design problems. When primary interest is in the precision of the parameter estimates, the algorithm computes optimal run orders by maximizing the information on the important parameters a whereas the parameters f3 of the time trend are treated as nuisance parameters. The resulting run order is referred to as the Vcoptimal run order 8vt and the optimality criterion equals F'F V t = F'G I I G'F G'G IG'GI = IF'F - F'G(G'G)-lG'FI . (4) In the absence of time trend effects, the V-optimal design 8v maximizes the determinant of the information matrix, i.e. V = IF'FI. The V-optimal design and the Vroptimal run order are then compared to each other through the trend factor TF(8 ) v, = {Vt(8v ,) } lip V(8v ) (5) The power lip assures that the trend factor is independent of the dimension of the model. This means for instance that a run order 8v , with trend factor 0.5 has to be replicated twice in order to be equally informative as the V-optimal design 8v . It follows readily that a trend-free run order, i.e. a run order with F'G = 0, attains the maximum value TF(8v ,) = 1. In situations where it is impossible to obtain completely trend-free run orders, TF(8vt ) will be less than 1. 4 An adjustment algorithm The exchange algorithm of Tack and Vandebroek (2001) computes 'Dcoptimal run orders by the allocation of n observations selected from a candidate list of d distinct design points to n out of h available time points in such a way as to maximize optimality criterion (4). The number n and the list of time points are user-specified. The candidate set of design points is also user-specified or can be computed as follows: the number of equally spaced levels per factor equals two, three, four or more depending on the fact whether the polynomial expansion f(x) is of first, second, third or higher order respectively. For instance, for a second-order polynomial the three factor levels under study are -1, 0 and 1. These three factor levels are commonly encountered in response surface methodology. This is due to two main reasons. Firstly, Farell et al. (1967) demonstrate that the factor levels of the optimal continuous design for a second-order polynomial over a hypercubic experimental region are -1, 0 and 1. Secondly, many practitioners prefer to use only a small number of distinct factor levels in order to keep the experimental effort within reasonable limits. 4 However, optimal design theory has shown that exact optimal designs found by searching over a continuous experimental region often contain other factor levels than -1, 0 and 1. For instance, in the absence of trend effects, Box and Draper (1971) and Donev and Atkinson (1988) show that for some factors other levels than -1, 0 and 1 lead to more efficient exact designs. Box and Draper (1971) analytically compute exact V-optimal designs for second-order response models whereas Donev and Atkinson (1988) present an adjustment algorithm to improve the efficiency of V-optimal designs with or without fixed block effects. This is done by moving the design points along their factor axes as long as an increase in the V-efficiency is achieved. For a second-order polynomial model with uncorrelated observations and the number of factors up to five, Donev and Atkinson (1988) and Atkinson and Donev (1992) show that the effect of their adjustment algorithm is largest when the number of observations is equal to or just larger than the number of parameters p. For instance, for a second-order polynomial in three factors and n = p = 10, they obtain improvements in the V-efficiency of about 3%. For a second-order model with fixed block effects and one or two factors, they show that substantial gains can be realized when there are only a few observations per block. Goos (2001) investigates the gains obtained by using an adjustment algorithm to improve the V-efficiency of a number of different designs with fixed or random block effects. For a full second-order response model, improvements of up to 10% are obtained. On the whole, the efficiency gains of his adjustment algorithm are superior to those obtained by searching over a fine hypercubic grid with 11 equally spaced levels for each factor. The best efficiencies are achieved by using both a fine grid and the adjustment algorithm. Generally speaking, the combined approach gives the best results for small n, small block sizes and highly correlated observations. The previous references have clearly illustrated the large efficiency gains that result from using an adjustment algorithm. For that reason, the next section will present an adjustment algorithm to improve the trend factors of the VI-optimal run orders computed with the exchange algorithm of Tack and Vandebroek (2001). 4.1 Description of the algorithm It is important to point out that the adjustment algorithm to be outlined in this section is a generalization of the adjustment algorithms of Donev and Atkinson (1988) and Goos (2001) because the time points can be adjusted too. The adjustment algorithm computes the effect on the Vcoptimality criterion of moving any design point a small amounthenceforth referred to as the step length-along its factor axes and of moving any time point a small step along the time axis. As a matter of fact, interest is only in design changes that lead to design points still lying within the experimental region and time points still belonging to [-1,1]. This means that at most 2nf modifications for the design points and 2n changes for the time points have to be investigated. The design modification that leads to the largest increase in the criterion value is carried out and the process is repeated until no further progress can be made. The step length is then halved and the 5 improvement process starts again. The algorithm attains its final iteration when the step lengths for the design points and the time points become smaller than their user-specified Illinimum value. The input to the algorithm consists of the number of factors f, the order and the number of parameters p of the response model, the polynomial expansion for the response model f (x), the order and the number of parameters q of the time trend, the polynomial expansion for the time trend g(t) and the number of observations n. The initial step length for the design points and the initial step length for the time points, henceforth denoted as 51 and 52 respectively, have to be supplied too. The minimum step length for the design points and the minimum step length for the time points, referred to as 51,min and 52,min respectively, are also user-specified. Another input parameter specifies the design changes to be evaluated: only design points considered to be moved, only time points to be moved or both design points and time points allowed to be changed. The first option is chosen when the time points are fixed by the design problem at hand whereas the second option refers to the situation in which the different factor level combinations have already been chosen. The experimenter has the additional possibility to impose a minimum distancein the sequel of this paper denoted as MD-to be maintained between two successive time points. Such a situation occurs when a fixed measurement time has to be taken into account. Finally, the input also contains the Vcoptimal run order to be used as the starting run order in the adjustment algorithm. The appendix gives a formal outline of the adjustment algorithm. 4.2 Update formulae A considerable reduction in the computation time is achieved when powerful update formulae are used to evaluate the effect of moving a time point or a co-ordinate of a design point along its axis. This is thanks to the fact that the matrices in the numerator and the denominator of optimality criterion (4) can be written as a sum of outer products H'H = L h(Xi' tk)h'(Xi' tk), (6) 'o'(x;.tkl where h(Xi, tk) = [f'(Xi) g'(tk)l' or h(Xi' tk) = g(tk) for the numerator or the denominator of (4) respectively. Let the current step length for the design points and the time points be written as 51 and 52 respectively. After moving the rth factor of design point Xj from level Xj,. to level Xj,. ± 51 and time point tl to t/ ± 52, (6) becomes where the f-dimensional column vector e,. contains element 1 at position T and zeros elsewhere. According to theorem 18.1.1 of Harville (1997), the determinant of (7) can be 6 written as !H'H![{l- h'(xj, tz)(H'H)-lh(Xj, tl)} x {I + h' (Xj ± 81 e r , tl ± 82)(H'H)-1 h(Xj ± 81 en tl ± 82)} +{h'(Xj ± 81er, tz ± 82)(H'H)-lh(Xj, tl)}2]. (8) Based on (8), the effect on the criterion value (4) of adjusting a design point or a time point can readily be calculated without the need for computationally intensive determinant operations. 5 Illustrative examples This section illustrates the utility of the adjustment algorithm described in section 4. The first example concerns quadratic regression in one explanatory variable. In the second example, adjusted run sequences for the 24 full factorial design are constructed. The full second-order polynomial in three factors is investigated in the third example. The final example traces whether the adjustment algorithm serves as a valid alternative to the exchange algorithm of Tack and Vandebroek (2001) when the total computation time has to be kept low. 5.1 Quadratic regression III one variable Consider the problem of designing an experiment for the estimation of a quadratic model f(x) = [1, X, x 2 ]' and coded factor levels -1, 0 and 1. Use is made ofthe exchange algorithm of Tack and Vandebroek (2001) to compute Vt-optimal run orders for postulated time trends of first, second, third and fourth order, denoted as gl(t), g2(t), g3(t) and g4(t) respectively. The design sizes equal n = 7, n = 8, n = 9 or n = 10 and there are as much equally spaced time points between -1 and 1. The trend factors of the computed Vt-optimal run orders bVt are given in table 1 in the columns with label TF(5v ,). For instance, for n = 7 and n = 9 the Vroptimal run orders are completely trend-free for a linear time trend gl(t). However, the Vroptimal run orders are poorly balanced for higher-order time trends. It will now be investigated whether an adjustment of the time points leads to an increase in the trend factors of the VI-optimal run orders or not. Therefore, the adjustment algorithm is used with the initial step length for the time points S2 = 2, the minimum step length S2,min = 10- 5 and the minimum distance MD = 10- 5 . The trend factors of the adjusted run orders badj are also given in table 1. It can easily be seen that considerable increases in the trend factors are achieved. For instance, adjusting the time points of the Vt-optimal run order for 17, = 7 and time trend g2(t) augments the trend factor from 0.712 to 0.752. This is an improvement of about 5.6%. The increase in the trend factor is largest for the run orders computed with time trend g4(t). In that case improvements of up to 30% are observed. 7 Table 1: Trend factors of the 'Droptimal run orders and the adjusted run orders for quadratic regression in one variable and different numbers of observations n=7 gl (t) g2(t) g3(t) g4(t) TF(ov,) 1.000 0.712 0.677 0.451 n=8 TF(oadj) 1.000 0.752 0.689 0.591 TF(ov,) 0.999 0.743 0.706 0.545 TF(oadj) 1.000 0.817 0.763 0.688 n=9 TF(ov,) TF(oadj) 1.000 1.000 0.818 0.753 0.763 0.705 0.559 0.689 n = 10 TF(ov,) 0.999 0.754 0.731 0.579 TF(oadj) 1.000 0.846 0.793 0.725 As an illustration, the time points of the adjusted run orders for n = 7, n = 8, n = 9 and n = 10 are given in figure 1 to figure 4 respectively. The adjusted time points are indicated by means of small vertical dashes. For a linear time trend gl(t) and any number of observations n, the time points are approximately uniformly spread over the entire time range [-1,1]. This observation closely resembles the usual assumption of equally spaced time points between -1 and 1. However, the figures reveal that this assumption is no longer optimal for higher-order time trends. For instance, it can clearly be seen that for time trend g4(t) the time points are more concentrated around t = 0 and the boundaries t = -1 and t = 1. (a) gl(t) -1 (b) g2(t) I I -1 0 -1 0 (c) g3(t) (d) g4(t) ! II ! I -1 I! 0 Figure 1: Time points of the adjusted run orders for quadratic regression in one variable and n = 7 8 (a) gl(t) -1 0 (b) g2(t) II II -1 0 -1 0 (c) g3(t) (d) g4(t) ! I! !I -1 0 Figure 2: Time points of the adjusted run orders for quadratic regression in one variable and n = 8 (a) gl(t) ! -1 0 (b) g2(t) II -1 0 (c) g3(t) (d) g4(i) II II -1 0 !I I! I -1 0 I! Figure 3: Time points of the adjusted run orders for quadratic regression in one variable and n = 9 9 (a) gl(t) -1 0 -1 0 -1 0 -1 0 (b) g2(t) II (c) g3(t) (d) g4(t) Figure 4: Time points of the adjusted run orders for quadratic regression in one variable and n = 10 It is interesting to note that in almost every case an adjustment of the design points -1, o and 1 of the 1\-optimal run orders also has a positive effect on the criterion value. For instance, whereas the V-optimal design for a quadratic regression in one variable and n = 9 observations has three observations at each of the levels -1, 0 and 1, the corresponding adjusted run order for time trend g4(t) possesses the levels -1, -0.297, 0.016, 0.063 and 1. As compared to the Vcoptimal run order for the candidate set -1, 0 and 1, the improvement in the trend factor amounts to 0.7%. Although the improvements are very small, the important conclusion is that trend-robust run orders do not necessarily have the same support as their V-optimal counterparts. 5.2 The 24 factorial design The aim of this example is to verify whether the benefit of the adjustment algorithm can be generalized to any initial step length 52 and minimum distance MD. This example starts with the computation of Vcoptimal run orders for the 24 full factorial design with polynomial expansion 24 = 16 observations, as much equally spaced time points between -1 and 1 and postulated time trends gl (t), g2 (t), g3 (t) and g4 (t). The trend factors of the Vcoptimal run orders are 1, 0.900, 0.849 and 0.758 respectively. It follows that for time trends of order two or higher complete trend-resistance cannot be attained. The adjustment algorithm is again used to trace whether improvements in the trend factors of the Vt-optimal run orders for time trends g2(t), g3(t) and g4(i) can be obtained. When only the design points of the Vt-optimal run orders are subject to adjustment, no improvements in the trend factors could be achieved. This result does not come as a 10 surprise since the 24 full factorial design is a 'V-optimal design. An adjustment of the time points is also investigated and the results are shown in table 2. Table 2: Trend factors of the 'Vcoptimal run orders and the adjusted run orders for the 24 full factorial design and initial step length S2 = 0.10 g2(t) g3(t) g4(t) TF( c5v ,) 0.900 0.849 0.758 0.00 0.903 0.871 0.808 0.02 0.903 0.868 0.803 TF(c5adj ); MD = 0.04 0.06 0.903 0.903 0.865 0.863 0.794 0.790 0.08 0.903 0.861 0.786 0.10 0.902 0.858 0.778 Each row refers to the adjusted run orders for a particular time trend and displays the trend factors for several minimum distances MD and initial step lengths S2 = 0.10. It is important to note that the results are independent of the initial step length S2' The minimum step length for the time points was set equal to S2,min = 10- 5 . The table shows that for all cases investigated an adjustment of the time points leads to an improvement in the trend factor. Consequently, it again follows that the usual assumption of equally spaced time points is not always the best one. For instance, for a second-order time trend g2(t) and minimum distance MD equal to 0.00, the upper left cell shows that the trend factor of the adjusted run order is 0.903. As compared to the trend factor 0.900 of the corresponding 'Vcoptimal run order, this is only a small improvement. The increases in the trend factor are more pronounced as the order of the postulated time trend grows larger. For instance, for time trend g4(t) and minimum distance MD = 0.00, improvements of more than 6.5% are achieved. As a matter of fact, the smaller the minimum distance to be maintained between two successive time points, the larger the improvement in the trend factor. As an illustration, for the postulated time trends g2(t), g3(t) and g4(t), the time points of the adjusted run orders are displayed in figure 5, figure 6 and figure 7 respectively. The adjusted time points are again indicated by means of small vertical dashes. As a matter of fact, the smaller the minimum distance to be maintained between successive time points, the less the adjusted time points are equally spaced. 5.3 The 3 3 factorial design In this section the effect of adjusting the design points will be investigated. Attention is restricted to experiments in which three factors are presumed to influence the response of interest, the polynomial expansion 11 1_ \ \ eL) 10. Jfl\ lV1V = n nn Ii i V.VV i i -1 (b) MD = 0.02 = 0 H II II I II I II I II II -1 (C) MD = 0.04 = H II -1 I I I I I I I I 0 (d) MD = 0.06 = H I I I I I -1 (e) MD = 0.08 = H I I I I I I -1 (f) MD = 0.10 I I I I 0 = f----l -1 Figure 5: Time points of the adjusted run orders for the 24 full factorial design, time trend g2(t) and initial step length 0.10 (a) MD = 0.00 -1 (b) MD = 0.02 = H III IIII IIII -1 (c) MD = 0.04 = H III III I -1 I I I I I I I -1 = H = 0.10 = III I I I I I I I 0 I I I I -1 (f) MD IIII 0 (d) MD = 0.06 =H (e) MD = 0.08 III 0 I I I I I I I I I I 0 f----l -1 Figure 6: Time points of the adjusted run orders for the 24 full factorial design, time trend g3(t) and initial step length 0.10 12 (a) MD = 0.00 II II II I -1 (b) MD = 0.02 = H ! 0 II I I I I II II -1 (c) MD = 0.04 = H II I I I (d) MD = 0.06 = H 0 I I I I I I I -1 (e) MD = 0.08 = (f) MD = 0.10 = I I I I I 0 I H II I I I I, I -1 I I I I -1 0 -1 0 I I I I I I I-l Figure 7: Time points of the adjusted run orders for the 24 full factorial design, time trend g4(t) and initial step length 0.10 33 = 27 observations and h = 27 equally spaced time points between -1 and l. The postulated time trends are again written as gl(t), g2(t), g3(t) and g4(t). The trend factors of the Vt-optimal run orders for the full 33 factorial design equal 0.9413, 0.8677, 0.8663 and 0.8230 for the respective time trends. In the sequel of this example, three different approaches will be investigated to improve the trend-resistance of the Vt-optimal run orders of the 33 full factorial design. Firstly, the adjustment algorithm is applied to the Vcoptimal run orders computed with d = 27 and h = 27. The initial step lengths for the design points and the time points are set equal to SI = 0.5 and S2 = 0.05 respectively. Since the design points of the 33 factorial design do not constitute a V-optimal design, improvements in the trend factor can be expected if the design points are considered to be adjusted. The results are given in table 3. Each panel refers to the run orders for a particular time trend. The last three lines of each panel are related to the adjusted run orders. The table shows that for a linear time trend gl(t) the adjustment algorithm does not lead to an improvement in the trend factor of the Vt-optimal run order with d = 27 and h = 27. For higher-order time trends, the adjustment algorithm yields an improvement. For instance, adjusting the time points of the Vt-optimal run order for time trend g2(t), d = 27 and h = 27 gives a run order with trend factor 0.9004. The best results are obtained when both the time points and the design points are considered to be adjusted. In that case, improvements of up to 7% are achieved. This once more demonstrates that the Vcoptimal run orders are in fact only optimal for the specified set of design points and time points. 13 Secondly, the use of a finer grid for the design points or the time points in the exchange algorithm is also studied. Three fine grids are investigated: one with d = 27 design points and h = 41 equally spaced time points between -1 and 1, one with d = 125 design points on a regular 5 x 5 x 5 grid and h = 27 time points and one with d = 125 design points and h = 41 equally spaced time points. The trend factors of the corresponding 'Dcoptimal run orders are given by the last three numbers in the fourth row of each panel. The largest improvement occurs for time trend g4(t), d = 125 and h = 41 and is slightly more than 5%. Finally, the best results are obtained when a search over a fine grid is combined with the adjustment algorithm. This is illustrated in the last three rows and columns of each panel. The most striking improvement is when both the time points and the design points are considered for adjustment. As compared to the trend factors of the 'Dcoptimal run orders computed over the coarse grid with d = 27 and h = 27, improvements of up to 13% are observed. 5.4 Computational aspects of the adjustment algorithm In this example focus is again on design problems for the second-order response function in three factors time trends gl(t), g2(t), g3(t) and g4(t), n = 27 observations and h = 27 equally spaced time points between -1 and 1. The design points are taken from the regular 3 x 3 x 3 grid or from the regular 5 x 5 x 5 grid. The number of different design points then equals d = 33 = 27 or d = 53 = 125 respectively. The exchange algorithm of Tack and Vandebroek (2001) is used to compute 'Dcoptimal run orders OVt for both grids and the results are given in table 4. The first panel of the table shows the trend factors of the computed run orders, whereas the second panel displays the average computation times per try. The best protection against time trend effects is obtained when a fine grid for the candidate points is used (i.e. d = 125 instead of d = 27). For instance, the 'Dcoptimal run order for a second-order time trend and d = 27 has trend factor 0.9217, whereas the trend factor of the corresponding run order for d = 125 has trend factor 0.9637. However, this increase in the protection against temporal dependence goes at the expense of the computation time; the computation of the 'Dt-optimal run order with d = 27 takes on the average 6.2 seconds, whereas the computation time of the 'Dcoptimal run order for d = 125 is about 100 times larger. It follows that for a fine grid of the design points (Le. for large values of d) the exchange algorithm becomes very slow. Consequently, when the total computation time is an important issue to be considered, there is a strong need for a worthy alternative to compute trend-robust run orders within a reasonable computation time. It is now interesting to investigate whether the adjustment algorithm can be used for that purpose or not. The 'Dt-optimal run orders for d = 27 are therefore used as the starting run orders in the adjustment algorithm. Only the design points can be adjusted and their initial step and minimum step are specified large enough to assure that the 14 Table 3: Trend factors of the Vroptimal run orders and the adjusted run orders for the full second-order response model in three factors and n = 27 and for different numbers of design points and time points gl (t) d h Vroptimal run order adjustment of time points adjustment of design points adjustment of time points and design points g2(t) d h Vt-optimal run order adjustment of time points adjustment of design points adjustment of time points and design points g3(t) d h Vroptimal run order adjustment of time points adjustment of design points adjustment of time points and design points g4(t) d h Vt-optimal run order adjustment of time points adjustment of design points adjustment of time points and design points 15 27 41 0.9413 0.9413 0.9413 0.9413 125 27 0.9519 0.9519 0.9964 0.9969 125 41 0.9519 0.9519 0.9956 0.9967 27 27 41 27 0.8677 0.8992 0.9004 0.9109 0.8677 0.8993 0.9275 0.9370 125 27 0.8773 0.9121 0.9139 0.9600 125 41 0.9097 0.9213 0.9461 0.9580 27 27 0.8663 0.8922 0.8667 0.9172 27 41 0.8895 0.9005 0.8899 0.9012 125 27 0.8765 0.9116 0.9053 0.9390 125 41 0.9005 0.9139 0.9402 0.9562 27 27 0.8230 0.8756 0.8238 0.8776 27 41 0.8552 0.8797 0.8561 0.8810 125 27 0.8317 0.8855 0.8575 0.9159 125 41 0.8652 0.8924 0.8919 0.9283 27 27 0.9413 0.9413 0.9413 0.9413 support points of the adjusted run orders form part of the regular 5 x 5 x 5 grid. This assures a valid comparison between the adjusted run orders and the Vroptimal run orders for d = 125. The trend factors of the adjusted run order::; Dadj are shown in the last column of the first panel. For instance, the adjusted run order for g2(t) has trend factor 0.9536, which is very close to trend factor 0.9637 of the Vroptimal run order computed with d = 125. Since the computation time for the adjustment algorithm is negligible to that of the exchange algorithm for large numbers of tries, the computation time of the adjusted run order is only a negligible amount larger than 6.2 seconds. This is much less than the 60l.2 seconds needed for the corresponding Vt-optimal run order. It follows that at least for h = 27 and a second-order time trend, the adjustment algorithm is very suitable to compute good run orders when the computation time is an important issue. A similar conclusion can be drawn for the other time trends. Table 4: Trend factors and computation times of the Vcoptimal run orders and the adjusted run orders for the full second-order response model in three factors, n = 27 and h = 27 DDt gl (t) g2(t) g3(t) g4(t) (d = 27) l.0000 0.9217 0.9202 0.8690 trend factor DDt (d = 125) l.0000 0.9637 0.9501 0.9189 Dadj l.0000 0.9536 0.9410 0.9045 DDt computation time (sec) (d = 27) DDt (d = 125) Dadj 3.1 6.2 7.6 7.7 274.8 60l.2 567.0 544.7 3.1 6.2 7.6 7.7 The next example compares the adjustment of the Vt-optimal run orders computed with d = 27 and h = 27 with the Vroptimal run orders for d = 27 and h = 53. The former ones are used as input to the adjustment algorithm and to make the comparison valid, only the time points are allowed to be adjusted. Besides, the parameters of the adjustment algorithm are specified such that the time points of the adjusted run orders form part of the set of time points of the Vroptimal run orders for h = 53. The results are given in table 5. For instance, the trend factor of the adjusted run order for a quadratic time trend is equal to 0.9536, which is again very close to trend factor 0.9638 of the corresponding Vcoptimal run order with h = 53. The computation time of the adjusted run order is much lower than that of the Vcoptimal run order, namely 6.2 seconds is almost negligible as compared to 131. 7 seconds. Similar conclusions hold for the other trend functions. As a final illustration consider the comparison between run orders with possibly replicated observations and run orders without replicated observations. The number of design points and time points equals d = 27 and h = 53 respectively. The results are given in table 6. Taking into account the possibility for replicated observations naturally involves an improvement in the balance for time trend effects. For instance, the trend factor of the Vt-optimal run order for gl(t) without replicated observations is equal to 0.9413, 16 Table 5: Trend factors and computation times of the 'Droptimal run orders and the adjusted run orders for the full second-order response model in three factors, n = 27 and d = 27 DDt gl (t) g2(t) g3(t) g4(t) trend factor (h = 27) DDt (h = 53) 1.0000 1.0000 0.9217 0.9638 0.9202 0.9501 0.8690 0.9193 Dadj DDt 1.0000 0.9536 0.9410 0.9075 computation (h = 27) DDt 3.1 6.2 7.6 7.7 time (sec) (h = 53) 38.2 131.7 99.4 100.3 Dadj 3.1 6.2 7.6 7.7 whereas allowing replicated observations leads to a completely linear-trend free run order. However, taking into account replicated observations goes at the cost of the computation time. Applying the adjustment algorithm-with the parameters set to allow replicated observations-to the 'Droptimal run order without replicates augments the trend factor to 0.9926, which approximates complete trend-resistance. The computation time of the adjusted run order is again much lower than that of the 'Droptimal run order with replicated observations. The conclusions can be generalized to the higher-order time trends. Table 6: Trend factors and computation times of the 'Droptimal run orders and the adjusted run orders for the full second-order response model in three factors, n = 27, d = 27 and h = 53 DDt gl (t) g2(t) g3(t) g4(t) trend factor (no rep!.) DDt (rep!.) 0.9413 1.0000 0.9090 0.9638 0.9501 0.8959 0.8682 0.9193 Dadj 0.9926 0.9544 0.9327 0.8978 DDt computation time (sec) (no rep!.) DDt (rep!.) Dadj 3.8 3.8 38.2 14.0 14.0 131.7 14.2 14.2 99.4 11.3 100.3 11.3 Summing up, this example has clearly illustrated that it is recommended to use the adjustment algorithm in addition to the exchange procedure when the computation time is an important issue. The resulting outperformance in computation time goes at the expense of the trend-resistance but the loss is only moderate. 17 6 Conclusion This paper has proposed an adjustment algorithm to improve the trend factors of the 'Droptimal run orders computed with the exchange algorithm. This is done by moving the design points and/or the time points of the predefined candidate set along their axes. Only design modifications that improve the 'Dt-criterion value are accepted. The computational results have shown that an adjustment of the design points and/or the time points considerably increases the protection against time dependent results. Besides, it has been illustrated that the assumption of equally spaced time points is in general not a good one and that the adjustment algorithm is very useful to compute the optimal time points. Finally, instead of using a conventional exchange procedure to compute run orders over a fine grid, the use of the adjustment algorithm in addition to an exchange procedure over a coarse grid is very recommendable in situations where the aim is to obtain good run orders at an acceptable computation time. 18 Appendix. The design algorithm In the outline of the algorithm, the starting run order is written as R = {(Xi, tjl) and the criterion value will be denoted as Q. The initial step lengths for the design points and the time points are denoted as Sl and S2 respectively. The minimum step lengths are written as 31,min and 32,min respectively. Finally, Vi = 1 indicates that design points are considered to be adjusted, whereas the opposite holds for Vi "# 1. In a similar way, the value of parameter V2 indicates whether the time points are considered for adjustment or not. The values of Vi and V2 are of course user-specified. The practitioner has the additional possibility to impose a minimum distance to be maintained between any two successive time points. To preserve clarity, this option is left out from the outline. The output of the algorithm consists of the adjusted run order Ropt and its corresponding criterion value Qopt. After reading the input, the algorithm proceeds as follows: 2. Compute the criterion value Q of the starting run order R. 3. Set Qopt = Q and R opt = R. 4. Adjust the run order: (a) Set 6 1 = 1. (b) If Vi = 1, then find the best change in design points: VXi E R, Vj E {I, ... ,J}: • If Xij - 31 ;::: -1, compute the effect 6 on Qopt of changing If 6> 6 1 then set 6 1 = L, Wll = -1, al = i and bl = j. • If Xij + 31 ::; 1, compute the effect 6 on If 6 > 6 1 then set 6 1 = 6, Wll = 1, al (c) Set 6 (d) If 2 V2 = Vt", E • If If • If If (e) Set 6 Qopt =i of changing and bl = j. Xij Xij to to Xij - Xij 31. + 31. = 1. 1, then find the best change in time points: R: -1, compute the effect then set 6 2 = 6, W22 = tk + 82 ::; 1, compute the effect 6 6 > 6 2 then set 62 = 6, W22 = tk - 32;::: 6> 6 3 = 2 6 on Qopt of changing tk to the - 32· -1 and C2 = k. on Qopt of changing tk to h + 32· 1 and C2 = k. 1. (f) If Vi = 1 and points: V2 = 1, then find the best change in design points and time V(Xi' tk) E R, Vj E {I, ... ,J}: 19 • If Xij-81 2: -1 and tk-82 2: -1, compute the effect I:::. on Qopt of changing Xij to Xij - 81 and tk to tk - 82· If I:::. > 1:::.3 then set 1:::.3 = 1:::., W31 = -1, W32= -1, a3 = i, b3 = j and C3 = k. • If Xij - 81 2: -1 and tk + 82 :'S: 1, compute the effect I:::. on Qopt of changing Xij to Xij - 81 and tk to tk + 82. If I:::. > 1:::.3 then set 1:::.3 = 1:::., W31 = -1, W32 = 1, a3 = i, b3 = j and C3 = k. • If Xij + 81 :'S: 1 and tk - 82 2: -1, compute the effect I:::. on Qopt of changing Xij to Xij + 81 and tk to tk - 82· If I:::. > 1:::.3 then set 1:::.3 = 1:::., W31 = 1, W32 = -1, a3 = i, b3 = j and C3 = k. • If Xij + 81 :'S: 1 and tk + 82 :'S: 1, compute the effect I:::. on Qopt of changing Xij to Xij + 81 and tk to tk + 82. If I:::. > 1:::.3 then set 1:::.3 = 1:::., W31 = 1, W32 = 1, a3 = i, b3 = j and C3 = k. 5. Compute I:::. = max(1:::. 1 , 1:::. 2 , 1:::.3)' 6. If I:::. > 1 then + W1l81 (a) If I:::. = 1:::. 1 , (b) If I:::. = 1:::. 2 , then replace te2 with te2 + W2282 and go to (d). = 1:::. 3 , then replace Xa3 b3 with Xa3b3 + W3181 and te3 with te3 + W3282. (c) If I:::. then replace Xa,b, with Xa,b, and go to (d). (d) Update Qopt. (e) Go to 4. 7. Reduce the step lengths: (a) If V1 = 1, V2 =f. 1 and i. Set 81 = 81/2. ii. Go to 4. (b) If V1 =f. 1, V2 = i. Set 82 = ii. Go to 4. (c) If V1 = 1, V2 Qopt and then 82/22: 82,min, then (81/22: 81,min or 82/2. = 1 and i. If s1/2 2: ii. If 82/2 2: iii. Go to 4. 8. Write 1 and 81/2 2: 81,min, Sl,min 82,min then set then set 82/2 2: 82,min), 81 = 81/2. 82 = 82/2. then Rapt. The algorithm is implemented in Fortran 77 and makes use of the library Netlib of Bell Labs. It uses predefined routines to factor a symmetric matrix, to compute the determinant of a factored symmetric matrix, to invert a symmetric matrix and to generate random numbers. 20 References Atkinson, A. and Donev, A. (1992). Optimum Experimental Des£gns, Oxford: Clarendon Press. Atkinson, A. and Donev, A. (1996). Experimental designs optimally balanced for trend, Technometrics 38: 333-341. Box, G. and Draper, N. (1971). Factorial designs, the matters, Technometrics 13: 731-742. IX/XI criterion and some related Cheng, C.-S. (1985). Run orders of factorial designs, In L. LeCam and R. Olshen (eds), Proceedings of the Berkeley Conference in Honor of Jerzy Neyman and Jack Kiefer, Wadsworth, pp. 619-633. Cox, D. (1951). Some systematic experimental designs, Biometrika 38: 312-323. Daniel, C. and Wilcoxon, F. (1966). Fractional 2P - Q plans robust against linear and quadratic trends, Technometrics 8: 259-278. Donev, A. and Atkinson, A. (1988). An adjustment algorithm for the construction of exact V-optimum experimental designs, Technometrics 30: 429-433. Farell, R., Kiefer, J. and Walbran, A. (1967). Optimum multivariate designs, Proceedings of the 5h Berkeley Symposium, University of California Press, Berkeley, pp. 113-138. Freeny, A. and Lai, W. (1997). Planarization by chemical mechanical polishing: a rate and uniformity study, In V. Czitrom and P. Spagon (eds), Statistical Case Studies for Industrial Process Improvement, Philadelphia: ASA-SIAM, pp. 251-264. Goos, P. (2001). The optimal design of blocked and split-plot experiments, Ph.D. dissertation, Katholieke Universiteit Leuven, Belgium. Harville, D. (1997). Matrix Algebra from a Statistician's Perspective, New York: SpringerVerlag. Joiner, B. and Campbell, C. (1976). Designing experiments when run order is important, Technometrics 18: 249-259. Tack, L. and Vandebroek, M. (2001). (V t , C)-optimal run orders, Journal of Statistical Planning and Inference 98: 63-80. 21