Robustness of response surface designs to missing data

In this paper, minimax loss response surface designs are constructed. These designs are more robust to one missing design point than the original designs. The proposed designs are compared with the designs in the literature, and they are better in terms of loss and number of runs. Moreover, the new suggestion for the value of α generates designs not only with less losses but also with higher D‐efficiency.


INTRODUCTION
Missing observations, even in well-planned experiments, can rarely be avoided. Researchers have discovered ways of handling problems associated with missing observations. These include methods either replace the missing values or to apply robust techniques to their analysis. The term "robust" that was first used by Box 1 is often encountered when studying the effects of missing observations. According to Wendelberger, 2 robustness is crucial in the real world of experimentation; thus, the need for experimental designs that are robust in the presence of missing observations arise.
Many researchers have discussed the robustness of statistical designs against missing observations. Hedayat and John 3 and Ghosh 4 considered robust balanced incomplete block (BIB) designs. Herzberg and Andrews 5,6 studied the problem of missing observations in a design from a different angle. They considered the probability of losing an observation at a design point. In a later study, Andrews and Herzberg 7 suggested maximizing the expected value of the determinant of the information matrix to compensate for possible missing observations. Ahmad et al 8 considered multilevel augmented second-order response surface designs and their robustness to missing observations.
In investigating the information contained in an observation, Ghosh 9 stated that in most designs, some design points were more informative than others and that the effects of one or more design points, when missed, attracted a higher loss in efficiency. Other economical designs requiring fewer design points were constructed and studied for their robustness to missing observations by Ahmad and Akhtar. 10 Angelopoulos et al 11 constructed new central composite designs (CCDs) that are more economical than the classical CCDs that were introduced by Box and Wilson 12 since the newer designs were built with fewer runs. Using their designs, they found that the minimum number of runs for a four-factor design was 18. Another economical design was constructed by Georgiou et al. 13 They constructed a class of composite designs with a higher D-value (see the study of Lucas 14 ) than Angelopoulos's designs but with the same number of runs.
Various robust-to-missing-observations criteria have been proposed in the literature. Akhtar and Prescot 15 developed the minimax loss criterion and applied it to the evaluation and generation of CCDs. They found that, when they applied minimax loss criterion on the classical CCDs with four factors, the number of runs was 26. Since then, this criterion has been adopted by many researchers including Ahmad and Gilmour, 16 who used it to study the robustness of subset designs. In addition, the hat matrix, also known as the projection matrix, performs a major role in quantifying the effect of removing one or more observations from a design. 17 Box and Draper 18 pointed out the connection between the diagonal elements of a hat matrix and a design's robustness and applied their related criterion to CCDs. Akram 19 utilized the hat matrix (through its compounds) in studying the robustness of CCDs up to three missing observations.
Each robust-to-missing-observations criterion has its limitation(s). One limitation is that they cannot be applied with the remaining runs if these runs are not enough to estimate the parameters of interest. Another limitation is that each criterion focuses on one area of interest such as loss, D-efficiency, estimation capacity, and so on. Therefore, more flexible robust-to-missing-observations criteria may be required to circumvent such restrictions. In the literature, there is nothing on the effect of applying minimax loss criterion on the latest CCDs that are constructed with a minimum number of runs. This paper integrates the minimax loss criterion with the latest CCDs proposed by Georgiou et al 13 denoted as GSA, and Angelopoulos et al 11 denoted as AEK, to achieve more robust designs in the case of missing observations. The GSA and AEK designs are used in this paper because these designs have less runs than the classical CCDs that were evaluated by Akhtar and Prescott. 15 In section 2, the construction of the GSA and AEK designs are discussed. In section 3, the minimax loss criterion is discussed. Section 4 details the methodology. The new results and a comparison of our findings are presented in section 5.

GSA AND AEK
In this section, we present the proposed construction methods of the GSA and AEK designs that are used in this paper.

