Teaching Optimal Design of Experiments Using a Spreadsheet.

Peter Goos
Universiteit Antwerpen

Herlinde Leemans
Katholieke Universiteit Leuven

Journal of Statistics Education Volume 12, Number 3 (2004), jse.amstat.org/v12n3/goos.html

Copyright © 2004 by Peter Goos and Herlinde Leemans, all rights reserved. This text may be freely shared among individuals, but it may not be republished in any medium without express written consent from the authors and advance notification of the editor.

Key Words:Linear Models; Microsoft Excel; Nonlinear models; Response surface experiment; Solver

Abstract

1. Introduction

Despite these attractive features, the optimal design theory has seldom received a considerable amount of attention in courses on experimental design, not even at the master's level. Often, the courses focus on the analysis of nicely balanced designs such as balanced complete and incomplete block designs and Latin square designs that can be analyzed using analysis of variance (ANOVA) methods. Typically, the experimental factors considered are qualitative. As such, a large spectrum of design problems and statistical models, having quantitative experimental factors and quadratic or higher order terms, is ignored. In addition, the design of experiments for estimating nonlinear statistical models is neglected.

Not only courses, but also textbooks on experimental design, (for example, Kuehl 2000; Montgomery 2000; Neter, Kutner, Nachtsheim, and Wasserman 1996; Oehlert 2000; Weber and Skillings 1999) pay little attention to the design of experiments involving quantitative variables. Typically, at most one chapter or section is spent on this kind of experiment, which is often referred to as a response surface experiment. The optimal design of experiments receives even less attention.

The purpose of this paper is to provide experimental design teachers with a simple way to introduce their students to the optimal design of experiments. It will be demonstrated how well-known spreadsheets and optimization software can be used to explain and solve the problem of designing efficient experiments for estimating both linear and nonlinear models. An important advantage of the approach is that it can be used in class interactively, avoids the "black box" approaches offered by statistical software programs and helps students to understand the problem and the solution. For these purposes, the students build the design matrix, the extended design matrix and the variance-covariance matrix of the estimators corresponding to simple concrete design problems using a spreadsheet. They can try to find a good experimental design by trial and error and for harder problems, by using the optimization software available in the spreadsheet. The examples employed in this paper are taken from different fields of study to demonstrate the usefulness of the optimal design approach in a wide range of research areas.

We have used the 2000 and 2002 versions of the Microsoft Excel spreadsheet and its standard add-in Solver for solving optimization problems. Spreadsheets containing all the examples described in this paper are available from Teaching Spreadsheets. We have included a short overview on the use of Solver for optimal experimental design in an Appendix A. The reader is referred to www.solver.com for more details.

2. An optimal design for a simple linear regression model

(1)

where x_i denotes the water temperature for the i^th child, is the design matrix of the experiment. This matrix contains the settings of the experimental variable, i.e. the water temperature, at each observation. Each of these settings can be chosen by the researcher. When this has been done, the reduction in pulse can be measured for each child (e.g., in beats/minute) and a regression line can be fitted to the data.

If the researcher fits a simple linear regression line, the model can be written as

(2)

where Y_i denotes the reduction in pulse for the i^th child. In matrix notation, this model can be written as

(3)

or

.

(4)

The matrix F is referred to as the extended design matrix of the experiment. Note that the (i,j)^th element of F can be computed as

.

(5)

The extended design matrix is needed to calculate the ordinary least squares estimator

,

(6)

the variance-covariance matrix of which is given by

,

(7)

where the constant represents the variance of the errors (i = 1, 2, ... ,6). The diagonal elements of the variance-covariance matrix are the variances of the ordinary least squares estimators for the intercept and the slope . In order to draw powerful conclusions about these parameters, it is desirable that the variances are small. The experimenter would therefore like to minimize the trace of (F^TF)^-1. A drawback of this approach is that it ignores the covariance between the estimators for the intercept and the slope . Ideally, the covariance between the estimators is close to zero, since in this case, the t-tests on the individual coefficient estimates can be made independent of each other. A design criterion which takes into account both the variances of the estimators and and their covariance is the determinant of the variance-covariance matrix (F^TF)^-1. For the six selected temperatures, the determinant of (F^TF)^-1 amounts to 0.000381. This can be verified by using a spreadsheet, constructed following the steps below:

1. Type in the design matrix X of the experiment. The dimension of X is n x m, where n is the number of observations and m is the number of experimental variables. As temperature is the only variable in this experiment which involves six observations, X is a 6 x 1 matrix containing the set of temperatures.

