UFDC Home  Search all Groups  UF Institutional Repository  UF Institutional Repository  UF Theses & Dissertations  Vendor Digitized Files   Help 
Material Information
Subjects
Notes
Record Information

Full Text 
A MEAN SQUARED ERROR OF PREDICTION APPROACH TO THE ANALYSIS OF THE COMBINED ARRAY By EILEEN M. O'DONNELL A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1994 ACKNOWLEDGEMENTS I would like to express my sincere gratitude to Dr. Geoff Vining for being my ad visor and for all the time and effort he put into helping me complete this dissertation. His patience, encouragement, and guidance during the sometimes slow and frustrat ing periods of independent study enabled me to get to the finish line. Dr. Vining has taught me much more than I could have learned on my own, and I look forward to continued collaboration with him during the years to come. I would also like to thank Dr. John Cornell, Dr. Ronald Randles, and Dr. Jack Elzinga for their support while serving on my committee. Thanks also go to Dr. Brett Presnell for attending my oral exams. There are many other teachers to whom I wish to express my gratitude for being role models and mentors during my education, from elementary school through graduate school. In addition, I would like to thank my family: Mom, Dad, Steve, Ed, and Tracey, who are responsible for the person I am today. My parents placed great value on our education and we have each achieved great success because of it. Special thanks are also owed to my classmate, Kannan, who at any time would drop whatever he was doing to help me with my latest "urgent" problem. There are many other friends who made my years in graduate school memorable. I would also like to thank Carol Rozear, who has been extremely helpful in dealing with the many nonacademic details during my years here and for letting me keep my jet skis at her home. Lastly, I wish to thank David Wolfe for his kindness, patience, and love and for providing me with muchneeded distractions. David will forever be my best friend. TABLE OF CONTENTS ACKNOW LEDGEM ENTS ..................................................... ii LIST OF FIGURES ........................... .... ......................... v A B ST R A CT ................................................................... viii CHAPTERS 1 INTRODUCTION .................................................... 1 2 REVIEW OF LITERATURE ......................................... 6 2.1 A Review of Taguchi's Parameter Design Methods ............ 6 2.2 The Combined Array: An Alternative to Taguchi's Crossed Array .......................................................... 13 2.3 Run Savings and Estimation Flexibility of the Combined Array 17 2.4 Analysis of a Combined Array Experiment ..................... 19 2.5 The Mean Squared Error Criterion ............................. 21 3 ANALYSIS OF THE COMBINED ARRAY WHEN THE ASSUMED MODEL IS CORRECT ................. 24 3.1 Introduction ..................................................... 24 3.2 Properties of the Estimated Process Variance .................. 25 3.3 Expressions for the Bias and Variance of the Estimated Process Variance ....................................................... 28 3.4 Comparing Designs Based on the Mean Squared Error Criterion 30 3.5 Results and Discussion .......................................... 40 3.6 Other Measures of Design Performance ....................... 46 3.7 Practical Application ........................................... 53 3.8 Conclusions ..................................................... 55 4 ANALYSIS OF THE COMBINED ARRAY WHEN THE ASSUMED MODEL IS UNDERSPECIFIED ....... 59 4.1 Introduction ..................................................... 59 4.2 Properties of the Estimated Process Variance .................. 60 4.3 Expressions for the Bias and Variance of the Estimated Process V ariance ....................................................... 66 4.4 The IMSE Criterion and Other Measures of Design Performance 67 5 SPECIAL CASES OF MODEL UNDERSPECIFICATION ......... 69 5.1 Introduction ..................................................... 69 5.2 Reduced Expressions for the Bias of the Estimated Loss ....... 70 5.3 Exam ple ......................................................... 73 5.4 Results and Discussion ....................................... 75 5.5 Summary and Conclusions ..................................... 113 5.6 Practical Application ................ .......................... 116 6 CONCLUSIONS AND FUTURE RESEARCH ...................... 126 APPENDICES A INVARIANCE OF MSE ............................................. 129 B BIAS DERIVATION .................................................. 132 C VARIANCE DERIVATION .......................................... 136 D DERIVATION OF COVARIANCE MATRIX ........................ 138 E SAS PROGRAM ..................................................... 143 REFERENCES ................................................................. 147 BIOGRAPHICAL SKETCH .................................................. 149 LIST OF FIGURES Figure 2.1 Illustration of a Crossed Array ............................... 3.1 16run Resolution V Fractional Factorial Design ................. 3.2 20run PlackettBurman Design .............................. 3.3 24run Plackett Burman Design .............................. 3.4 24Run Orthogonal Array ................................... 3.5 The Effect of Correlation 3.6 The Effect of Correlation 3.7 The Effect of Correlation 3.8 The Effect of Correlation 3.9 The Effect of Correlation 3.10 The Effect of Correlation 3.11 The Effect of Correlation 3.12 The Effect of Correlation (Unweighted Weighted) 3.13 The Effect of Correlation (Unweighted Weighted) P the Average Bias (Unweighted Weighted) the Average MSE .................. the Weighted Average MSE .......... the MaximumAverage MSE .......... the Weighted MaximumAverage MSE . the MaximumMaximum MSE ........ the Weighted MaximumMaximum MSE the Average Bias and Average MSE .............................. the Average Bias and Average MSE the Average Bias and Average MSE Case 1 (C2) The Effect of Correlation on the Average Bias ........ Case 1 (C2) The Effect of Correlation on the Weighted Average Bias Case 1 (C2) The Effect of Correlation on the Average MSE ....... Case 1 (C2) The Effect of Correlation on the Weighted Average MSE Case 2 (N2) The Effect of Correlation on the Average Bias........ age 9 34 38 39 41 43 44 45 48 49 51 52 56 57 77 78 79 80 82 5.6 Case 2 (N2) The Effect of Correlation on the Weighted Average Bias 5.7 Case 2 (N2) The Effect of Correlation on the Average MSE ....... 5.8 Case 2 (N2) The Effect of Correlation on the Weighted Average MSE 5.9 Case 3 (C2 x N) The Effect of Correlation on the Average Bias .... 3 (C2 x N) 3 (C2 x N) 3 (C2 x N) 4 (C2 x N2) 4 (C x N2) 4 (C x N2) 4 (C x N2)  The Effect of Correlation on the Weighted Average S. . . . . 88  The Effect of Correlation on the Average MSE ... 89  The Effect of Correlation on the Weighted Average S. . . . 90  The Effect of Correlation on the Average Bias .... 92  The Effect of Correlation on the Weighted Average . . . . . 93  The Effect of Correlation on the Average MSE 94  The Effect of Correlation on the Weighted Average . .. .. .. .. .. .. 95 5.17 Case 5 (C2, N2) The Effect of Correlation on the Average Bias ..... 5.18 Case 5 (C2, N2) The Effect of Correlation on the Weighted Average Bias .................................................. 5.19 Case 5 (C', N2) The Effect of Correlation on the Average MSE . 5.20 Case 5 (C2, N2) The Effect of Correlation on the Weighted Average MSE ................................................... 5.21 Case 6 (C2, C2 x N) The Effect of Correlation on the Average Bias . 5.22 Case 6 (C2, C2 x N) The Effect of Correlation on the Weighted A average Bias ........................................... .. 5.23 Case 6 (C2, C2 x N) The Effect of Correlation on the Average MSE 5.24 Case 6 (C2, C2 x N) The Effect of Correlation on the Weighted A average M SE ............................................ 5.25 Case 7 (N2, C x N2) The Effect of Correlation on the Average Bias . 5.26 Case 7 (N2, C x N2) The Effect of Correlation on the Weighted A average B ias. ............................................ 5.27 Case 7 (N2, C x N2) The Effect of Correlation on the Average MSE 5.10 Case Bias 5.11 Case 5.12 Case MSE 5.13 Case 5.14 Case Bias 5.15 Case 5.16 Case MSE 5.28 Case 7 (N2, C x N2) The Effect of Correlation on the Weighted Average M SE ............................................ 110 5.29 Case 8 (C2, N2, C2 x N, C x N2) The Effect of Correlation on the Average Bias ......................................... 111 5.30 Case 8 (C2, N2,2 x N,C x N2) The Effect of Correlation on the W eighted Average Bias ..................................... 112 5.31 Case 8 (C2,N2, 2 x N,C x N2) The Effect of Correlation on the Average M SE ............................................ 114 5.32 Case 8 (C2, N2, C2 x N,C x N2) The Effect of Correlation on the W eighted Average M SE .................................... 115 5.33 Case 1 (C2) The Effect of Correlation on the Average Bias and Average MSE (Unweighted Weighted) ........................ 118 5.34 Case 2 (N2) The Effect of Correlation on the Average Bias and Average MSE (Unweighted Weighted) ........................ 119 5.35 Case 3 (C2 x N) The Effect of Correlation on the Average Bias and Average MSE (Unweighted Weighted) ........................ 120 5.36 Case 4 (C x N2) The Effect of Correlation on the Average Bias and Average MSE (Unweighted Weighted) ........................ 121 5.37 Case 5 (C2, N2) The Effect of Correlation on the Average Bias and Average MSE (Unweighted Weighted) ........................ 122 5.38 Case 6 (C2, C2 x N) The Effect of Correlation on the Average Bias and Average MSE (Unweighted Weighted) .................... 123 5.39 Case 7 (N2, C x N2) The Effect of Correlation on the Average Bias and Average MSE (Unweighted Weighted) .................... 124 5.40 Case 8 (C2, N2, C2 x N,C x N2) The Effect of Correlation on the Average Bias and Average NISE (Unweighted Weighted) .......... 125 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy A MEAN SQUARED ERROR OF PREDICTION APPROACH TO THE ANALYSIS OF THE COMBINED ARRAY By Eileen M. O'Donnell August 1994 Chairman: G. Geoffrey Vining Major Department: Statistics The combined array is an experimental design used for quality improvement within a robust parameter design framework. Robust parameter design is concerned with how to reduce the variation of a product's performance. In this framework, the factors that influence the quality characteristic of interest are of two types: control factors and noise factors. Control factors are those that are under the direct control of the experimenter, both in the experiment and in the process. Noise factors can be controlled for the purposes of an experiment, but are difficult to control in the process. The goal of a robust parameter design experiment is to select the values of the control factors that make the product insensitive to changes in the noise factors. The combined array approach represents an effort to integrate Taguchi's parameter design principles with more wellestablished and efficient statistical techniques. This approach uses a single design matrix and models the response directly as a function of the control and noise factors. The fitted model is used to determine settings of the control factors that will minimize an appropriate loss function involving the noise factors. Often, the most appropriate loss function is simply the resulting process variance, recognizing that the noise factors, although fixed in the experiment are actually random effects in the process. When the focus of an experiment is to find control factor settings that achieve a target with minimum variability, it is vital to understand the properties of the process variance estimator. This dissertation represents the first formal attempt at deriving such properties. The derivation of the properties of the process variance estimator is dependent upon the model associated with the combined array. The mean squared error for the estimated process variance for the combined array approach is developed for two cases: (1) when the model is correctly specified and (2) when the model is underspecified. Specific combined arrays are evaluated graphically through use of a mean squared error criterion. Simulation results and additional practical examples outline how this approach may be used to select an optimal combined array within a particular experimental situation. CHAPTER 1 INTRODUCTION The quality control movement in the 1990s has recognized the need to do more than find operating conditions that achieve a target value for the mean in order to improve product quality. Today's industrial statisticians are encountering a need for developing experimental strategies for achieving some target condition for a charac teristic of interest while simultaneously minimizing its variance. To accomplish this, it is assumed that the experimental factors under consideration separate into a set of control factors and a set of noise factors. Control factors are those variables under the direct control of the experimenter. Noise factors are often environmental variables such as ambient temperature and humidity, properties of raw materials, product ag ing, etc. In certain applications they measure how the consumer uses or handles a product. These noise variables can be controlled by the experimenter during the pilot plant phase or at the research and development level, but they are difficult to control at the production level or out in the field. To motivate this type of optimization problem, let us consider an example dis cussed by Box and Jones (1992). Suppose a manufacturer is seeking the best cake recipe for a cake mix to be sold in a box at the supermarket. Preliminary experi mentation has identified three important control factors to be flour (F), shortening (S), and egg powder (E) and the ranges of F, S, and E that are likely to produce a goodtasting cake if the mix is baked at a specified oven temperature (T) for a given length of time (t). Now, it is anticipated that the consumer may have an oven in which the temperature is biased up or down and that the consumer may overcook or undercook the cake. Thus, the manufacturer desires a robust recipe so the cake will taste reasonably good even when the time and temperature of baking, which are noise variables in this case, differ somewhat from the recommended levels on the box. Hence, the job of the statistician during the product design stage of the cake mix is to provide an experimental strategy for finding the settings of the three control variables that are most robust to changing levels of the noise variables. If we assume the response will be taste measured on a scale from 1 to 7, with 7 being best, then we aim to find settings of E, S, and F which achieve a score near 7 with minimum variability. The cake mix optimization problem falls under the experimental design method ology known as robust parameter design, in which the experimenter attempts to an ticipate the maximum amount of noise expected in the process and then finds those settings of the control variables most robust to the maximum noise. We note that this robustness to noise variables can only be incorporated into the product design stage, a fact emphasized by Genichi Taguchi, a pioneer in the methodology for robust parameter design. Taguchi was among the first practitioners to formalize the idea that as process variables change, process performance variability changes, and this change must be incorporated into the statistical design and analysis. In order to im prove products and processes, Taguchi recognized that it was no longer sufficient to concentrate on the mean response alone. In addition, the practitioner must include product variability as a performance response. In a first approach to experimentation using both control and noise factors, the experimental plans introduced by Taguchi have become known as "crossed arrays." These plans cross one experimental design for the control factors with a separate design for the noise factors. Crossed arrays have been criticized for requiring a large number of experimental runs, often wasting resources on the estimation of higher order controlbynoise (C xN) interactions, and not allowing estimation of control bycontrol (C x C) or noisebynoise (Nx N) interactions. To address these and other criticisms and disputes concerning Taguchi's approach, a number of authors have proposed alternative methods which seek to unite Taguchi's principles with more traditional statistical design and analysis techniques. One of the alternative approaches to the construction of a parameter design ex periment utilizes a single experimental design in both the control and noise factors. This approach was introduced by Welch et al. (1990) and discussed further by Shoe maker et al. (1991), Lucas (1991), and Box and Jones (1992). Underlying the single experimental design, now referred to as a "combined array," is a single linear model in which the response is modeled directly as a function of the control factors and the noise factors. Shoemaker et al. (1991) noted that the combined array approach resolves some of the criticisms of Taguchi's crossed array. Specifically, combined ar rays can be more economical in terms of run size or number of runs and by proper construction can allow estimation of all CxC, N xN, and CxN interactions of interest. In this dissertation we aim to study how the combined array approach can be used as an experimental strategy for obtaining a target condition on the mean re sponse while simultaneously minimizing its variance. The combining of the control and noise variables into a common array and treating the noise variables as fixed in the experiment enables us to employ modifications of traditional response surface methodology. Specifically, optimum control factor settings may be determined by constructing response surfaces for the process mean and variance. The response sur face for the process mean is obtained once we fit a single model in the control and noise variables. The response surface for the process variance is obtained by applying the variance operator to the response model in which it is assumed the noise variables are truly random. Then we substitute the coefficient estimates and denote this esti mated response surface for the process variance by f(x). Now the response surface derived from y(x) is used to locate a contour in the experimental region where the target is met. Along this contour, we can then obtain the settings in the control variables that minimize variability. For this purpose we optimize the response surface derived from (x), the estimated process variance. Since methods for exploring the response surface generated by the mean of the response and obtaining reliable choices of control factor settings that optimize (ax) are well documented in the literature, we will instead focus on T(x), the response surface derived from the estimated process variance. This dissertation provides an important contribution to the field of robust parame ter design. Combined arrays are currently being utilized by industrial practitioners for the purpose of designing a quality product or process. This objective is accomplished by determining settings in the control factors that minimize a measure of loss. Since this loss is the process variance, a technique is needed for evaluating which combined array provides "good" estimates of the process variance. However, little consideration has been devoted to the properties of combined arrays. This research represents the first attempt at formalizing the properties of the estimated process variance when a combined array is used. We develop a mean squared error criterion for the estimated process variance. This criterion will be used to compare the ability of commonly used designs to support the estimation of r(x) over a region of interest. Thus, if an experiment is to be run during the product development stage to determine those settings of the control factors that yield a product that performs well in the presence of variability, our mean squared error criterion will enable us to recommend to the practitioner a combined array that minimizes the mean squared error of T(x). The literature review that follows in Chapter 2 is intended to further acquaint the reader with the main aspects of Taguchi's methods for robust parameter design and the criticisms of his methods and to discuss the development of the combined array approach. In addition, since we will utilize a mean squared error criterion to compare designs, some background material regarding mean squared error will be discussed. In Chapter 3 we derive the prediction properties of the combined array and develop a mean squared error criterion for evaluating combined arrays that are used to support a model that is presumed to be adequate. In Chapter 4 we do the same for the case where the fitted model is presumed to be underspecified. In both chapters we also evaluate combined arrays graphically through the use of a mean squared error criterion. In Chapter 5 we examine various special cases of the underspecified model and make recommendations concerning the best combined array to use when protecting against specific higher order models. Chapter 6 contains conclusions and potential areas of further research into the combined array and robust parameter design. CHAPTER 2 REVIEW OF LITERATURE 2.1 A Review of Taguchi's Parameter Design Methods The focus of this dissertation is on the properties of the combined array, an ex perimental design that was developed to be an alternative to Taguchi's crossed array design. However, it is important to review some of Taguchi's methods for the de sign and analysis of robust parameter design experiments in order to understand the particular aspects of these methods that prompted the development of the combined array approach. Kackar (1985) provides a very good review the concept of parameter design and the quality control methods of Genichi Taguchi. Much of the material that follows in this section is summarized from Kackar's article. The evolution of a manufactured product consists of three stages: product design, process design, and manufacturing. It is the designs of both the product and the process that play crucial roles in the resulting quality of the product. A good quality product is one that performs at a target level consistently throughout its life span and over a wide variety of operating conditions. Define performance variation to be the amount by that a manufactured product's performance deviates from its target value under all different operating conditions. Parameter design is an offline quality control technique that is used to make a product robust against sources of variation by identifying settings that minimize performance variation. Parameter design is based upon classifying the variables that affect a product's performance into two categories: control factors and noise factors. The settings of the control factors (design/process parameters) are the product de sign characteristics whose nominal values can be easily controlled or manipulated by the product designer. A vector x' = (x1, ..., xki) of settings of the control variables defines a product design specification. Noise factors are those that affect product performance but are difficult or expensive to control. Noise can be classified into two categories: external and internal. Examples of external sources of noise are ambient temperature and humidity, suppliers of raw materials, and human variations in oper ating the product. Internal sources of noise are those that cause the settings of the control variables in a manufactured product to deviate from the nominal settings. The main sources of internal noise are manufacturing imperfections and product de terioration. Note that not all sources of noise can be included in a parameter design experiment because of physical limitations and because of lack of knowledge about the sources. Let z' = (zi,..., zkn) be the vector of noise variables that are included. Parameter design experiments are possible because the effects of noise factors may change with the settings of the control factors. There could be many settings of a for which the system can perform, on the average, at desired target levels. Among these, there will be some settings that make the system less sensitive to variation in the noise variables z. The basic idea in parameter design is to identify, through exploiting interactions between control variables (C) and noise variables (N), appro priate settings of control variables that make the system's performance robust to uncontrollable variation in z. A parameter design experiment defined by Taguchi consists of crossing two or thogonal arrays: the control (inner) array and the noise (outer) array. Each row or test run of the matrix for the control variables represents a specific product de sign. The rows of the noise matrix represent different combinations of the levels of the noise factors. Note that by choosing extreme values for the levels of the noise variables the system is purposely made to be as noisy as possible. The rationale for this is that settings in the control variables that are determined to perform well in terms of extreme variability should also perform well under normal conditions. Each of the n, rows of the control matrix is crossed with the nz rows of the noise matrix resulting in a crossed array with nxn, rows. For example, suppose there are four con trol variables, each having two levels (low = 1, high = 1), and three noise variables each having two levels (low = 1, high = 1). Attention should be given at this point regarding the coding of variables. Throughout this dissertation it should be noted that for convenience all computations will be made assuming that the factors under study have been transformed from their natural units to design units through the use of appropriately chosen centering and scaling constants. In this case, we choose a control array that is 8 x 4 and represents a 241 fractional factorial design in the control variables. We also choose a noise array that is 4 x 3 and represents a 231 fractional factorial design in the noise variables. Crossing each of the eight rows in the control array with the four rows in the noise array yields a 32run crossed array. This crossed array is illustrated by Figure 2.1. In general, Taguchi proposes the use of saturated or nearly saturated two and threelevel experimental designs for both the control and noise arrays. These designs have been labeled in the literature L4, L8, L9, L12, L16, L18, L20, L24, L27, and L32 OAs. The 4, 8, 9, etc. refer to the number of design runs, and OA stands for orthogonal array. Orthogonal arrays are those n x k matrices that have the following two properties (Rao, 1947): 1. All levels appear an equal number of times in every column. 2. In every pair of columns, all combinations of the levels must occur, and they must occur an equal number of times. The notation A(n, k, s, d) has also been used to denote an orthogonal array where n is the total number of runs, k is the total number of columns, s is the number of levels X1 X2 3 X4 Z1 Z2 Z3 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 Figure 2.1. Illustration of a Crossed Array used, and d is the "strength." An orthogonal array of strength d is an n x k matrix where all sd combinations of levels corresponding to any d factors chosen out of the k occur an equal number of times. All of Taguchi's designs are orthogonal arrays of strength 2. A twolevel design of strength 2 will allow estimation of the main effects of the k factors but not pure higherorder terms or interactions between the factors. Examples of twolevel orthogonal arrays of strength 2 are the L4, L8, L12, L16, L20, L24, and L32 designs. The L4, L8, L16, and L32 designs are the familiar 22, 23, 24, and 2s resolution III fractions of 2k factorial designs. The L12, L20, and L24 designs are twolevel PlackettBurman designs. The L9, L18, and L27 designs are threelevel PlackettBurman designs. Note that these threelevel orthogonal arrays of strength 2 will allow estimation of main effects and pure quadratic terms but not interactions. Lastly, note that Taguchi's experimental designs are formed by crossing the control array with the noise array, and this crossed design is such that all first order and higher C x N interactions are estimable. The crossed array displayed in Figure 2.1 is an L8 control array crossed with an L4 noise array. The analysis of a crossed array uses the nz values of the response variable (or qual ity characteristic) obtained for each of the n, rows of the control matrix to compute n, values of a performance statistic that Taguchi calls a signaltonoise (SN) ratio. A performance statistic is a function of the control variables that estimates the effect of noise factors on the characteristic of interest. The values of this statistic are then used to predict optimal settings of the control variables. The identification of optimal control factor settings requires specification of a criterion to be optimized. One such criterion is the expected loss. Let Y be a value of the performance characteristic of interest and suppose its target value is t. Further, assume that Y is a function of the control and noise variables, i.e., Y = f(x, z). The goal of parameter design is to choose the setting of x that minimizes the expected loss caused by deviations of Y from the target value t. This expected loss is given by R(x) = EzL(Y, t) where L is a loss function. Taguchi recommends that the loss be measured by the quadratic loss function L(Y, t) = K(Y t)2 where K is some constant. In practice, however, Taguchi maximizes SN ratios instead of minimizing R(x) and generally gives no connection between these two optimization problems. Leon et al. (1987) show that for certain underlying models for the response, maximization of the SN ratio does indeed correspond to minimization of the expected quadratic loss. The forms of the SN ratios suggested by Taguchi depend on whether "targetis best," "smallerisbetter," or "largerisbetter." The SN ratio for the targetisbest case takes into account the relationship between the mean and the variance. If the variance increases linearly with the mean, then the appropriate SN ratio is given by SN, = 101oglo(y72/s) where y, and s5 are the sample average and sample variance of the nz response values at the ith control setting, i = 1,2,..., n If the mean and variance are independent, then the appropriate SN ratio is SN, = logo s,. For the cases of smallerisbetter and largerisbetter, the SN ratios are respectively given by n, n,z /y SN= 101oglo(Ey /nz) and SN,=101og0o('). j=1 j=1 nz In each case described above it is desired to maximize the appropriate SN ratio. To carry out this optimization, Taguchi introduces the concept of adjustment factors. These are special control variables whose values can easily be changed to reduce an adjustable component of the total variation about the target without affecting the nonadjustable component of the total variation. In other words the adjustment factors are those control variables that only affect the mean. The remaining control factors affect both the mean and variance or just the variance. The division of control factors into two groups x = (d, a) enables the analyst to minimize the variability using the main control factors d and then adjust the mean to its appropriate target using the adjustment factors a. In summary, Taguchi recommends the following twostep optimization procedure: Step 1: Find the setting d = d* that maximizes the SN ratio. Step 2: Adjust a to a* while d is fixed at d* in order to minimize the bias. Leon et al. (1987) provide the link between maximizing the SN ratio and minimiz ing the loss R(x). They show that if Y = f(x, z) has the multiplicative form given by Y = y(d,a)e(z,d) where E[Y] = p(d,a) is a strictly monotone function of each component of a for each d, then the use of Taguchi's SN ratio in a twostep procedure leads to minimization of quadratic loss. To see this, the authors note that the mul tiplicative model implies that the squared coefficient of variation, var[Y]/E2[Y], does not depend on a (or approximately, for log Y, a affects location but not dispersion). Also, to link the SN ratio with quadratic loss R(x), the authors argue as follows: 1. First note that we can find (d*, a*) such that R(d*, a*) = minda R(d, a). 2. Define the general twostep optimization procedure: Step 1: Find d* that minimizes P(d) = mina R(d, a). Step 2: Find a* that minimizes R(d*, a*). 3. If we use quadratic loss and the multiplicative model, then P(d) is a decreas ing function of Taguchi's SN ratio and thus the above twostep procedure is equivalent to Taguchi's. Taguchi's method for conducting a robust parameter design experiment and sub sequently maximizing a signaltonoise ratio can be summarized as follows: 1. Identify initial and competing settings of the control variables, dividing them into adjustment and nonadjustment variables if possible. Identify the noise variables and specify the ranges of the noise variables within which the product performance is desired to be insensitive. 2. Use orthogonal arrays to construct the control and noise variable matrices. 3. Conduct the parameter design experiment and evaluate the performance statis tic (SN ratio) for each test run of the control matrix. 4. Use the computed values of the SN ratio to predict new settings of the control variables that maximize the SN ratio. 5. Conduct a confirmatory run to check that the new settings do indeed improve the performance statistic. Lastly, note that the foregoing discussion was concerned with the product design stage in the evolution of a manufactured product. In the process design context, parameter design experiments are used to identify settings of controllable process variables at which the effect of manufacturing variation is minimum. The control variables are thus process variables whose operating standard can be chosen by the process designer, and manufacturing variations are sources of noise. 2.2 The Combined Array: An Alternative to Taguchi's Crossed Array Recall that Taguchi's approach uses two experimental designs: one for the con trol factors and one for the noise factors; and for each row in the control array the "replicate" observations generated by the noise array are reduced to a loss statistic or related SN ratio that can be used to compare specific product designs defined by the settings of the control factors. The observed losses are then modeled. A common criticism of this approach is that the model often includes only the main effects of the control variables and little attention is paid to the effect of interactions among the control variables on the underlying quality characteristic Y. Welch et al. (1990) proposed a new approach to parameter design that models a quality characteristic Y(x, z) as a function, typically polynomial, of both x and z. The single experimental design for both types of parameters generally requires far fewer runs than crossed arrays. The method of Welch et al. (1990) involves the following six steps: 1. Design a single experiment to predict the response of interest as a function of the control factors x and the noise factors z, where the response model is given by The fjs are assumed known functions of x and z and the /3js are unknown model coefficients to be estimated from the data. Hence the combined array approach selects designs on the basis of their ability to estimate the appropri ate coefficients in the presumed model for the response. Attention should be given at this point regarding the application in this article. The authors used computer simulation to model a clock skew circuit design and thus the error term above represents systematic departure from the assumed model because the computer response Y is deterministic. 2. Predict Y(x, z) using the ordinary least squares estimator Y(x,z) = E f(x, z). 3. Once the terms in the model are appropriately estimated, the experimenter can predict the value of the loss statistic R(x) for a given x. For the example in which a clock skew circuit design is modeled using computer simulation, the target for the simulation output Y(x, z) is zero, and following Taguchi, the expected squared error loss is given by R() = JY2(mx,z)g(z)dz, which averages a measure of loss over the distribution g(z) of the noise factors. Denote the predicted value of this loss statistic by R(x) fY2(x, z)g(z)dz. Welch et al. (1990), for the purposes of their example, assume a uniform distri bution for the noise factors. 4. Minimize the estimated expected loss !R(x) as a function of x. 5. Conduct a confirmatory experiment to evaluate R(x) by fixing x at the values) found in step four and varying z according to g(z). 6. Iterate if necessary. The optimization step may suggest for example that the design region needs to be shifted. Shoemaker et al. (1991) suggest complementing the formal analysis of the expected loss in step four above with graphical dataanalytic techniques. Specifically, the authors recommend the use of controlbynoise interaction plots to help determine the control factor settings that dampen the effects of the individual noise factors. To illustrate and compare their method with others, Welch et al. (1990) pro vided an example in which the goal was to determine the widths wl,..., w6 for the transistors that gave the smallest clock skews (si and s2) in the presence of process variability. Three strategies were employed. The first was a 60run experiment for the skew models from which two specified loss statistics were predicted indirectly (this is the new method of Welch et al., 1990). The second was a (40 x 5)run crossed array (i.e., separate control and noise array as advocated by Taguchi) for modeling the loss statistics directly. The third was a hybrid strategy that used the (40 x 5) experiment from the second strategy but modeled the skews to predict the losses indirectly as in the first strategy. A comparison of the predicted values of the loss statistics with the actual values obtained in confirmatory runs for the best three vectors of control factor settings resulted in the following conclusions: 1. Strategy one predicted the skews accurately enough to give predictions very close to the actual losses determined by confirmatory experiments. The pre dictions were more reliable and the actual losses were smaller than in strategy two. 2. Strategy one gave virtually identical results as strategy three, indicating that small experiments can be adequate. 3. Modeling the skews (i.e., the response) appeared to be more accurate than modeling the loss statistics directly. 4. Strategy one gave the same bestthree circuit designs for the two loss statistics considered. 5. The authors point out that choosing models for complex loss functions, as in the second strategy, is difficult and often not very intuitive. When the response was modeled directly, certain control factor interactions were justifiably omitted, but the same interactions may not be negligible when modeling the loss statistics. A full secondorder model in the widths was chosen for both loss statistics. Also, modeling the response generally admits simpler models. Furthermore, collapsing the data to a loss statistic could hide important relationships in the data. The method of Welch et al. (1990) can be summarized simply as follows: combine control and noise factors into a single array; model the response rather than expected loss; and approximate a prediction model for the loss based on the fitted response model. The single array for the control and noise factors has since become known as the combined array. 2.3 Run Savings and Estimation Flexibility of the Combined Array Shoemaker et al. (1991) further developed the responsemodel/combinedarray approach suggested by Welch et al. (1990). The authors' contributions to parameter design are motivated by what they consider to be the following disadvantages of Taguchi's loss model approach and corresponding crossed array designs: 1. A very large number of runs may be required because the noise array (NA) is repeated for every row in the control array (CA). 2. Many degrees of freedom (df) are dedicated to estimation of all interactions between control factors and noise factors. 3. The focus is on modeling the loss, or related SN ratio, which is often a nonlinear, manytoone transformation of the response (i.e., performance characteristic) Y. Even when Y follows a linear model in the control and noise factors, it is less likely that the loss can be modeled well by a low order linear model even if data transformations are employed. 4. Modeling the loss may also hide some of the relationships between individual control and noise factors. To elaborate on the second disadvantage stated above, Shoemaker et al. (1991) show that for the general crossed array case, the effects that are estimable in a CA x NA design correspond to those that are estimable in CA, those that are estimable in NA, and the generalized interactions between the estimable effects in the control array and the estimable effects in the noise array. This fact illustrates one of the limitations of the crossed array approach. If the control array has n., runs and the noise array has n. runs, then the CA x NA design has nnz runs, and of the nnz 1 df in CA x NA, (n, 1)(nz 1) df are used for estimating twofactor and higher order controlbynoise (C x N) interactions. In practice, however, some of the C x N interactions may be judged a priori to be unimportant; so, the large run size of the crossed array is often unnecessary. Also, there is no flexibility to reallocate some of the df to estimate controlbycontrol (C x C) interactions that may be important. Taguchi's threelevel designs are also unable to estimate (C x C) interactions. The inflexibility and large run size of the crossed array can be avoided by combin ing the control and noise factors in a single array. This allows for a smaller experiment since a more highly fractionated design may be chosen. Also, with suitable choice of a defining relation, a fractional factorial design can be constructed in which some C x C interactions will be estimable if certain C x N interactions are assumed neg ligible. To illustrate this, Shoemaker et al. (1991) provide an example to compare the estimation flexibility of a crossed array design with a combined array design. Three control factors and three noise factors, each with two levels, were used. To construct the crossed array, a 2R'1 CA was crossed with a 231' NA. The resulting 16 run, six variable crossed array is equivalent to a 262 fractional factorial design. On the other hand, by combining the six factors into a single array, a 2'62 fractional factorial design can be constructed using defining relations that allow estimation of certain C x C interactions. From the example, it was observed that the combined array is superior to the crossed array because it allows more effects to be estimated under weaker assumptions. In conclusion, Shoemaker et al. (1991) illustrate that the run savings allowed by the combined array format comes from the flexibility to estimate those effects that are most likely to be important, rather than forced estimation of a large number of controlbynoise interactions. 2.4 Analysis of a Combined Array Experiment The development of optimization strategies for robust parameter design appli cations is important for the statistician in industrial practice. Recall that robust parameter design experiments are possible because the effects of noise factors may change with the settings of the control factors. Since the variance fluctuates over the region of interest, Taguchi was the first to emphasize the possibility of modeling the variance. The basic idea in parameter design is to identify, through exploiting interactions between control variables (C) and noise variables (N), appropriate set tings of control variables that make the product or process performance robust to uncontrollable variation in z. In other words, although there may be many settings of x for which the product or process can perform, on the average, at desired target levels, the goal is to determine the setting that makes the product or process less sensitive to variation in the noise variables z, and this is done by exploiting C x N interactions. Myers, Khuri, and Vining (1992) show that if the model for the response contains C x N interactions, then the process variance is a function of the levels of the con trol variables; on the other hand, if C x N interactions are not included (i.e., they are deemed to be unnecessary or insignificant), then the process variance does not depend on x, and in this case the "robust" parameter design problem is nonexistent. Consequently, when interactions between control and noise variables are found, the process variability depends on the specific values of the control variables through the C x N interactions included in the underlying model. The interaction structure then can provide a reasonable basis for estimating the process variance where the noise factors are truly random variables. As an example, suppose that for a particular combined array experiment, there is one control variable and one noise variable. If there is no interaction between the control and noise variables, then the fixed effects model for the experiment is given by y= o + xi 71z +yii c. Now, in reality, z1 is a random variable with a variance ao, and thus the process variance is given by var(y) 1= 72( + C2. Notice here that the variance does not depend on the design levels in x\. Hence there is no need to incorporate the process variance as a part of the performance criterion for parameter design. On the other hand, consider the fixed effects model for the combined array experiment where it is assumed that there is interaction between the control variable and noise variable. This model is given by y = /3o + iXi + 71zi + Aixxzi + e. (2.1) Again assume that in reality the levels of the noise variable are random in the process so the process variance is given by var(y) = (71 + Aixi)2 a + a (2.2) The process variance is now a function of the levels of the control variable, and thus it is important that the process variance be involved in the response surface performance criterion. The example above has provided motivation for including the process variance in the choice of optimum settings for the control variables. Note that the variance function (2.2) is essentially a quadratic response surface. The underlying model (2.1), where the noise variables are treated as fixed effects, provides a basis for estimating '71 and A,, and thus allows estimation of the process variance. The estimated variance response surface is given by var(y)= (71+ +iXx)2 + ^2 where we assume that or is known. Now, assume that E(e) = 0, and without loss of generality assume that E(zi) = 0. Then E(y) = 3o + A/1x, and so the estimated response surface for the mean is given by (x) = 3o + /xi. As a result, optimum conditions may be obtained through simultaneous exploration of these two response surfaces. One method for exploring the response surfaces is through the construction of contour plots. Joint consideration of these contour plots can help determine the optimum conditions where the target value for the response is obtained with minimum variability. More formal methods for dual response opti mization can be found in Myers and Carter (1973), Vining and Myers (1990), and Del Castillo and Montgomery (1993). 2.5 The Mean Squared Error Criterion Let us suppose that the true functional relationship between a response r7 and k continuous variables x1, x2,... Xk is given by yf = g(x, 0). When the form of the true functional relationship is unknown, the objective is to approximate within a given region R of the kdimensional space of the variables, the function g(x, 0) by some function f(a, /3). Since the fitted model may be an inadequate representation (i.e., an underestimate) of the true response surface, there is a need to consider two sources of error. Error can arise both from sampling variance and bias introduced by model inadequacy. In other words, we explain the discrepancy between the true function r(x) and the fitted function y(x) as being due to "variance error" and "bias error." In this problem setting, Box and Draper (1959, 1963) propose that the design D, i.e., an n x k matrix whose (ui)th element, xu,, is the value of the ith input variable used in the uth experimental run (u = 1, 2,..., n; i = 1,2,... k), should be chosen so as to minimize the quantity nK J JR E[y(x) r(x)]2dx where dx = dxedx2 ... dxk, K = [fR dx]', E[ (x) i(x)]2 is the mean squared error of the predictor y(x), o2 is the experimental error variance, and n is the number of design runs. Here we interpret J as being the mean squared deviation of y(x) from the true response 7r(x), averaged over the region R and normalized with respect to the number of designs runs and the experimental error variance. It can easily be shown that J separates into contributions from the "variance error" and the "bias error" according to J=nK var[ (x)]dx + ))2d= V + B. This sum of the average variance and the average squared bias over the region R has become known as the BoxDraper "JCriterion" or the integrated mean squared error (IMSE) criterion. The utility of the JCriterion is that it represents a compromise between purely variance based or purely bias based criteria that had dominated the literature on the construction of optimal designs prior to Box and Draper's work. Box and Draper (1959, 1963) apply the JCriterion to the setting in which f(x, 3) is assumed to be a polynomial of degree d, while the true functional form g(x, 0) is a polynomial of degree d2 = d1 + 1 where di was taken to be equal to 1 and 2. Designs having minimum integrated mean squared error were obtained for various ratios of V/B and compared with "all variance" (B = 0) and "all bias" (V = 0) designs. The designs were constructed to satisfy certain moment conditions derived from the minimization of J. It was concluded that if the practitioner is concerned with protecting against model inadequacy, then the design that minimizes J is close to the "all bias" design. In the combined array approach, recall that we assume a model for the response and then approximate a prediction model for the loss based on the fitted response model. In practical circumstances, it is reasonable to believe that the assumed model is always to some extent incorrect. We will develop a JCriterion for the estimated process variance, in light of Box and Draper's work, not for constructing a new class of optimal designs, but rather as a criterion for choosing between already existing designs that have been used for constructing combined arrays. Again we emphasize that in the robust parameter design problem that we are studying, the focus is not on minimizing the bias and variance of the predicted response y(x). Rather, the focus is on the estimated process variance, var(y), which we denote by (ax). The rationale behind this is that we assume that a set of operating conditions for which the mean response achieves a target value has been found, and the objective is to determine from this set the setting x that achieves the target value with minimum variability. These conditions are such that (a) is minimized. Hence, we focus on the properties of N(a). It is hoped that this review has set the stage for the following chapters in which the reader will be shown how the abundance of literature in traditional response surface methodology can be utilized and extended to develop properties of the combined array. Once these properties are established, the theory will be laid for the industrial practitioner who must choose the most appropriate experimental design to run for his or her particular robust parameter design problem. CHAPTER 3 ANALYSIS OF THE COMBINED ARRAY WHEN THE ASSUMED MODEL IS CORRECT 3.1 Introduction In the combined array approach designs are selected on the basis of their abil ity to support the estimation of the coefficients in the presumed model for the re sponse. Once these terms are estimated, the analyst is able to estimate the appropri ate squared error loss. Typically, an additive error structure is assumed, in which case the process variance is the appropriate loss function. In a robust parameter design problem the predicted response and the predicted variability are equally important for determining the optimum control factor settings. We assume the experiment will be run in a region where the predicted response is likely to achieve a target value, and our goal is to determine where in the region the variability around the target value is minimized. Hence the estimated process variance becomes the performance criterion that we focus on. In this chapter, we will determine the properties of the estimated process variance under the assumption that the fitted model is correct. The statistician would like to choose an experimental design that will produce accurate and reliable estimates of the process variance over the entire experimental region. Hence, we will develop a mean squared error criterion in order to evaluate the properties of the estimated process variance obtained from a given design. For a given model, a series of sim ulations will be conducted to illustrate the use of the mean squared error criterion for comparing commonly used designs to fit the given model. An additional practical example outlines how this approach can be used to select an optimal combined array among a given set of experimental designs. The mean squared errors will be displayed graphically in order to make conclusions and recommendations regarding competitor designs. 3.2 Properties of the Estimated Process Variance For a combined array experiment suppose that the model for the process is given by y = 3o + x', + z'y + x'Az + e (3.1) where /o is the intercept, x is a p x 1 vector representing the appropriate polynomial expansion of the control variables for the assumed model, 3 is the p x 1 vector of unknown parameters associated with the control factors, z is a q x 1 vector representing the appropriate polynomial expansion of the noise variables for the assumed model, 7 is the q x 1 vector of unknown parameters associated with the noise factors, A is the p x q matrix of unknown parameters associated with the controlby noise interactions (some elements of A may be zero), and e is the random error. We shall assume that the es are independent and identically distributed with mean zero and variance ucy. For the purposes of the particular planned experiment, we shall treat the noise variables as having fixed effects, and thus there is no distinction between the roles of the control and noise variables. In the process, however, the noise factors are truly random variables with E(z) = p and var(z) = V, where we shall assume V is known. Furthermore, the error terms and noise variables shall be assumed to be independent. Assume the appropriate combined array experiment has been run to estimate the linear model given by equation (3.1). The analysis focuses on both the response and an appropriate loss function. Leon et al. (1987) show that if the goal of experimentation is to find conditions for which a target value is achieved with minimum variability and if an additive model as in (3.1) is assumed, then the appropriate loss function is the process variance, which we shall denote by r(x). Since the noise variables are considered to be random in the process, the value of the process variance at a point x inside a region of interest R is, for model (3.1), given by T(x) = var(y) = var[(Y' + x'A)z + e] = [7' + x'A]V [7' + x'A]' + a2 where V is the q x q variancecovariance matrix for z. Now if we define O(x) = 7 + A'a, we can write r(x) = var(y) = 0'(x)V (ax) + a~l. (3.2) Since O(x) is unknown, we use the results of the combined array experiment to estimate O(x) with 4(x). The estimated process variance is then given by f(x) = ' (x)V<>(xs) + <2.. The aim is to find control factor settings x that minimize (ax), subject to the constraint that the predicted response is at the target value. The derivation of prediction properties of f(x) will make use of the following more convenient notation. Let A' = [y, A'] and x'a = [1,x']. Then O(x) may be reexpressed as O(x) = A'xa. Also for convenience, let A = vec(Aa), where the vector of the matrix Aa is the q(p + 1) x 1 column vector formed by stacking the columns of Aa one under the other. Using this notation we can rewrite the model given by (3.1) in the more general form y = P +E where y is the n x 1 vector of responses, P is the n x m model matrix (m = (p + 1)(q + 1)) whose elements are known functions of the settings of k = kc + k, input variables determined by an n x k design matrix, i.e., P = [1, X, Z, XZ] where X and Z are n x p and n x q design matrices in the control and noise variables, respectively, and XZ is an n x pq matrix of cross products, P' = [/o, /3', A'] is the m x 1 vector of unknown parameters, and e is the n x 1 random error vector. If we assume that E[E] = 0 and var[e] = olI and use ordinary least squares (OLS) to estimate the model coefficients, the resulting estimate of the coefficient vector is , = (P'P)1'Py. By the properties of OLS estimation and under the assumptions that the model is correct and the noise variables have fixed effects in the experiment, E[ ] = 4 and var[ ] = cy(P'P)'. Note that /(x) = ,+A 'a is a linear combination of the elements of ,, and so E[ (x)] = O(). Thus, E[(x)] = E[k'(x)V(x) + ] = tr(VE) + 0'(x)V((x) + cU (3.3) where E = var[4(x)]. From (3.2) and (3.3) we note that bias[f(xa)] = E[f(ax)] r(x) = tr(VE). (3.4) Note that V and E are both variancecovariance matrices. Consequently, V and E are positive definite. Thus, it follows that tr(VE) > 0 and so bias[(xa)] > 0 (Graybill, 1983, p. 397). We thus obtain an unexpected result, which is that ?(ax) is a biased estimator of r(x) even though the presumed model is believed to be correct. Next if we assume the random error vector e is distributed as multivariate normal with mean 0 and variance a2I then the variance of f(xa) is given by 22 var[(a)] = var[ck (x)Vo(x) + &,] = 2tr(VE)2 + 40'(x)VYEV(ax) + 2u (3.5) n m where n is the total number of design runs and m is the number of unknown coefficients in the model. 3.3 Expressions for the Bias and Variance of the Estimated Process Variance To derive a more explicit expression for the bias, assume the correct model is given by (3.1). From least squares estimation of the equivalent form of the model given by y = Pi + e, we noted above that ^ = [#o, 3', A']' = (P'P)lP'y with var[^] = a2(P'P)'. Let C be the appropriate q(p + 1) x q(p + 1) submatrix of a7(P'P)1 corresponding to the variancecovariance structure of A. Note that C is symmetric and can be partitioned into q2 blocks of (p + 1) x (p + 1) matrices of the form C11 C12 CIq' q C 2 C12 C22 C2q (3.6) C.q 2q C qq Since (ac) = A'am, and A = vec(Aa) it follows that (ax) = (Iq 0 X'a)A, where denotes the Kronecker cross product. Hence E = var[(B(x)] = (Iq 0 X')C(Iq 0 )' = (Iq 0 X'a)C(Iq 0 a), where we used the general result that (A B)' = A' 0 B'. It is easy to see that S {aYxaCijXa} for i,j = 1,2,...,q and let V = {ciij} for i,j,= 1,2,...,q. Now recall that for two q x q matrices A and B, tr(AB) = Zf Z. aijbji. Thus, it follows that q q bias[(.(x)] = tr(VE) = a 2 orijsXaCjiXa i=1 j=l q q1 q T 7, u ( Cjx.+ 2 2 ^ C3 j a (3.7) i=1 i where we used the fact that aij = Oji and Cj = C'.. Further simplification of the bias expression depends upon the form of the Cij matrices which in turn depends on the model matrix P for a given design. Next, consider the variance of f(x). Utilizing the result that O(x) = (Iq 0 x')A, equation (3.5) can be expressed as 274 var[( (x)] = 2tr(VE)2 + 4A'[Iq 0 Xa]VEV[Iq 0 x'a]A + 2 (3.8) Our goal is to evaluate var[f(x)] for any given xa, but observe that var[(x)] is a function not only of the particular location in the design region but the unknown parameter vector A as well, which presents obvious difficulties. Our approaches for overcoming these difficulties will make extensive use of the Euclidean norm of the vector A, which we shall denote by r = JAI. We shall assume that we have some reasonable idea about the value of K. Next we use a result proven by Stroud (1971). This result states that if U is the surface of the unit sphere and ds8 a differential on this surface, then the average value of the quadratic form A'AA over all possible A's with length K is given by K2 fu A'AAdsA (A dtr(A) fu ds \ nN where n\ is the number of nonzero elements of A. Applying this result to equation (3.8), the average value for var[,a(x)], denoted by varA[f (X)], over all possible As with length n, is )] fu var[ (x)]dsA f dvarA = 2tr(VVE)2 + 4 tr[(Iq 0 a)V V(I, 0 X')] + 2u4 (3.9) nx n m where nA is the number of elements in A not presumed a prior to be 0. From (3.4) and (3.9) we can obtain MSE[f(x)] = [bias(f(x))]2 + varA[(X)] 4c2 = [tr(VE)]2 + 2tr(VE)2 + tr[(Iq xa)VEV(Iq 0 x')] nA 2,4 +  (3.10) which can be readily calculated at any point Xa. 3.4 Comparing Designs Based on the Mean Squared Error Criterion A reasonable method for comparing designs is to evaluate the average or integrated mean squared error (IMSE) of f(x) for each design. The IMSE is given by IMSE[f(a)] = K I MSE[f(x)]dx = K J[tr(VE)]2ddx f 4r2 2U4 + K [2tr(V )2 + tr[(q )Vx V(Iq, 0 ')] + ]d (3.11) R nA n m where R is the appropriate region of interest and K = [fR dx]1. Thus, analogous to the BoxDraper (1959) J = B + V criterion, we note that the average mean squared error of (xa) is the sum of the average squared bias of (ax) and the average variance of (x), where average means averaged over the region R. The IMSE can be found using Monte Carlo integration techniques. A series of simulations will be conducted to illustrate the proposed IMSE criterion for comparing designs. We will consider a robust parameter design experiment con sisting of three control variables and two noise variables. Suppose we use the model given by (3.1), and assume the model is first order in the control and noise variables and contains first order C x N interactions. This model, in expanded form, is given by y = /o + /lXl + 4/2x2 + /3X3 + 71'YZl + 7'22 + Alx1lz1 + A21X2Z1 + A31X3Z1 + A12xiZ2 + A22X2Z2 + A32X3z2 + e. (3.12) From equation (3.11), it can be noted that computation of the integrated mean squared error of a(x) at any point x depends on V, the assumed known variance covariance matrix for the noise variables, and the value of K2. For this example with two noise variables, let orn, o12, and 022 denote the elements of the V. Since our method assumes the noise variables follow a normal distribution, it seems reason able to estimate all and a22 using (Rfi/4)2, where Ri denotes the estimated range (specified by the experimenter) of the it noise variable, i = 1, 2. If the noise vari ables are each studied at two levels in the experiment, conveniently denoted by 1 and +1 in coded (i.e., centered and scaled) design units, then 1, = 2, i = 1, 2, and aO' = 622 = (2/4)2 = 0.25. (Recall that throughout this dissertation all computa tions are made assuming that the factors under study have been transformed from their natural units to design units through the use of appropriately chosen center ing and scaling constants. We prove in Appendix A that the mean squared error of (x) is invariant to our choice of scaling in the noise variables.) Next, since 012 = cov(zl,z2) = p12(0o11622)2 = P12/4, where p12 is Pearson's coefficient of corre lation, we have transformed the problem so that the experimenter now only needs to have knowledge of the correlation structure of the noise variables. In other words, the matrix V can be expressed as a constant times the correlation matrix: V = 1 1 p12 4 p12 1 J To investigate the effect of the correlation structure of the noise factors on the bias and mean squared error, we shall consider a sequence of values P12 between 1 and 1. Next, note that if iR? is misspecified, then it can easily be shown that var[zi] = l(Ri/Ri)2. Hence, to examine the impact of misspecifying the range on the bias and mean squared error, we shall consider three cases: 1. 2. II = 1; 722 = .25 < (Ri = R1/2,R2 = R2); 3. oai = 2.25; 022 = .25 = (Ri = R1/3, R2 = R2) In cases two and three the corresponding V matrices are respectively given by V = 1 [2 P12 and V= [3 3 P12 2 p12 1/2 and 4 p12 1/3 Note that we need not consider the possibility of misspecifying the range of both noise variables by the same constant k. Clearly, if R1 = kR1 and R2 = kR2 then V = var[z] will be multiplied by k and the resulting mean squared error in (3.10) will be multiplied by k2. Thus, although the magnitude of the mean squared error will change, the relative ranking of the competitor designs with respect to the mean squared error will not change. Corresponding to each of the above three cases, we shall consider two values of K2 = 1Al2. In this example A' = (71, A, A21, A31i, 72, A12, A22, A32) where the elements of A are the coefficients of the coded noise factors and controlbynoise interactions. Without loss of generality assume that a, = 1. Now if we assume that the magnitude of each element of A is four times the magnitude of ao, then n2 = na(16au) = (8)(16oa) = 128. Secondly, if we assume that the magnitude of each element of A is twelve times the magnitude of ao, then c2 = (8)(144U2) = 1152. Consequently, the constant 402/nx in equation (3.11) assumes the values 64 and 576. Thus, the three V matrices and two values of K2 represent a total of six different cases of interest. These six cases intentionally cover a wide variety of states of nature that will allow for a fair comparison of different combined arrays. In addition, this will enable us to investigate whether the methodology developed in this paper for comparing designs is "robust" to the various conditions that the practitioner may come across. To illustrate the use of the IMSE criterion, we will compare four combined arrays that support the estimation of model (3.12): Design 1: 16run Resolution V Fractional Factorial Design Design 2: 20run PlackettBurman Design Design 3: 24run PlackettBurman Design Design 4: 24run Orthogonal MainEffect Plan. These designs are easy to construct, and for strict firstorder models they are or thogonal. Since model (3.12) contains interactions between the control and noise variables, designs 24 are nonorthogonal, but they still support the estimation of all of the parameters in the model. 3.4.1 Design 1: 16run Resolution V Fractional Factorial Design A full 2k factorial design, expressed in design variables, is every possible combina tion of 1's for the k factors of interest. The full 2k factorial design allows estimation of all main effects and first order and higher interactions among the variables. How ever, when fitting a model such as (3.12), we are not interested in estimating higher order interactions; so a fraction of a 25 factorial design may be used. A Resolution V fraction of a factorial design is an orthogonal array where main effects are aliased with fourfactor interactions, and first order interactions are aliased with threefactor interactions. By "aliased" we mean that the main effects cannot be estimated inde pendently of the fourfactor interactions, and the first order interactions cannot be estimated independently of the threefactor interactions. However, if model (3.12) is correct, these higher order interactions are assumed to be negligible. For the five variable example that we are considering, the 251 Resolution V design matrix is illustrated in Figure 3.1. x1 X2 X3 zl z2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 Figure 3.1. 16run Resolution V Fractional Factorial Design Computation of the mean squared error of r(ax) given by equation (3.10) is sim plified for the case when the combined array is a Resolution V or higher, twolevel fractional factorial design. For such a fractional factorial design that is appropriately constructed to estimate the C x N interactions in the model, the matrix C given by equation (3.6) has elements Cij = 0 for i 7 j and Ci = 4I. Consequently, the bias of T(s) given by equation (3.7) reduces to 2 q bias [f(a)] = ' E cr n i=1 where r2 = (x'xa,) = 1 + EZ1 xz. Thus, note that the bias depends the variance of the zi's, but not on the correlation between zi and zj. In other words, the bias is constant over values of pij. Under case 1 in the example above, where p = kc = 3, q = kn = 2, orn = var(zi) = .25, 0C22 = var(2) = .25, n = 16, and o2 = 1, we have bias[f(x)] = 6(.25 + .25)= . Since Xa = (1, x1,x2, X3)', it follows that r2 = 1 + X2 + x2 + x2. Thus, the average bias can be determined as follows: K bias[V(a)]dx = K2 [(1 + x + x2 + x)d S1l f1l f11(1 + + x2 + x)dxidx2dx3 32 f l fi fdxldx2dx3 24 1 2 ) = = .0625. 32(23) 16 Similarly, it can be shown that the average squared bias portion of the IMSE is equal to 0.0041667. Next consider the maximum value of bias[a(x)]. In general, since we assume the control factors have been coded such that 1 xi : 1 i = 1,..., k, the maximum value of r2 is kc + 1, and so the maximum bias is S(kc + 1) q n i=1 For the example being discussed it is easy to show that the maximum bias is equal to 0.125. It is not surprising to note that the bias is maximized at the perimeter of the region of interest. Also note that in general, the bias of a(x) can be large if the design is highly fractionated, if the inherent process variability o2 is large, or if there is substantial variability in the noise variables. Next consider the variance of "(x). For the Resolution V fractional factorial design, recall that the matrix C given by equation (3.6) simplifies to 2I, and thus E = var[4(x)] = (Iq 0 X')C(Iq 0 Xa) 2 n 2 2 2 X1 X I n n Thus varA[T(x)] given by equation (3.9) reduces to 294r4 4 + 2 2 2 4 vaA[r(x)] = r tr(V2) + 4 tr[(I, x)V2(I, )] +  n nA na nm 24r 442 2 2 224 = 2ortr(V2) + 4r, tr [V2 0 XaX'a + n nA n n m S2o4r4 V2) + 4K2 o2r2 (V2,)( + 2o74 = ( + )r tr(V2) + ( n n nn + ( m) 2 4 24 n n, n i=1j=1 (n m) where we used the general results that tr(xa'aX) = r2, tr(A 0 B) = tr(A)tr(B), and tr(V2) a= =1 ,=J oj for a q x q symmetric matrix V having elements aij, i,j = 1,2,..., q (Graybill, 1983, p. 301). Under case 1 in the example above, where p = 3, q = 2, o = 1, n = 16, 4K2/n = 64, and m = 12 and if V=[ .25 .25 S.25 .25 ' then 2 64 2 129 varA[f(X)] 2 + 4)r4(.0625 + .0625 + .0625 + .0625) + = 12r + .5. Recalling that r4 (4 X a)2 = (1 + x1 + X2 + X2)2, we can now determine the average variance portion of the IMSE as follows: K f1 11l ft1,j [( + X2 + X2 + x2)2 + .5]dxdx2d3 K 4 varA [ a()]dx = . l 1 1298"1 1 2 3 + A JR f1f1 fl 1, dxedx2dx3 129 = 129 (4.2666667) + .5 = 4.775. 128 The computations above are meant to illustrate that for this 16run Resolution V design, the IMSE's can be analytically computed. However, these analytic compu tations are not feasible for the other designs being considered so we will numerically compute the IMSE's using Monte Carlo integration. The results will be displayed graphically in Section 3.5. 3.4.2 Design 2 and Design 3: 20run and 24run PlackettBurman Designs Plackett and Burman (1946) demonstrate how to construct Resolution III orthog onal arrays involving 4L design runs, L = 2, 3,..., for up to 4L 1 variables. For the fivevariable example, the design matrices for the second and third combined array experiments were obtained from the first five columns of PlackettBurman designs for n=20 and n=24, respectively. These designs are displayed in Figures 3.2 and 3.3 respectively. Note that we did not investigate whether there exists a better choice of columns from the 20 and 24run PlackettBurman designs in their entirety. For X1 X2 X3 Z1 z2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 Figure 3.2. 20run PlackettBurman Design Xl X2 X3 Z3 Z2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 Figure 3.3. 24run Plackett Burman Design the purposes of this dissertation, a design was chosen only for its ability to estimate the model. The question of whether a better set of columns exists for estimating the models of interest (Draper and Lin, 1990) will be left for future research. For these designs, the P'P matrices are invertible but not diagonal and thus the Cij matri ces are not 0, making it impractical to reduce equations (3.7) and (3.9) into simpler forms. 3.4.3 Design 4: 24run Orthogonal MainEffect Plan Orthogonal maineffect plans represent a class of designs for planning industrial experiments with mixed levels. Methods for construction of these designs can be found in many articles including Wang and Wu (1991), Cheng (1989), and Nigam and Gupta (1985). The 24run orthogonal array design that we used was obtained from columns 2, 4, 6, 18, and 20 of the matrix [B B2 B2] 3B1 B2 B2 where B = [B1 : B2] is a Hadamard matrix of order twelve (H12) with the first column of ones removed, B1 is the second column of H12, and B2 represents the remaining 10 columns of H12. The design matrix is illustrated in Figure 3.4. 3.5 Results and Discussion For each of the four designs and six states of nature corresponding to the example discussed above, the IMSEs will be calculated using Monte Carlo integration with points generated over the cuboidal region R, defined in coded variables by R = {x = (x, x2, x3)': 1 < xi < 1, i = 1,2,3}. Since the designs employ different numbers of design runs, the IMSEs shall be displayed first unweighted and then weighted X1 X2 X3 Z1 Z2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 Figure 3.4. 24Run Orthogonal Array (multiplied) by the number of runs in order to provide fair comparisons across the designs. Figure 3.5 considers the bias contribution to the mean squared error for cases 13 discussed above. (Cases 46 need not be considered since the bias contribution is independent of the value of the constant K2 in equation (3.10).) The left hand column displays the unweighted average bias while the right hand column displays the weighted average bias. The graphs in Figure 3.5 illustrate the effect of correlation on the average bias for the example and designs discussed above. By comparing the vertical axes of the graphs, it can be noted that the average bias increases as the value of oar = var(zl) increases, although the relative ranking of design performance is not affected by the value of (11. Note that for the Resolution V design and the 24run PlackettBurman design, the average bias is constant over the values of p12 between 1 and 1. Also, when multiplied by the number of runs, the average bias for the Resolution V design is smallest. With regard to the PlackettBurman designs, the 24run has the smallest average bias when unweighted by the number of runs, but after weighting, there is little difference between the performance of the 20 and 24run designs across the three cases. Interestingly, when weighted, the 20run is worse than the 24run design only for large positive correlation. The orthogonal maineffect plan had the largest average bias for each case, regardless of weighting. One interesting feature of the orthogonal maineffect plan however is that the average bias decreases as the P12 goes from 1 to 1. This suggests the possibility of the existence of a nonorthogonal design with smaller average bias than an orthogonal design. Figures 3.6 and 3.7 illustrate the effect of correlation on the average MSE. In 3.6 where the average MSE is unweighted by the respective numbers of design runs, observe that the 24run PlackettBurman design performs best (i.e., has the smallest IMSEs). Interestingly, in the first and second plots in the left hand column, observe that the value of aon has an impact on the performance of the Resolution V design 0 = a2 = 0.25 1.0 0.5 0.0 0.5 1.0 P 12 a11 =1.0; 022=0.25 1.0 .5 0.0 05 1.0 Pl2 S,,= 2.25; oa =0.25 1.0 0.5 0.0 0,5 1.0 P12 ,11 = 1.0; 22= 0.25 1,0 05 0.0 05 10 P 12 G11=2.25; a22 =0.25 1 1.0 5 0.0 0,5 LO 1.0 .5 0.0 05 1O Legend: RESV .......... PB20 PB24 __ OA24 Figure 3.5. The Effect of Correlation on the Average Bias (Unweighted Weighted) o = = 0.25 <0 S,,=O 22= 0.25; K = 1152 10 05 00 05 1.0 o, =1.0;o = 0.25; K= 128 \ / 10 05 00 0,5 O P12 T = 2.25; a2= 0.25; K2 = 128 0 05 00 / N \ / ., ..  o' ^ . \ / 10 OS 00 05 10 P 12 a = 2.250 =0.25; c 2= 1152 10 5 00 0.5 0/ 10 05 0,0 05 10 Legend: RESV ............ PB20  PB 24 OA24 Figure 3.6. The Effect of Correlation on the Average MSE o, = Y = 0.25; Kx = 128 an= a22= 0.25; K= 128 \ \ / S\ / \ / ", N / .., "',. ___ ,7 10 0.5 00 0.5 1.0 P 12 ao = 1.0; 22 =0.25; 2= 128 \ \ / \/ \ / .1.0 0.5 00 0.5 1.0 P12 N\ \ / \ / / "'. , .. .. " 10 0.5 00 05 1.O0 ao = 22= 0.25; 2= 1152 .\ / S12 1.0 0.5 0,O 0.5 10 P 12 ao =1.0;o5 = 0.25; 12= 1152 0\ / \ \ / 7 .1,0 0.5 0.0 0.5 10 P12 oa=2.25;o22 =0.25; K2 1152 \ / 1.0 0,5 0.0 0.5 1.0 Legend: RESV ............ PB20 PB24 _ OA24 Figure 3.7. The Effect of Correlation on the Weighted Average MSE relative to the 24run orthogonal array. Also, a comparison of the first plot on each of the left and right side reveals that the value of K2 has an impact on the performance of the Resolution V design relative to the 24run orthogonal array and the 20run PlackettBurman design. Otherwise, the remaining plots show that the values of an and K2 have little impact on relative design performance. Lastly, note that the average MSE increases as the value of p12 deviates from zero in either direction. Overall, the correlation structure of the noise variables has an impact in an absolute sense but only a small impact in a relative sense. In Figure 3.7, where the average MSEs are multiplied by the respective numbers of design runs, observe that the Resolution V design now performs uniformly best while the 24run orthogonal array performs uniformly worst. Also, observe that the values of an and K2 no longer have an impact on relative design performance. The reason for this is that regardless of the value assumed by K2 = A12, the variance contribution to the MSE dominates the bias contribution. Thus, the actual value of JAI seems to have little practical influence on our recommendations. In some sense, these results are not very surprising. The superiority of the Res olution V design in terms of the mean squared error criterion is enhanced by the fact that it supports the estimation of the first order C x N interactions in the model. The latter three designs, however, are strictly main effect plans whose columns were chosen to yield a nonsingular P'P matrix. The choice of columns allowed the model to be estimated by the designs, but not as well as the Resolution V design. 3.6 Other Measures of Design Performance The integrated or average mean squared error criterion combines two aspects of design performance: average squared bias the average squared distance of f(x) from the true value of r(x), and average variance the average precision with which (x) estimates Tr(a), where average means averaged over the region of interest R. Rather than considering average design performance, the practitioner may be interested in the maximum mean squared error of (ax) within the region of interest for a given design. In this case we seek a design that minimizes the maximum mean squared error of T(x). From equation (3.10) the maximum mean squared error can easily be determined: maxMSE[f (x)] = max{ [bias((B(x))]2 + varAi['(x)]} XER 2ER = max{[tr(VE)]2 + 2tr(VE)2 XER 4K2 2r4 + tr[(I 0 x,,)VEV(Iq 0 X')] +  }. nA n m For the simulation study and four designs discussed above, the maximum MSE's were determined by randomly generating points over the cuboidal region defined by R = {x = (x1, X2, X3)' : 1 < x, < 1, i = 1,2,3}. The resulting maximums are displayed in Figures 3.8 and 3.9. In equation (3.10) recall that the variance portion varA[T(x)] was determined by averaging the variance of f(x) given by equation (3.8) over all possible As with length K. Hence the maximum mean squared error given above can be thought of as a maximumaverage, i.e., the maximum MSE over the region R after averaging over the parameter space of A. Instead of averaging var[af(x)] over all possible A's with length K, consider maximizing var[f(x)] over all possible A's with length K (Noble and Daniel, 1977): varM[T(x)] = max var[f(ax)] 2o4 = max {2tr( VE) + 4A'[Iq x]VEV[Iq0 A +  IAI=n n m 0,= o22= 0.25; K 2= 128 \  \ / 00 \ / L \ / S10 0.25; K / 128 0\ 5 .. .10 05 00 05 10 P12 al 2.250; o = 0.25; K= 128 0 05 / 10 .05 00 0.5 10 P,2 o, = 2.25; o0= 0.25; K2 128 \ \\ /. 10 05 00 05 10 S=0 22= 0.25; K2= 1152 1\ / \ / \ / \ / \ / ,.. \ / ..* .10 05 00 05 10 P 12 a, = 1.0;o2= 0.25; K2= 1152 \ 7/ N \ / ."." *1.O 05 00 0.5 10 P12 0O= 2.25;o02 = 0.25; Ic2= 1152 \ / ^\ / \ / / \ / ,, ...* '7 ..\ .10 0. O0 05 t0 SLegend: RESV ............ PB20 . PB24 _ OA24 Figure 3.8. The Effect of Correlation on the MaximumAverage MSE , ,=a22= 0.25; K 2= 1152 \/ \ / / ,/ 1 S = 1.0 a =2 025 .= 128 1.0 05 00 0.5 10 P 12 =2.250;o = 0.25; K= 128 \ / \ / \ \ // 1.0 0.5 00 0.5 10 P,2 0,,= 2.25;o22=0.25; 12= 128 \  S/ \ / \ / \ / \ . 1.0 0.5 0.0 0.5 1.0 P12 / \ /0 P /1 ai 1.0; = 0.25; 2= 1152 1.0 0.5 0.0 0,5 10 P 12 o,, = 1.0;o=220.25; K2=1152 \ F \ / \ / N / J J 1O 0.5 0.00 10 P 12 o,,=2.25;022 = 0.25; K 2= 1152 /  I  P,2 Legend: RESV ............ PB20 PB24 OA24 Figure 3.9. The Effect of Correlation on the Weighted MaximumAverage MSE o = 22= 0.25; K = 128 I where em denotes the largest eigenvalue of (Iq 0 )a)VEV(Iq x'). This yields a maximummaximum mean squared error criterion: maxMSE[4 (x)] = max{[bias( (x))]2 + varM[(xa)]} XER 2XR 24 = max{[tr(VE)]2 + 2tr(VE)2 + 42em + }. nrn For the simulation study and four designs discussed above, the maximum MSE's were determined by randomly generating points over the cuboidal region defined by R = {x = (Xl,X2,X3)' : 1 < xi < 1,i = 1,2,3}. The resulting maximums are displayed in Figures 3.10 and 3.11. Figures 3.8 3.11 illustrate that when weighted by the number of design runs, the Resolution V design provides the lowest values of the maximumaverage and maximummaximum mean squared errors. Interestingly, when unweighted by the number of design runs, the Resolution V design still provides the lowest values of the maximumaverage and maximummaximum mean squared errors except when P12 = cov(zl, z2) = 0. When the noise variables are independent the 24run Plackett Burman design provides the lowest values. Overall, the Resolution V design per forms best under these maximumaverage and maximummaximum criteria whether weighted or unweighted by the number of design runs. This differs from the conclu sions obtained with regard to the integrated mean squared error criterion. In that case, the Resolution V design was recommended under weighting while the 24run PlackettBurman design was the clear winner when unweighted by the number of de sign runs. Lastly, a comparison of the left and right hand columns of plots in Figures 3.8 3.11 indicates that once again the value of r2 (128 or 1152) has little impact on the relative performance of the four designs. 1.0 0.5 0.0 0.5 10 P12 al = 2.25; a22 = 0.25; 2= 128 N / / / S. *\. ..*" 10 05 0.0 0.5 10 11 = 22= 0.25; K'= 128 \ \ / \ / \ / \ / ..; 10 0.5 00 0.5 1.0 P 12 a0, = 1.0; a22 = 0.25; K2= 128 / ; // \ / "';.. \ ..' ^ ~~\ ': / **'/ > Legend: RESV .......... PB20  PB24 _. OA24 Figure 3.10. The Effect of Correlation on the MaximumMaximum MSE 0r11= 0 22= 0.25; K 2= 1152 S = 1.0;22 0.25; ic 2= 1152 \ \ / \ / \/\ / ., \ / .. " 10 05 00 05 10 S\/ / \ .' / 1.0 0.5 00 0,5 1.0 P.2 Oau = 2.25; o= 0.25; c2 = 1152 \ // \ / \. , 1,0 0,. 0.0 0.5 10 a, = o~= 0.25; = 128 \/ S12 11 1.0; 22 0.25; 2 128 10 05 00 05 10 P 12 O11 =1.O;o2z = 0.25; Kc 128 \/ \ / \ / \ / i i i 1 0 0. 00 05 10 P12 a11=2.25;o=o0.25; K= 128 \ / N / 0'05 00 1) 10 1 0 .0.5 0.0 0.5 1.0 S11= 22= 0.25; K2= 1152 10 _05 00 0.5 10 \ / S225 0.25 1152 00 0 / S. \ I ,/ 10 0,5 00 05 10 P,2 ao =21.0;o =0.25; K2 1152 X 7I \ / / \. /  10 .05 00 0.5 0 r 12 r 12 Legend: RESV ............ PB20 . PB24 . OA24 Figure 3.11. The Effect of Correlation on the Weighted MaximumMaximum MSE 3.7 Practical Application In this section we will illustrate our methodology for an example discussed in Box and Jones (1992). They consider a setting in the food manufacturing industry in which an experimenter is seeking the best recipe for a cake mix to be sold in a box at the supermarket. For this example, we suppose that there are four control factors, flour (F), shortening (S), sugar (G), and egg powder (E). Next suppose the experimenter anticipates that a consumer's oven temperature may be biased up or down and that a consumer may overcook or undercook the cake. Thus, there are two noise variables, time (t) and temperature (T), and it is of interest to determine a recipe, i.e., settings of F, S, G, and E that yield a goodtasting cake even if the consumer strays somewhat from the baking time and temperature recommendations on the box. For this combined array experiment having four control and two noise variables, assume we are interested in fitting the model y = 3o N+ O13iX + /2X2 + #3X3 + /4X4 + /l1Zl 7222 + Alxlzli + A21x2z1 + A31X3Z1 + A41x4z1 + A12XlZ2 + A22X2Z2 + A32X3Z2 + A42X4Z2 + 6. Let us consider six reasonable experimental designs that can be used to fit this model: Design 1: 32run Resolution VI Fractional Factorial Design (261), Design 2: 32run Taguchi Crossed Array Design (241 22), Design 3: 22run Doptimal design recommended by Shoemaker et al. (1991) (this design was intended to be a smaller competitor to Taguchi's 32run crossed array and was generated using RS/DISCOVER to fit a model that also included the interactions X1X2, x1x3, x1x4, and zlz2), Design 4: 22run Doptimal Design generated by SAS's OPTEX procedure, assuming Shoemaker's model, Design 5: 22run Doptimal design generated by SAS's OPTEX procedure, assuming our model, and Design 6: 16run Doptimal design generated by SAS's OPTEX procedure, assuming our model. In order to compare these designs utilizing the integrated mean square error cri terion developed in Section 3.3, the experimenter must specify V and K'2 in equation (3.10). For this example, let us assume that the recommended baking time and tem perature are 40 minutes and 325F, respectively. We anticipate that the consumer may bake the cake for 35 to 45 minutes at an oven temperature between 300F and 3500F. Hence we may estimate all = var(zi) using (range/4)2 = (2.5 min)2 and a22 = var(z2) using (range/4)2 = (12.5F)2. Assuming each variable is studied at two levels, denoted by 1 and +1 in coded design units, then this corresponds to taking ano = 22 = (.5)2 = .25. Also, it would seem reasonable to assume that in this example the noise variables are negatively correlated. This is because we be lieve that a consumer is likely to know for instance that his or her oven temperature is biased upward and will try to compensate by cooking the cake for less than the recommended amount of time. Thus we will consider P12 equal to 0.5, 0.25, and 0 which respectively indicate strong negative correlation, weak negative correlation, and independence between cooking time and temperature. To specify K2, we recall from the simulation study that the results were robust to the actual choice of JAI, so we will assume that the magnitude of each element of A is four times the magnitude of o,. If we take a, = 1, then 02 Al2 = nA(16ac) = (10)(16) = 160. The IMSEs for the competing designs were calculated using Monte Carlo integra tion with points generated over the cuboidal region R, defined in coded variables by R = {a = (xi,x2,X3,X4)' : 1 < xi < 1, i = 1,2,3,4}. The results are displayed in Figure 3.12. Since the designs have different numbers of runs, we display the graphs first unweighted and then multiplied by the number of runs to allow for fair com parisons. The Taguchi design is not present in the graphs since its P'P matrix is equivalent to the P'P matrix from the fractional factorial design. Thus, the designs have equivalent C matrices and so the MSEs calculated from equation (3.10) would be equivalent apart from simulation error. From Figure 3.12 we note that the 32run resolution VI fractional factorial design has the smallest bias and MSE in both the weighted and unweighted comparisons. Interestingly, Shoemaker's Doptimal design generated using RS/DISCOVER has sig nificantly larger bias and mean squared error than the 22run Doptimal design we generated using the OPTEX procedure in SAS. Hence, it seems that the algorithm for generating Doptimal designs can have some impact on our criterion. Another interesting feature is that the 22run Doptimal design that was generated to fit our smaller model performed worse than the two 22run designs used to fit Shoemaker's larger model. Lastly, note that the 16run Doptimal design, the smallest design that is allowed by our criterion since n m must be greater than zero, performed significantly worse than the other four designs, regardless of weighting. After remov ing the 16run Doptimal design from consideration, Figure 3.13 provides a better comparative display of the remaining four designs. 3.8 Conclusions The combined array is an experimental design used for finding conditions that achieve a target condition for a response of interest while also minimizing squared error loss. Under the assumption of an additive error structure, the appropriate squared error loss involves the estimated process variance. O = 22= 0.25 .   0, 04 03 402 01 00 5 04 ,3 02 01 010 P12 In = o22 = 0.25 S 04 3 02 01 00 P 12 a I= 22= 0.25 05 04 0 02 4O1 00 P 12 Legend: DOPT16 ............ SHOE22  DOPT(SHOE22) DOPT22 RESVI Figure 3.12. The Effect of Correlation on the Average Bias and Average MSE (Un weighted Weighted) 0 =11 C 22= 0.25   ....    0. 0 03 0.2 01 0,0 P12 011= 022= 0.25 . . 0.5 0.4 03 02 01 0.0 P12 0.=0 22= 0.25  .. ... .. . . 05 04 03 02 0.1 00 P 12 n= g 22= = 0.25 N.  N ". " 2... .. ^ '''' ^ 05 04 03 02 01 00 Legend: ........... SHOE22 . DOPT(SHOE22) DOPT22 RESVI Figure 3.13. The Effect of Correlation on the Average Bias and Average MSE (Un weighted Weighted) 5, ~ Our research establishes that the combined array yields a biased estimate of the loss. This bias is maximized at the perimeter of the region of interest and can be relatively large if the design is highly fractionated, the noise factors are highly variable in the process, or the inherent variability is large. In terms of the mean squared error of the estimated loss, our research indicates that the variance of the estimated loss tends to dominate the bias. We have shown that the nature of the correlation structure can have a significant impact on the mean squared error. However, the relative comparison of competitor designs through our mean squared error criterion was impacted very little by the choice of V, the variance covariance matrix of the noise factors, or n, the length of the vector of coefficients. On the whole, when the number of design runs is taken into consideration, Resolution V and higher fractions of 2k factorial designs appear to produce estimated losses with the smallest mean squared error. CHAPTER 4 ANALYSIS OF THE COMBINED ARRAY WHEN THE ASSUMED MODEL IS UNDERSPECIFIED 4.1 Introduction In Chapter 3, the properties of the estimated process variance and the resulting conclusions regarding various combined arrays were derived under the assumption that the presumed model for the response was correct. These derivations represented the first formal attempt at considering the properties of the combined array. Com bined arrays are currently being used in practice but little attention has been devoted to the properties of the estimator of the process variance. The importance of consid ering the properties of the process variance estimator lies in the fact that the analyst considers predicted values of the process variance to determine the robust control factor settings. We showed in Chapter 3 that the estimator of the process variance is biased even though the presumed model was correct. This led to the derivation of a mean squared error criterion which was used to determine which combined arrays among a given set produced process variance estimators having the smallest mean squared error for a particular experimental situation. It is widely recognized that the fitted model serves only as an approximation of the true model and the practitioner prefers the fitted model to be as simple as possible. Hence the fitted model is often times an underrepresentation of the true response surface. Thus, it is of interest and importance to determine the impact of underspecifying the model on the properties of the estimated process variance. Using the foundational notation and methodology that have been set up for the correct model case, in this chapter we will extend our results to the case where the model is presumed to be underspecified. Our goal once again is to derive the mean squared error of the estimated process variance in order that an integrated mean squared error criterion can be used to evaluate the performance of combined arrays. Ultimately, we want to specify which combined arrays yield estimates of the process variance that are least affected when the underlying model may have been underspecified. We will begin by formulating the problem of model misspecification in its most general setting, and in the next chapter discuss special cases of interest within the details of various examples. 4.2 Properties of the Estimated Process Variance Suppose an experimenter fits the following model to a quality characteristic of interest: y = 00 + x'i31 + z'y71 + x 'lAlli + e (4.1) where /#o is the intercept, x, is a pi x 1 vector representing the appropriate polynomial expansion of the control variables of order dc, for the assumed model, $1 is the pi x 1 vector of unknown parameters associated with the control factors, zl is a ql x 1 vector representing the appropriate polynomial expansion of the noise variables of order ds, for the assumed model, y is the q, x 1 vector of unknown parameters associated with the noise factors, Aln is the pi x ql matrix of unknown parameters associated with the control bynoise interactions (some elements of All may be zero), and e is the random error. However, let 7r be the true response represented by 77= 03o+x~ +x2I32+z~7l1+z2'yz2+'Anl+xlA12z2+a2A21z1+zx2A22z2+E (4.2) where /o, a1, i1, z3 ~y1, A11, and e are as defined above and X2 is a p2 x 1 vector representing the appropriate polynomial expansion of the control variables of order dC2, f32 is the associated P2 x 1 vector of unknown parameters, z2 is a q2 x 1 vector representing the appropriate polynomial expansion of the noise variables of order dn2, Y2 is the associated q2 x 1 vector of unknown parameters, A12 is the matrix of unknown coefficients of the controlbynoise interactions between the control variables in xa and noise variables in z2, A21 is the matrix of unknown coefficients of the controlbynoise interactions between the control variables in x2 and noise variables in z1, and A22 is the matrix of unknown coefficients of the controlbynoise interactions between the control variables in a2 and noise variables in z2. Note that the results that follow in this chapter do not depend on the orders of the polynomial expansions of the control and noise variables in the fitted and true models. However, the smallest model underlying a combined array experiment must contain first order terms in the control (C) and noise (N) variables as well as first order controlbynoise (CxN) interactions. Only by exploiting the CxN interactions can the experimenter determine the effect that changing the levels of the noise variables has on the control variables. Hence, some but not all of the elements of Aln may be assumed a priori to be zero in the fitted model, while some or all of the elements of A12, A21 and A22 may be taken to be zero in the true model. In Chapter 5 we will take d,1 = dc, = 1 and d, = d,, = 2 in order to investigate the impact of underspecification in the control variables alone, noise variables alone, controlby noise interactions alone, and some combinations of the three. For models (4.1) and (4.2) we shall assume that the es are independent and iden tically distributed with mean zero and variance or. For the purposes of a particular planned experiment, the noise factors shall be treated as having fixed effects. How ever, in the process we shall assume zl follows a multivariate normal distribution with E(zi) = p, and var(zi) = V1. Likewise we assume z2 follows a multivariate normal distribution with E(z2) = 12 and var(z2) = V2. We shall assume aui = var(zi) and oij = cov(zi, Zj) are known and depending on the order of the vectors zi and z2 additional elements of Vi and V2 may be derived through differentiation of the moment generating function of the multivariate normal distribution. These compu tations are shown for a special case in Appendix D. Furthermore, the error terms and noise variables shall be assumed to be independent. To develop the prediction properties of the estimated process variance, it will be helpful to establish a set of notation which follows: A'/Ia =[71,,A'1] S= (1,'1) A11 = vecAlna 1=1 (00,"1,Al) A12a = A12, A2] A12 = vecAi2a A21 = vecA21 A22 = vecA22 '2 = ( 2 12 ) 211 \22 02(Xi) = 2 + A'121 = A'2aXia = (Iq2 0 Xja)A12 A212 = ( 0 2)21 A2 = (q2 0 X)A22 n = total number of design runs mi = total number of unknown coeficients in model (4.1) = (1 + pi)(l + qi) M2 = total number of unknown coeficients in model (4.2) = (1 + p1 + P2)(1 + q, + q2) Using the above notation, models (4.1) and (4.2) can be expressed in the more general forms given by: y = P101 + e (4.3) and ri = P PI 2 + 2 2 + e (4.4) where P1 is the n x mi (mi = 1 + pi + ql(pi + 1)) model matrix whose elements are known functions of the settings of k input variables determined by an n x k design matrix, and P2 is the n x (m2 m1) (m2 ml = p2 + q2(P2 + 1)) matrix associated with terms not present in the fitted model but present in the true model. Assume the appropriate combined array experiment has been run to estimate the linear model given by (4.3). Using ordinary least squares to estimate the model coefficients, the resulting estimate of the coefficient vector is i1 = (PP1i)P1y Now if the true model is given by (4.4) and we assume the noise variables to be fixed in the experiment, then E[y] = P101 + P202. Thus the prediction properties of 4i are as follows: E[] = E[(P'P)1 P'y] = (P'Pi)'PIE[y] (PIP)lP(PI'1 + P212) = (PlPI)1PlPI,1 + (PIP1) PIP2b2 = 01 + A 2 (4.5) where A = (P'P1) P'P2 is called the alias matrix and var[ ] = P)1 Closer consideration of equation (4.5) will enable us to derive the prediction properties of A11: E[I]= E(qI) = (3 A A12 E(AI) A11 A 2 1 \ 22 / Thus, it follows that E[An1] = All + A*02 (4.6) where A* is formed by deleting the first p, + 1 rows of the alias matrix A. Now, since k1(Xi) = (Iqi 0 x31)A1, it follows from equation (4.6) that E[ 1(xi)] = (Iq 'ia)(A11 + A*02) = 41(z1) + (Iq 0 .ia)A*I2. (4.7) We are now in a position to derive the prediction properties of the estimated process variance which we denote by ,(xl). If the fitted model is given by (4.1) and the noise variables are truly random in the process, then analogous to equation (3.2) the estimated process variance may be expressed as (Mxi) = var[y] = 4~(xl)Vpk(xi) + &o. If the true model is given by (4.2), the derivation of the expected value of Ty(xa) follows from (4.7) E[,ry(ax)] = tr(ViSi) + E[ j(xl)]'ViE[1l(xl)] + o' = tr(vi i) + [i('i) + (Iq 0 ,o)A* '']V1[1(,1i) + (Ij xi)A*l2] + 02 = tr(Vix) + 0Y1(Xi)V1i1(ai) + 204((1)V1(Iq, 0 a'j)A*' + 2A*(IJ, 0 Xla)Vl(Iql, la)A*02 + 0 (4.8) where E, = var[ 1(xi)]. Also, from (4.2) we can derive the true process variance, which we denote by r,((x), x = (x'1,x')', as follows: ,(x) = var[y] = var[(i' + xA11 + 'A21)z1 + (72 + ,'A12 + x'A22)z2 + E] = (7' + xAnll + x'2A21)Vi(7 + BAllA + '2A21)' + (7 + 'A12 + x'A22)V2(Y72 + A12 + x'A22)' + 2(7'i + An11 + x2A21)V12(72 + 1xA12 + x'A22)' + o = 4'(xi)V11(xi) + 20'(x1)VIAa2 + X1A21VA'21 2 + '2(al)V202(al) + 24'(aX)V2A2x2 + 2xA22V2A'222 + 2(0((xi) + x' A21)V12( 2(X) + A'222) + U2 (4.9) where V12 = cov(z1,z2) and the terms 0'(xa) and '2(x1) are explained in the list of notation given on pages 62 and 63. 4.3 Expressions for the Bias and Variance of the Estimated Process Variance From (4.8) and (4.9) note that the estimated process variance is biased. In Ap pendix B it is shown that the bias can be conveniently expressed as bias[fy(axi)] = E[f,(xi)] r,(x) = tr(V1 E) + a'Ba (4.10) where a' = (A,, ip) and the matrix B is is given by equation (B.3). In order to develop a mean squared error criterion for estimating the loss incurred by using ,(axi) to estimate r,(ax), we must also derive the variance of fy(xi). If in model (4.4) we assume that e is distributed as multivariate normal with mean 0 and variance oIr, then in Appendix C it is shown that the variance of the estimated loss can be conveniently expressed as var[y,(xi)] = 2tr(ViEl)2 + 4a'Wa + 2a (4.11) n m1 where a'= (A', i') and the matrix W given by equation (C.2). Both the bias and the variance of Ty(x1) given by equations (4.10) and (4.11) respectively depend on the unknown parameter vector a. We proceed here as in Chapter 3. Let us assume that the length of the vector a is known where K = lal. Define U to be the surface of the unit sphere and ds, to be a differential on this surface. Then by averaging the bias and variance over all possible as with length K it is possible to eliminate the dependence on a. The averaged bias and variance are given by 2 KK2 biasA[y(xl)] = K2K bias[f,(xai)]ds, = tr(ViSi) + trB (4.12) JU n, and varA[y(Xi)] = 2K J var[fy(xi)]dsa JU 4K2 2U4 = 2tr(ViEi)2 + trW + (4.13) no n mi where K = [fu dsa]', and n, is the number of elements in a not assumed a priori to be zero. Lastly, from equations (4.12) and (4.13) we obtain the following equation for the average mean squared error of Ty(xi) which can be readily calculated at any point X1: K2 MSEA['y(xl)] = [tr(ViE) + trB]2 4K2 294 + 2tr(ViEi)2 + tW +  (4.14) n, n mi 4.4 The IMSE Criterion and Other Measures of Design Performance As discussed in Chapter 3, a reasonable method for comparing designs is to eval uate the integrated mean squared error (IMSE) of Ty(xl) for each design. The IMSE is given by IMSE[fy(a)] = K MSEA[?y(x1)]dx where R is the appropriate region of interest, K = [fRdx~]', and MSEA[y1(xa)] is given by equation (4.14). It may also be of interest to compare designs by evaluating the maximum mean squared error of N(xE) which is given by max MSEA[f'y(Xl)] = max{[biasA(,y(Xl))]2 + varA[fy(X1)]}. XiER X1ER Since MSEA[y1(xl)] was determined only after averaging bias[f',(xi)] and var[fy(xi)] over all possible as with length K, this maximum mean squared error can be thought of as a maximumaverage. An alternative maximum mean squared error criterion can be determined by first maximizing the bias and variance of fy(ax) over all possible as with length K as follows: biasM[y(xl)] = max[tr(ViEi) + a'Ba] = tr(ViEl) + K2eM(B) where eM(B) denotes the largest eigenvalue of B, and 24 varM[Ty(i)] = max[2tr(ViEi)2 + 4a'Wa +  laI=K n mi = 2tr(Vl~l)2 + 4K2eM(W) + 2 n m1 where eM(W) denotes the largest eigenvalue of W. This yields a maximummaximum mean squared error criterion given by max MSEM['y(xl)] = max{[biasM(jy(xl)]2 + varM[y (a )]}. X1ER XiER We are now ready to extend the case study and practical example that were discussed in Chapter 3 to the situation where the presumed model is underspecified. CHAPTER 5 SPECIAL CASES OF MODEL UNDERSPECIFICATION 5.1 Introduction In this chapter we consider special cases of model misspecification. For all cases we will begin by assuming that the fitted model is first order in the control (C) and noise (N) variables and contains first order controlbynoise (C x N) interactions. In other words, in model (4.1) we take x1 = (xi,..., ,Xp)' to be a vector of order d,, = 1 and z1 = (1, ,. .. zq)' to be a vector of order d, = 1. We will consider eight different "true" models containing the same terms in the fitted model plus additional higher order terms not present in the fitted model. Essentially, we assume the fitted model is underspecified. The additional higher order terms that will be incorporated into the true model are as follows: Case 1 (C2): The true model contains pure quadratic terms in the control factors and first order controlbycontrol interactions. Case 2 (N2): The true model contains pure quadratic terms in the noise factors and first order noisebynoise interactions. Case 3 (C2 x N): The true model contains interactions between noise factors and second order control effects. Case 4 (C x N2): The true model contains interactions between control factors and second order noise effects. Case 5 (C2, N2): The true model contains pure quadratic terms in the control factors, first order controlbycontrol interactions, pure quadratic terms in the noise factors, and first order noisebynoise interactions. Case 6 (C2, C2 x N): The true model contains pure quadratic terms in the con trol factors, first order controlbycontrol interactions, and interactions between noise factors and second order control effects. Case 7 (N2, C x N2): The true model contains pure quadratic terms in the noise factors, first order noisebynoise interactions, and interactions between control factors and second order noise effects. Case 8 (C2, N2, C2 x N, C x N2): The true model contains pure quadratic terms in the control factors, first order controlbycontrol interactions, pure quadratic terms in the noise factors, first order noisebynoise interactions, interactions between noise factors and second order control effects, and interactions between control factors and second order noise effects. These eight cases cover a variety of experimental situations that the practitioner may be concerned about. Under differing circumstances the experimenter may fear the presence of higher order terms in the control variables alone (Case 1), noise variables alone (Case 2), controlbynoise interactions alone (Case 3 and Case 4), or some combination of higher order terms in both the control and noise variables (Cases 58). 5.2 Reduced Expressions for the Bias of the Estimated Loss For each special case, the general expression for the bias of the estimated loss that was derived in Appendix B can be reduced into a simpler form. The simplification is a direct consequence of setting appropriate terms equal to zero in the bias expression. The elements of the B matrix in equation (B.3) depend on Q1, Q2,..., Q1 which are matrices of the quadratic and bilinear forms described in Appendix B. For the fitted and true models that we are considering, note that Q8, Q9, Q10, and Q11 are null matrices. This first simplification is a consequence of the fact that the fitted model contains only first order noise variable terms and the true model contains at most second order noise variable terms. In this case V12 = cov(zi, z2) = 0 under the assumption that z = (z1, z2)' follows a multivariate normal distribution with mean JA and variancecovariance matrix v=(VIl V12) V'12 V2 " Since V12 = 0, clearly Q8, Q9, Q10, and Q11 are null matrices. 5.2.1 Case 1 Reductions In case 1, the true model contains C, N, C x N, and C2 terms. Thus, in equation (4.2) we set 'Y2, A12, A21, and A22 equal to zero. This model is given by = /#0+ X'131 + x'232 + z1Y1 + x'1Alz + (5.1) Under these conditions, any terms involving M1, M2, and M3 would be eliminated, and so the matrices B1 and B2 in Appendix B reduce to B1 = Q1 and B2 = Q2. 5.2.2 Case 2 Reductions In case 2, the true model contains C, N, C x N, and N2 terms. Thus, in equation (4.2) we set 132, A12, A21, and A22 equal to zero. This model is given by ] = 00 + X111 + z'1iY + z1y2 + 'AnAzi + e. (5.2) Under these conditions, any terms involving M2 and M3 would be eliminated, and so the matrices B1 and B2 in Appendix B reduce to B1 = Q1 and B2 = Q2 Q5. 5.2.3 Case 3 Reductions In case 3, the true model contains C, N, C x N, and C2 x N terms. Thus, in equation (4.2) we set /2, Y2, A12, and A22 equal to zero. This model is given by S= 0o+ x31 + z'iy1 + xAlzzi + x'2A21z + e. (5.3) Under these conditions, any terms involving M1 and M3 would be eliminated, and so the matrices B1 and B2 in Appendix B reduce to B1 = Q1 Q3 and B2 = Q2 Q4. 5.2.4 Case 4 Reductions In case 4, the true model contains C, N, C x N, and C x N2 terms. Thus, in equation (4.2) we set /2, 'Y2, A21, and A22 equal to zero. This model is given by 7 = o + a'/31 + z'y + a'Anzi + x'A12z2 + e. (5.4) Under these conditions, any terms involving M2 and M3 would be eliminated, and so the matrices B1 and B2 in Appendix B reduce to B1 = Q1 and B2 = Q2 Q5. 5.2.5 Case 5 Reductions In case 5, the true model contains C, N, C x N, C2, and N2 terms. Thus, in equation (4.2) we set A12, A21, and A22 equal to zero. This model is given by q1 = /o + x/31 + X'202 + z'721 + z2'72 + x'lAliz + e. (5.5) Under these conditions, any terms involving M2 and M3 would be eliminated, and so the matrices B1 and B2 in Appendix B reduce to B1 = Q1 and B2 = Q2 Q5. 5.2.6 Case 6 Reductions In case 6, the true model contains C, N, C x N, C2, and C2 x N terms. Thus, in equation (4.2) we set 72, A12, and A22 equal to zero. This model is given by =0 + X'131 + x2'32 + z[1 + xaAniz + x2A21 z + (5.6) Under these conditions, any terms involving M1 and M3 would be eliminated, and so the matrices BR and B2 in Appendix B reduce to B1 = Q1 Q3 and B2 = Q2 Q4. 5.2.7 Case 7 Reductions In case 7, the true model contains C, N, C x N, N2, and C x N2 terms. Thus, in equation (4.2) we set /32, A21, and A22 equal to zero. This model is given by 0= /0 + X'11 + z' Y1 + z2'2 + zAnzllz + zxA12 z+ (5.7) Under these conditions, any terms involving M2 and M3 would be eliminated, and so the matrices B1 and B2 in Appendix B reduce to B1 = Q1 and B2 = Q2 Q5. 5.2.8 Case 8 Reductions In case 8, the true model contains C, N, C x N, C2,N2, C2 x N, and C x N2 terms. Thus, in equation (4.2) we set only A22 equal to zero. This model is given by i = o/ + 1301 + X'202 + Z7 11 + z"'y2xAzi + x'2A21Z1 + a 'A12z2 + e. (5.8) Under these conditions, any terms involving M3 would be eliminated, and so the matrices B1 and B2 in Appendix B reduce to BI = Q1 Q3 and B2 = Q2 Q4 Q5. 5.3 Example In Chapter 3, we noted that even though the fitted model was presumed to be adequate, the process variance estimator was biased. The development of a mean squared error criterion for the estimated process variance and subsequent simulation study enabled us to compare a set of designs. In this section, where the fitted model is presumed to be underspecified, we will conduct another series of simulations to investigate the additional impact of model misspecification on the bias and variance of (xzi). For each of the eight special cases we will again consider a robust parameter design experiment consisting of three control variables and two noise variables, each being studied at two levels. For this example, the fitted model given by (4.1), in expanded form is: y = 3o + /1x1 + #2x2 + 03X3 + 71zi + 7y22 + A11xzi + A21x2zI + A31 x3z1 + A12X1Z2 + A22X2Z2 + A32X3Z2 + E. We will compare the same four designs utilized in Chapter 3 to fit this model: a 16run Resolution V fractional factorial design, 20 and 24run PlackettBurman designs, and a 24run orthogonal main effect plan. We will compare these designs by evaluating the integrated mean squared error of B(xi), and the maximumaverage, and maximum maximum criteria that were derived in Section 4.4. These criteria will be calculated using Monte Carlo techniques with points generated over the cuboidal region defined by R = {x= (x,x2,x3)' :1 < xi < 1,i= 1,2,3}. Recall from equation (4.14) that computation of the mean squared error of 'y(x() at any point xla depends on V1, the assumed known variancecovariance matrix for the noise variables in the fitted model, and the value of K2/n,. For this example with two noise variables in the fitted model, let j1i, a12, and a22 denote the elements of the V1. Using the same argument from Chapter 3 based on specifying the range of the noise variables, we will consider three cases: 1. anl = a22 = .25 (R1 = Ri, R2 = ); 2. all = 1; a22 = .25 A (t1 = Ri/2, R2 = R2); 3. all = 2.25; 022 = .25 # (iR = Ri/3, R2 = R2); for which the corresponding VI matrices are given by 1P 1 1 P12 I 1 2 P12 and V, = 3 [ 3 P12 V=4 P12 1 j' 12 P12 1/2 and 4 P12 1/3 To investigate the effect of the correlation structure of the noise factors on the bias and mean squared error, we will consider a sequence of values of P12 between 1 and 1. Lastly, note that in special cases 2, 4, 5, 7, and 8, where the vector z2 of second order noise effects is present in the true model, we must determine V2, the variance covariance matrix of these second order effects. The general form of V2 is derived in Appendix D. 5.4 Results and Discussion For each special case of the underspecified model we made use of the simplifications discussed above and used SAS/IML to compute the average bias and average mean squared error of T,(xai) for the four designs under consideration. Eight special cases, two values of K2, and three Vi matrices produced a total of forty eight sets of values. As in Chapter 3, we will display the biases and mean squared errors graphically, both unweighted and weighted (multiplied) by the respective numbers of design runs, in order to make conclusions and recommendations regarding the four designs. The graphical displays will illustrate the impact of: (1) the variances, all and U22, of the noise factors; (2) the correlation, P12, between the noise factors; and (3) the length of the vector a (2 = a 2). These plots are given by Figures 5.1 through 5.32. When observing the graphical displays, bear in mind the threefold intent of this simulation study: 1. To investigate whether the average bias and average mean squared error are "robust" to the values of '11, 0"22, P12, and K2. 2. To investigate the impact of model underspecification on the bias and mean squared error of (axi). 3. To determine which design performs best within each case and across cases of the underspecified model. By "performs best", we mean that the design has the smallest bias and/or mean squared error of fy(xi). Figures 5.1 and 5.2 display the average bias for case 1 in which the true model contains additional pure quadratic and interaction terms in the control variables only. In Figure 5.1, where the average bias is unweighted by the number of design runs, observe that the 16run Resolution V fractional factorial design (RESV16) performs best by a small margin over the 24run orthogonal array (OA24). It is important to note that the difference in the average bias associated with the RESV16 and OA24 designs is significant based on a 95% confidence interval computed using simulated standard errors. The scale on the vertical axis causes the difference to appear non significant. The 24run PlackettBurman design (PB24) is third best followed by the 20run PlackettBurman design (PB20) in last place. One interesting feature to notice is that the average bias for the PB20 decreases as p12 ranges from 1 to 1, while the average bias is nearly constant for the other three designs over the values of P12. Comparing the left and right columns of plots we note that apart from a change in scale the relative performance of the four designs does not depend on the value of K2. In Figure 5.2 where the average bias is multiplied by the number of design runs, we see that necessarily the RESV16 still performs best with the OA24 coming in second. Interestingly, the PB20 design performs worst when the noise variables are negatively correlated while the PB24 performs worst when the noise variables are positively correlated. Figures 5.3 and 5.4 display the average mean squared error for case 1. Interestingly, in Figure 5.3, where the average MSE is unweighted by the number of design runs, we see that the RESV16 performs best with one exception. When K2 = 128, acr = o 1 = a = 0.25; K = 1152 10 0.5 0.0 0.5 10 S12 a = 1.0; a22 = 0.25; K1= 128    1.0 05 00 05 10 P12 a0I = 2.25; a22= 0.25; i =_ 128 ............*.... 1.0 0.5 0.0 0.5 10 P 12 0 = 1.0; 22= 0.25; Ki 1152 ............ ....... Ii 1o0 405 0.0 05 10 P 12 o ,= 2.25;Y22 =0.25; K = 1152 ... .  10 0.5 00 P 12 05 1.0 .10 0.5 00 05 10 Legend: RESV ............ PB20  PB24 OA24 Figure 5.1. Case 1 (C2) The Effect of Correlation on the Average Bias o, = y,2 = 0.25; K =_ 128 0 = 0=2 0.25; = 128 . . : _  10 .0 00 0.5 10 S12 0,, = 1.0; o22 = 0.25; )= 128   10 05 00 05 10 P 12 0 I= 2.25;022= 0.25; x = 128 .. . . . S= ao22= 0.25; ic = 1152 .1.0 05 00 0.5 10 P 12 o1, = 1.0; 22= 0.25; i= 1152  T .  10 05 00 05 10 P 2 0a1,=2.25;022 = 0.25; ie= 1152 10 0.5 0,0 P12 0,5 10 1,0 05 00 05 10 OA24 Figure 5.2. Case 1 (C2) The Effect of Correlation on the Weighted Average Bias Legend: RES ............ PB20 . . PB24 o = 22= 0.25; = 128 " .......... . LO 0,5 0.0 0.5 10 P 12 a = 1.0;22 = 0.25; i = 128 .10 .05 00 05 10 P 12 1 = 2.25; a = 0.25; ? = 128 ........ ......... I . o0. 0.5 1L 0, = 022= 0.25; K 2 1152 10 0.5 00 05 10 P 12 ao1 = 1.0; 2= 0.25; i = 1152 O 05 0.0 05 10 P 12 a, = 2.25;oa = 0.25; e = 1152 .i n5 on 05s O Legend: RESV ............ PB20  PB24 OA24 Figure 5.3. Case 1 (C2) The Effect of Correlation on the Average MSE i0      0 = o,, = 0.25; c 1152 10 0.5 00 05 10 P 2 0,, =1.0; 0 =0.25; = 128 \ .10 05 00 05 10 P 2 ao, 2.25; 22= 0.25; i= 128 .. ... . 10 05 00 05 10 10 05 00 05 LO P12 ,, =1.0;o22= 0.25; Kd= 1152 10 I05 00 05 L Pl2 o,,= 2.25;o22 =0.25; = 1152  10 05 00 05 10 P,2 Legend: RESV ............ PB20 PB24 OA24 Figure 5.4. Case 1 (C2) The Effect of Correlation on the Weighted Average MSE S,,= 0,= 0.25; i = 128 " *  . a22 = .25, and .1 < P12 < .7 (top plot in left column), note that the OA24 performs best. Also, the PB24 is in third place for all situations except when r,2 = 1152 and all = U22 = .25 and all = 1.25, a22 = .25 (top and middle plots in right column) in which case the PB20 outperforms it for large positive correlation. In Figure 5.4 where the average MSE is multiplied by the number of design runs, the RESV16 performs best in all situations with the OA24 performing second best. The PB20 performs worst except for large positive correlation in which case the PB24 performs worst. The superiority of the performance of the RESV16 in case 1 can be explained by considering the alias matrix A*. For the RESV16, A* = 0. Consequently, it follows from Section 5.1.1 that B1 = Q1 = 0, B2 = Q2 = 0, W2 = 0, and W3 = 0, where the terms B1 and B2 are explained in Appendix B and W2 and W3 are explained in Appendix C. As a result, the bias and mean squared error of Ty(xl) are equivalent to the bias and mean squared error of (xa) obtained for the case in which the fitted model was presumed to be correct. In other words, the second order control effects present in the true model have no additional impact on the bias and mean squared error of the estimated process variance when the RESV16 is utilized. For the OA24, we also have A* = 0, and so the average bias and average MSE are equivalent to that found in Chapter 3 for the correct model. However, the A* matrices for the PB20 and PB24 designs are nonzero. This serves to explain the increase in bias and MSE observed for these designs as compared to the case of the fitted model being correct. It also explains why the OA24 performs better than the PB20 and PB24. This differs from the results obtained in Chapter 3 where the OA24 was observed to perform worst. Figures 5.5 and 5.6 display the average bias for case 2 in which the true model contains additional pure quadratic and interaction terms in the noise variables only. In Figure 5.5, where the average bias is unweighted by the number of design runs, note that the average biases are negative and also that a very small difference separates S,,= 22= 0.25; 2= 1152 10 0 5 00 0,5 10 S12 oi =1.0;22 = 0.25; K= 128 1.0 05 00 05 10 P 1.2 oa = 2.25; 0a22= 0.25; K 2= 128 05 00 P 12 05 10 L0 .0,5 00 05 1.0 P 12 oa = 1.0; 22= 0.25; K = 1152 10 05O 0,0 P 12 01= 2.25;22 = 0.25; K = 1152 10 05 00 P 12 05 10 0.5 L Legend: RESV ............ PB20  PB24 OA24 Figure 5.5. Case 2 (N2) The Effect of Correlation on the Average Bias o = 022= 0.25; 2= 128 oa= 022=0.25; Ki= 128 .     / \ 10 0.5 0.0 05 1.0 P 12 11 =1.0;22 = 0.25; C2= 128 10 0,5 00 05 1.0 S12 1 = 2.25; 22= 0.25; K= 128 0 05 00 0.5 10 Legend: RESV ............ cH = 02= 0.25; K = 1152 ..       10 05 00 0.5 1.0 P12 1 = 1.0; o = 0.25; K= 1152 / ". 10 05 0.0 0.5 10 Pl2 S=2.25; a 22= 0.25; K2= 1152 a, = 2.25; a22 = 0.25; K< 2 = 1152 10 .05 0.0 05 10 PB20 .. PB24 OA24 Figure 5.6. Case 2 (N2) The Effect of Correlation on the Weighted Average Bias ... ... ..... ... ... ..... .. .. . the best and worst designs. However, after multiplying by the number of design runs, one can observe in Figure 5.6 that the RESV16 performs best followed by the PB20. The OA24 and PB24 have nearly equal weighted average biases. Figures 5.7 and 5.8 display the average MSE for case 2. In Figure 5.7 where the average MSE is unweighted by the number of design runs, note that a very small difference separates the best and worst designs although it appears that the Plackett Burman designs are performing better then the RESV16 and OA24. However, after multiplying by the number of design runs, Figure 5.8 shows that the RESV16 has the smallest average MSEs followed by the PB20 and PB24 and the OA24 comes in last. In case 2 we observed again that A* = 0 for the RESV16 and OA24. Consequently, it follows from Section 5.2.2 that Bi = Q, = 0 and B2 = Q2 Q5 = Q5, where Q5 is positive definite and not dependent upon the design. Hence the additional impact of the presence of second order noise effects in the true model is influenced by Q5. The A* matrices for the PB20 and PB24 are nonzero. However, A* contains only one nonzero column and so most of the elements of the matrices Q1 and Q2 are zero, which serves to explain the negative biases that were observed for all the designs considered. Figures 5.9 and 5.10 display the average bias for case 3 in which the true model contains additional interactions between noise factors and second order control effects. In Figure 5.9, where the average bias is unweighted by the number of design runs, one can observe that the PB24 performs best except for large negative correlation in which case the PB20 performs best. The OA24 performs uniformly worst. However, after multiplying by the number of design runs, Figure 5.10 shows that the RESV16 performs uniformly best and the OA24 performs uniformly worst. Figures 5.11 and 5.12 display the average MSE for case 3. In Figure 5.11, where the average MSE is unweighted by the number of design runs, we see that again the PB24 performs best except for large negative correlation in which case the PB20 o =o 22,= 0.25; 2= 1152 10 05 0,0 0.5 10 P2 12 11 =1.0; o22 = 0.25; K 2= 128 12 12 ~ 10 0,5 00 0.5 10 P t2 a = 2.25;a022= 0.25; K 2= 128 1.0 05 00 P 2 05 1.0 Legend: RESV 1.0 0.5 00 05 10 P 12 o = 1.0;a22= 0.25; K 2= 1152 1.0 0.5 00 05 1 P 12 S= 2.25; 22 = 0.25; 2= 1152 LO 05 0,0 05 10 ............ PB20  PB24 OA24 Figure 5.7. Case 2 (N2) The Effect of Correlation on the Average MSE (o = 2,= 0.25; K 2= 128 a = o22= 0.25; K 128 P 12 0 = 1.0; a = 0.25; K 2= 128 \ /' '\ 5 P t2 a = 2.25; 22= 0.25; K= 128 \ 0 10 .05 00 05 10 0 =022= 0.25; K = 1152 Lo 45 00 05 10 P12 1.0 = 0.25; K 1152 '/ .10 05 00 05 10 P12 oG = 12.25; 22= 0.25; K= 1152 10 405 00 05 10 Legend: RESV ............ PB20  PB2 OA24 Figure 5.8. Case 2 (N2) The Effect of Correlation on the Weighted Average MSE 1,,= Oa= 0.25; i2= 128 10 05 00 05 10 p12 a = 1.0;a =0.25;K2= 128 1.0 0,5 0.0 0.5 1,0 P 12 ao, =2.250;o= 0.25; K 2= 128 LO 05 00 0,5 1.0 10 .0.5 00 05 1) o =o=22 0.25; K2= 1152 . .  10 05 0.0 0.5 1.0 P 12 a = 1.0; a 22= 0.25; = 1152 10 05 0.0 0,5 1.0 P 12 a =2.25; 22 = 0.25; i = 1152 10 05 0.0 05 10 .10 .1)5 0.0 0.5 10 Legend: RESV ............ PB20  PB24 OA24 Figure 5.9. Case 3 (C2 x N) The Effect of Correlation on the Average Bias 1 =( =22 0.25; C 2= 128 10 05 00 05 10O P 12 0, =1.0;22 = 0.25; K2= 128 10 05 00 0,5 10 P 2 a, = 2.25; a22= 0.25; K 2= 128 a1 =a22= 0.25; K = 1152 10 405 00 05 10 S12 C,, = 1.0; 022= 0.25; K = 1152 10 4)5 00 05 10 p,2 0, = 2.25; a22 = 0.25; K 2= 1152 00 05 10 1.0 0.5 0.0 05 10 r 12 r 12 Legend: RESV ............ PB20  PB24 _ OA24 Figure 5.10. Case 3 (C2 x N) The Effect of Correlation on the Weighted Average Bias 1.0 05  ................. ......... ....... ...... 11, =a o= 0.25; K 2= 128 \ \ \ 1,0 05 0.0 0,5 1,0 P 12 011 =1.0; a = 0.25; K2= 128 \ 10 05 00 0.5 1.0 a I= 2.25; a22= 0.25; I 2= 128 \ 10 05 00 05 10 N\ N~ 10 0.5 00 0.5 10 a = 22= 0.25; K2= 1152 10 0.5 00 05 1,0 P 12 Ca1 = 1.0;a = 0.25; K 2= 1152 \ \ 1,0 0.5 00 0.5 1.0 (F= 2.25; 022 = 0.25; K 2= 1152 1'0 0.5 0.0 0.5 10 Legend: RESV ............ PB20 PB24 OA24 Figure 5.11. Case 3 (C2 x N) The Effect of Correlation on the Average MSE N N N o, = oa= 0.25; K2= 128 10 05 00 015 10 P12 0,, = 1.0; =0.25; K 2= 128 10 05 00 05 10 P 12 11 = 2.25; a22= 0.25; K2= 128 1.0 405 0.0 05 10 Legend: RESV ............ ,=1 22= 0.25; K= 1152 \ \ \ 10 05 00 05 1.0 S12 0,, =1.0; a = 0.25; K2= 1152 \ \. .10 05 00 0(5 1 I P12 01 = 2.25; 022 = 0.25; K2= 1152 10 05 00 05 10 PB20  PB24 OA24 Figure 5.12. Case 3 (C2 x N) The Effect of Correlation on the Weighted Average MSE N 'N. 'N. performs best. The OA24 performs uniformly worst. However, after multiplying by the number of design runs, Figure 5.12 shows that the RESV16 performs best except for large negative correlation in which case the PB20 performs best. The OA24 performs uniformly worst. In case 3, all of the A* matrices were nonzero. The overall conclusions regarding design performance are close to the recommendations obtained for the correct model case with the exception of observing that the PB20 performs best for large negative correlation. This was the first time we observed the PB20 to be the best design. Figures 5.13 and 5.14 display the average bias for case 4 in which the true model contains additional interactions between control factors and second order noise ef fects. In Figure 5.13, where the average bias is unweighted by the number of design runs, one can observe that the average biases are negative and the PB20 performs best except for large negative correlation in which case the PB24 performs best. The RESV16 and OA24 average biases are nearly equal and these designs perform uni formly worst. However, after multiplying by the number of design runs, Figure 5.14 shows that the RESV16 performs best with two exceptions. The top plots in the left and right columns show that the PB20 performs better than the RESV16 when the noise variables are positively correlated. The OA24 performs uniformly worst. Figures 5.15 and 5.16 display the average MSE for case 4. In Figure 5.15, where the average MSE is unweighted by the number of design runs, observe again the PB20 performs best except for large negative correlation in which case the PB24 performs best. The RESV16 and OA24 average MSEs are nearly equal and these designs perform uniformly worst. However, after multiplying by the number of design runs, Figure 5.16 shows that the RESV16 performs best with two exceptions. The top plots in the left and right columns where oll = .25 show that the PB20 performs better than the RESV16 except for large negative correlation. The OA24 performs uniformly worst. 