GSA designs
The GSA designs are constructed as a three-part design matrix X that would be suitable for fitting the second-order response surface as follows: 1. For the factorial part, select a suitable number of balanced columns from a known two-level design with good properties. For example, one may use columns from a Hadamard matrix of order n f , from a fractional factorial design with n f runs, from a supersaturated design with a run size of n f , or from other suitable two-level design where n f is to be selected. They called this part of the design F. 2. Instead of using the traditional axial points, they suggested using vertices of a hypercube that circumscribe the sphere of zero center and radius . This approach includes the axial points as a special case. They called this part of the design V. 3. Select a number of center points n c , and this part of the design is denoted as C.
The desirable GSA design matrix is as follows: X = (F', V', C',-V) Example 1. Suppose we wish to construct a six-factor GSA design. Select six columns for part F from the Hadamard matrix of order 20, and for part V, select six columns of the weighing matrix of order 8 (W 8 ) multiplied by . More details on weighing matrices can be found in Koukouvinos and Seberry. 20

AEK designs
The AEK designs are constructed by an exhaustive computer search. The aim was to find the minimum number of suitably selected points that are needed for the fractional factorial part (four factors to nine factors) so that the designs estimate all the model parameters. The classical method of CCDs is used for the axial part. Table 1 present the general structure of the AEK designs.

Example 2.
Using a computer search for the factorial part, the authors provide the following design of six factors that maximize the criterion.

MINIMAX LOSS CRITERION
Akhtar and Prescott 15 studied the reduction in the determinant of the information matrix because of missing observations. As previously mentioned, they developed the minimax loss criterion to find designs robust to missing observations. The effect of missing observations is quantified using the criterion of the loss function proposed by Akhtar and Prescott. 15 More specifically, the loss function gives a measure of the relative reduction in the determinant of the information matrix associated with the complete design in the presence of missing observations. More studies on the loss function consider second-order CCDs differing in the number of control variables, configurations of the factorial, and axial and centre portions. For a p parameters, p = ((K + 1)(K + 2))∕2, second-order model written in matrix form as follows: where K is the number of factors, y denotes the vector of observations, X denotes the model matrix, denotes the vector of unknown parameters that are estimated based on N uncorrelated observations and where denotes the random error component assumed to be normally and independently distributed with a mean of zero and a constant variance. The N×K design matrix, D, comprising N design points, is expanded into N×p model matrix X. The columns of X are obtained using the design matrix and the full second-order model, where the number of model parameters p equals ((K + 1)(K + 2))∕2.
When m observations in the design are missing, the model matrix reduces by m rows. As expected, the information matrix, X ′ r X r , for the reduced design will be different from the information matrix, X ′ X, for the complete design. If X m is the m × p matrix of the m missing rows corresponding to the m missing observations, then the information matrix may be expressed as follows: Post multiplying by (X ′ X) yields so the determinant will be Akhtar and Prescott 15 call the diagonal element of the m th compound of (I − R) where R is the hat matrix according to the m missing observations and I is an identity matrix of the same dimension as R.
For a square matrix H of order n, the m th compound of H is a square matrix of order the elements of which are the minors of determinant |H| of order m. For a particular number of design points n and number of parameters p to be estimated, the sum of diagonal elements of the m th compound is constant. For more detail, see Aitken and Rutherford. 21 Andrews and Herzberg 7 presented a relation between |X ′ r X r | and |X ′ X| by using identities on the expansion of a bordered determinant (see the study of Rao 22 ) as follows: where A uv … is the corresponding diagonal element of the m th compound of (I − R) matrix and R is the symmetric matrix X(X ′ X) −1 X ′ of order N. Specifically, A u is the u th diagonal element of the first compound of (I − R) matrix, A uv … is the vu th diagonal element of the second compound of (I − R) matrix and so on. For m missing observations, the determinant of the complete information matrix may be obtained as follows: The relative loss because of any set of m missing observations is defined as follows: It was found that Under D-optimality, the determinant |(X ′ X)| is maximized. It is well established for CCDs that the determinant |(X ′ X)| is an increasing function of the axial distance and is maximized at = ∞. The value of specifies the location of the axial point and may be chosen to satisfy various design conditions such as rotatability and orthogonal blocking.
In the presence of missing observations, the determinant of the information matrix |(X ′ X)| of a complete design is reduced to | | (X ′ r X r ) | | . It is essential that this reduction is as small as possible. Thus, the relative loss, l uv … needs to be minimized, and this is possible if tends to one. The most practicable way to minimize the losses is to minimize the maximum loss because of a set of m missing observations. Minimizing the maximum loss is the idea behind the minimax loss criterion of Akhtar and Prescott. 15

The choice of for the new designs
Alpha ( ) is used to satisfy or maximize some conditions such as the rotatability of the design or to maximize the determinant of the information matrix. Fang and Mukerjee 23 discussed the consequences of rescaling by restricting themselves to a cuboidal region [−1, 1]. They concluded that = 1 is the best choice under the D-criterion. Here, the condition of the cuboidal design region was relaxed, and the designs were constructed by finding a new value of using the minimax loss criterion presented by Akhtar and Prescott. 15 This criterion is propitious for the construction of designs that are more robust to missing observations. The value of was chosen so that the maximum loss of a missing design point is minimized.

The detailed procedure for constructing the proposed designs
We now present the detailed procedure for constructing the proposed designs. Algorithm Step1. Choose an AEK or GSA design X = , where F is the factorial part and V is the axial part. Note that center points may be added to the chosen design.
Step2. Calculate the loss function L f , L a , and L 0 based on X for the second order model (1), where L f is the loss for a missing factorial point, L a is the loss for a missing axial point, and L 0 is the loss for a missing centre point. Step3. Find that minimizes the maximum loss of { L , L a , L 0 } Step4. Substitute the new in the chosen design.

THE NEW RESULTS
In this section, we describe the new results we obtain by exhaustively searching the values of that minimize the maximum loss.

Example 3.
An example of how to find that minimizes the maximum loss for 6 factors with one missing observation is given in Table 2. The general structure of the design matrix is presented in Table 3. It can be seen in Figure 1 that L f is decreasing and L a is increasing due to the function of with a negative and positive slope on the graph. In order to compute the value of at which the maximum loss is minimized, function L f is equated with L a , and it was observed that 1.335 was the minimum value of at which L f is equal to one of the other functions. Thus, the maximum loss of missing one design point will be minimized at = 1.335 with loss = 0.89.  Wilson, 12 and their values were investigated by Akhtar and Prescott. 15 Small composite designs introduced by Draper and Lin, 24 and their values are studied here for the first time. Augmented pairs designs introduced by Morris, 25 and their values were investigated by Ahmad et al. 8 These designs introduced by Angelopoulos et al, 11 and their values are studied here for the first time. These designs introduced by Georgiou et al, 13 and their values are studied here for the first time.  Table 4 shows the values of maximum loss for one missing design point with the number of runs in parentheses and the values of that we used to construct new designs that are robust against a missing observation. From Table 4, it is obvious that the newly proposed designs, AEKM, which are the minimax designs derived from AEK and GSAM that are the minimax designs derived from GSA, perform better in comparison to at least one of the three classical designs that are used in the literature on CCD, small composite design (SCD), and augmented pair (AP) designs. For example, the five-factor design of GSA with 26 runs and one centre point has the lowest value of maximum loss for one missing point of all the designs in Table 4. Furthermore, it had a fewer number of runs, which makes it an economical and robust design for one missing observation. Although the CCD is superior in most cases, the newly proposed designs show improvement over one or more of the other two designs. The designs given in bold face provide the minimum number of runs for estimating p parameters with less loss than the designs in the literature in the case of one missing observation. In summary, the new  results on Table 4 are K = 4 with 21 runs, K = 5 with 27 runs and 31 runs, and K = 6 with 35 and 41 runs. Also, there is new improved design for K = 7 with 41 runs. Table 5 shows the relative D-efficiency of the robust designs with the new values that are presented in bold in Table 4. Following the literature, comparison is presented in the cases of one centre point and five centre points. Since we investigated all possible values of to find the one that gives the less minimax loss when a design point is missing, then it makes sense to see whether we lost D-efficiency by trying to optimize the design based on the loss function. Thus, we think that a comparison of the relative D-efficiency of designs with the values suggested in the literature and the new values suggested here is relevant. Table 5 shows that the designs with minimum losses, AEKM and GSAM, also perform better in terms of the relative D-efficiency when they are compared with the classical AEK and GSA designs (with the used in the literature).
The D-efficiency of a design is defined as follows: where X is the corresponding model matrix that is defined in (2) . The relative D-efficiency of GSAM is defined to be The relative D-efficiency of AEKM is defined as follows: The values in Table 6 are selected so that the design gives the minimax loss. For the values of CCD and AP that give the minimax loss, see Akhtar and Prescott 15 and Ahmad et al, 8 respectively. For SCD, AEKM, and GSAM, we investigated exhaustively and found the values of that give the minimax loss. So Table 6 compares the D-efficiency of different designs using the values that make them robust against a missing observation. The purpose is that all the designs in Table 6 are robust against a missing observation, so D-efficiency gives a secondary evaluation of these robust designs. It resulted in most of the generated robust designs, in particular the AEKM designs, having higher D-efficiency than the designs in the literature. For instance, the six-factor design of the AEKM appears to have higher D-efficiency, even though it has less runs. The six-factor design is an improvement over all the known designs in the literature with same number of runs. Moreover, this design overtakes the other known designs in the literature even if those have more runs. In addition, if we compare the five-factor design of AEKM with AP and SCD, we can see that AEKM's D-efficiency is higher than both AP and SCD even although its number of runs is equal to SCD and less than AP design. GSAM designs have less D-efficiency than the other designs except in the case of five factors and 26 runs where this design is superior to all the designs in the literature.

DISCUSSION
Robust designs are essential in the case of missing observations as even the most well-planned experiments are susceptible to missing data. Choosing a response surface design to fit certain kinds of models is a difficult task. However, a design that performs well when an experimental observation is missing is always the most desirable. In this article, new designs were presented as a variation of the designs provided in GSA and AEK, which are more robust in the case of missing observations. The modified designs used the minimax loss criterion to choose a value of for a whole set of axial parts in the experimental region. The newly proposed designs showed better performance than the GSA and AEK designs and occasionally performed better than CCD, SCD, and AP when the restriction to a cuboidal experimental region was relaxed. These designs are also comparable with standard CCDs because of their five levels for each factor. CCDs are still superior in most cases, but because of a significant reduction in the design size, for larger K, the newly proposed designs may be more preferable in some cases. Even though the newly proposed designs are introduced as single-stage designs, their structure suggests their potential for use in sequential experimentation as well.
Additional interested problems in this field would be to investigate the changes in losses and value in the case of more than one missing observations. That interesting and related problem might be of particular interest and will be a generalization of this work. Furthermore, rotatability is another additional criterion that can be impended in studying response surface designs but this will requires extended effort as different missing observations may resulted at the same losses bur provide different rotatability properties.
Generally, when researchers are dealing with missing observations, they need to be extremely careful trying to handle simultaneously optimization of several criteria. The multiple stage analysis will be affected by the initial screening stage and the overall efficiency will also be modified.