2. Construct the n x p extended design matrix F, where p is the number of unknown model parameters, using (5). In this example, F is a 6 x 2 matrix because there are two unknown parameters: and . The first column of F, which corresponds to , contains ones and the second column, which corresponds to , contains the set of temperatures.

3. Transpose the extended design matrix to obtain F^T.

4. Calculate (F^TF).

5. Calculate (F^TF)^-1.

6. Calculate the design criterion value, e.g. the D-criterion value, i.e. the determinant of (F^TF)^-1.

For those readers, who are not familiar with matrix computations in Microsoft Excel, we provide a succinct tutorial in Appendix B. Now, it is easy to investigate the impact of changing the six temperatures on the design criterion value. Manually finding the best design for the experiment, that is the set of temperatures that produces the smallest determinant, is in general not easy. For the present example, large reductions in the determinant can be realized by changing a few temperatures. In order to find the best possible design, Solver can be used (see Appendix A). If no temperatures below 45 or above 70 degrees are allowed, e.g. because the water would then be too cold or too warm, the smallest possible determinant is 0.000178 and is obtained by using the temperatures 45 and 70 three times each. This design is called D-optimal. The relative efficiency of the the original design compared to the D-optimal one equals (0.000178/0.000381)^1/2=68.35%. This means that the original design has to be replicated 1/0.6835=1.4630 times in order to achieve a confidence ellipsoid about as small as in the D-optimal design. The D-optimal design is therefore 46.30% more D-efficient than the original design. In general, the relative D-efficiency of two designs is defined as the ratio of the two determinants raised to the power 1/p, where p is the number of unknown model parameters.

3. Optimal designs for multiple linear regression models

3.1 Quadratic regression on one variable

(8)

where Y_i denotes the threshold current observed at the i^th trial and x_i is the corresponding Al mole fraction. In matrix notation, this model can be written as

.

(9)

As in the medical experiment in Section 2, a design consisting of equally spaced points can be chosen. For instance, if the Al mole fraction is allowed to vary between 0.15 and 0.60, the values 0.15, 0.24, 0.33, 0.42, 0.51 and 0.60 can be used in the experiment. It can, however, be verified applying Solver that this design is not optimal. It turns out that the best design possible for estimating the quadratic model (8) has two observations at 0.15, two at 0.375 and two at 0.60.

Depending on the starting design chosen, the optimal design might not always be found by Solver. Instead, Solver might converge to a locally optimal design. If, for example, the points 0.15, 0.15, 0.15, 0.40, 0.60 and 0.60 are used as a starting design, Solver will produce a design with 0.15, 0.15, 0.15, 0.375, 0.60 and 0.60. If six equidistant values between 0.15 and 0.60 are used, the global optimum is obtained. The probability of ending up with locally optimal designs increases with the number of observations and with the number of experimental variables investigated. It is therefore recommended that more than one starting solution, e.g. one for each student in the class, is tried and that the resulting designs produced by Solver are compared. How starting designs and (locally) optimal designs can be stored is described in Appendix A.

It may happen that Solver is not able to improve the initial starting design. This is often due to the fact that a singular design is used as a starting design. If, for example, only two different Al mole fractions were used in the above example, then the matrix (F^TF) would be singular and not invertible. This causes numerical problems and causes Solver to remain stuck in this starting design. In order to avoid such singular starting designs, at least as many distinct Al mole fractions are needed as there are unknown model parameters in the model. For the present example, three parameters need to be estimated, so that at least three Al mole fractions should be used in the starting design. Another reason why Solver may not improve upon the starting design is that the default precision and convergence values used are too large. Therefore, it may be necessary to adjust these values (see Appendix A).

3.2 A quadratic model in two variables

Consequently, the 9 x 2 design matrix and the 9 x 6 extended design matrix for this problem are given by

and

respectively.

The optimal design for this experiment is a 3² factorial design, which is a design with observations at the factor level combinations (10%,10^oF), (10%,20^oF), (10%,30^oF), (20%,10^oF), (20%,20^oF), (20%,30^oF), (30%,10^oF), (30%,20^oF) and (30%,30^oF). Again, it is possible that a suboptimal solution is obtained. It may therefore be necessary to try several starting designs. Our experience suggests that in this example only a couple of tries are needed to find the global optimum. As illustrated in Figure 1, scatter plots can be used to display the nine design points and to compare the starting designs and the optima found. The first starting design in Figure 1a leads to the locally optimal design in Figure 1c, whereas the second starting design in Figure 1b leads to the globally optimal design in d. In this way it can be verified that a starting design that does not fill the design space well usually yields a suboptimal result.

(a) Starting design 1	(b) Starting design 2
(c) Locally optimal design	(d) Globally optimal design

Figure 1. Two starting designs and the resulting optimal designs produced by Solver.
(Note that two points are duplicated in panels (a) and (c).)

4. Constrained mixture experiment

The electric resistivity of the powder did not depend on the total amount of the mixture but only on the relative proportions of the three components. Each component is therefore restricted to lie between 0 and 100%, i.e. 0 x_i 1. In addition, the proportions in a mixture experiment have to add up to 100%, so that

It was also required that

0.2	x1	0.8,
0.2	x2	0.8,
0.0	x3	0.6.

Now, assume that the model is given by the first order Scheffè polynomial

(10)

and that nine observations are available for estimating this model. For the first order Scheffè model, the design matrix and the extended design matrix are identical:

Again, a spreadsheet can be used to evaluate different design choices. For example, the design displayed in the left panel of Table 1 produces a determinant of 4.2409. Although this design covers the design region well, it is far from being D-optimal. The D-optimal design for estimating the model under investigation is displayed in the right panel of Table 1. It has three observations at each of the corner points of the design region and yields a determinant of 0.2858. As a result, it is {(4.2409/0.2858)^1/3-1}100% = 145.74% more D-efficient.

5. Optimal designs for nonlinear regression models

(11)

Here, Y_i represents the concentration of a substance B in two consecutive first order chemical reactions

(12)

where and represent the rates of conversion, and t_i is the time at which Y_i is observed. An example of two consecutive first order chemical reactions is the transition of O₂ to H₂O via H₂O₂. Another example of a process of consecutive first-order reactions is the nuclear decay of ₉₂U²³⁸ which eventually leads to ₈₂Pb²⁰⁶.

Suppose, for example, that four observations are available to estimate the unknown parameters and in (11). The design matrix can then be written as

(13)

The elements of the extended design matrix F are now also given by

where

(14)

and

.

(15)

These expressions show that the extended design matrix F depends on the unknown parameters and . As a consequence, the experimenter needs prior guesses of the parameter estimates. Suppose, as in Atkinson and Hunter (1968), that the experimenter believes that is around 0.7 and that is around 0.2. In that case, the derivatives become

(16)

and

.

(17)

It can be verified using a spreadsheet that the determinant of (F^TF)^-1 equals 1.024 if we would observe the concentration of substance B at time t₁=1 in the first experimental run, at time t₂=2 in the second run, and at times t₃=3 and t₄=4 in the third and fourth run respectively. By using Solver, a design that produces a determinant of 0.3807 can be found. This design had two runs in which the concentration of substance B is observed at time t=1.2295 and two runs in which the response is observed at time t=6.8577.

This design is D-optimal for the parameter values =0.7 and =0.2. It is easy to verify how sensitive the design is to the prior guess of and . Choosing =0.8 and =0.1 leads to a design with two observations at time t=1.1467 and two observations at time t=11.2963. Choosing =0.6 and =0.3 leads to a design with two observations at time t=1.3042 and two observations at time t=5.7254.

6. Discussion

In small groups (about 10 to 15 students), the students quickly learned by hands-on experience (in a pc-room) to solve example problems similar to those described in this paper. In the cases with multiple local optima, they each started with different starting design points and exchanged their solutions so as to find the global optimum. In doing so, they grasped that it can be difficult to find a globally optimal design and why the optimal solution is better than alternative solutions. These insights are only possible because the spreadsheet approach is not a 'black-box' method: every step in the solution is visible in the spreadsheet and has to be carried out by the students themselves.

The spreadsheet approach can also be used to illustrate a number of advanced design issues, such as the effect of including center points on design efficiency and the principle of design augmentation. In addition, the use of design criteria other than the D-optimality criterion can be demonstrated, as well as the design of experiments in the presence of heteroscedastic or correlated errors. For example, the design of small blocked and split-plot experiments can be tackled as well. The optimal design of experiments involving qualitative experimental variables can be illustrated as well. Unlike with quantitative variables, standard Solver is however often unable to find the optimal design for these kind of problems.

The optimization software that is readily available in spreadsheets is however not a panacea to solve real-life optimal design problems. For more complex and larger design problems, the use of specialized software such as Design Expert, Genstat, JMP or SAS proc optex is needed. The purpose of the approach presented in this paper was mainly to familiarize students with the ideas of optimal design and with the large number of applications in which they can be employed.

Appendix A

Solver allows for the optimization of an objective function, subject to a number of constraints. Standard Solver is an add-in in Excel. If it has not yet been activated, it can be done by selecting the Add-Ins... item in the Tools menu.

Figure 2

Figure 2. Activating Solver in Excel.

After checking the Solver Add-in box (see Figure 2), Solver can be opened from the Tools menu, through the item Solver... A window as in Figure 3 appears. The parameters the user should provide to Solver are:

• The target cell, i.e., the cell containing the optimality criterion.

• The adjustable cells. These are the cells containing the design points.

• The constraints, i.e., the upper and lower values for the design points.

Figure 3

Figure 3. The Solver window.

The constraints are entered for each variable separately by clicking the Add button, which makes the window in Figure 4 appear. The easiest way to enter the constraints is by using (named) ranges of cells.

Figure 4

Figure 4. The Add Constraint window.

As indicated in Section 3.1, it may be necessary to change the precision and convergence values used by Solver. This can be done by clicking the Options button in Figure 3. When the parameters have been entered, the optimization problem is solved by hitting the Solve button. If Solver finds a solution, the original values of the independent variables are automatically replaced by the optimal values. Through a communication window (as in Figure 5), the user may keep the Solver solution or may choose to restore the original values of the design points. The solution can be saved for later use (Save scenario). This might be useful for cases where Solver returns a local minimum, as in Section 3. In such cases, different starting design points might yield different solutions, which should be compared in order to find the global minimum. Each of the saved scenarios can later be shown through the Scenario... item in the Tools menu.

Figure 5

Figure 5. The Solver Results window.

Solver generates reports on demand. For each report that is selected in the Solver Results window, a worksheet is added to the workbook. The Answer report lists the target cell and the adjustable cells with their original and final values, constraints, and information about the constraints. The Sensitivity report provides information about how sensitive the solution is to small changes in the objective function or the constraints. Finally, the Limits report lists the target cell and the adjustable cells with their respective values, lower and upper limits, and target values. The Sensitivity and Limits reports are not generated for models with integer constraints.

Appendix B

In order to compute the transpose of the extended design matrix, one option is to select a 2 x 6 matrix in the worksheet, for example the cells A8:F9, to type the function

=transpose(A1:B6)

and to conclude the operation by pressing Ctrl-Shift-Enter. This way, the cells A8:F9 are defined as a matrix and can no longer be modified cell by cell.

The next step in the computation of the D-criterion value is the computation of the matrix product F^TF. To do so, select a 2 x 2 square, for example the cells A11:B12, type

=mmult(A8:F9,A1:B6)

and conclude by pressing Ctrl-Shift-Enter. Next, F^TF is inverted by selecting another 2 x 2 square, for example the cells A14:B15, entering the function

=minverse(A11:B12)

and pressing Ctrl-Shift-Enter. Finally, the determinant is obtained by entering

=mdeterm(A14:B15)

in any available cell.

Acknowledgments

References

Atkinson, A.C., and Donev, A.N. (1992), Optimum Experimental Designs, Oxford U.K.: Clarendon Press.

Atkinson, A.C. and Hunter, W.G. (1968), "The design of experiments for parameter estimation," Technometrics 10, 271-288.

Fedorov, V.V. (1972), Theory of Optimal Experiments, New York: Academic Press.

Kiefer, J., and Wolfowitz, J. (1959), "Optimum designs in regression problems," Annals of Mathematical Statistics, 30, 271-294.

Kuehl, R.O. (2000), Design of Experiments: Statistical Principles of Research Design and Analysis, New York: Duxbury.

Meyer, R.K., and Nachtsheim, C.J. (1995), "The coordinate-exchange algorithm for constructing exact optimal experimental designs," Technometrics, 37, 60-69.

Montgomery, D.C. (2000), Design and Analysis of Experiments, New York: Wiley.

Neter, J., Kutner, M.H., Nachtsheim, C.J., and Wasserman, W. (1996), Applied Linear Statistical Models, London: Irwin.

Oehlert, G.W. (2000), A First Course in Design and Analysis of Experiments, New York: Kluwer Academic Publishers.

Savova, I., Donev, T.N., Tepavicharova, I., and Alexandrova, T. (1989), "Comparative studies on the storage of freeze-dried yeast strains on the genus saccharomyces," Proceedings of the 4th International School on Cryobiology and Freeze-drying, 29 July - 6 August 1989, Borovets, Bulgaria, pp. 32-33.

Weber, D.C., and Skillings, J.H. (1999), A First Course in the Design of Experiments: A Linear Models Approach, New York: CRC Press.

Peter Goos
Department of Mathematics, Statistics and Actuarial Science
University of Antwerp
Prinsstraat 13
2000 Antwerpen
Belgium
Peter.Goos@ua.ac.be

Herlinde Leemans
DOC - Research Coordination Office
Katholieke Universiteit Leuven
Naamsestraat 22,
B-3000 Leuven
Belgium
Herlinde.Leemans@doc.kuleuven.ac.be