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), www.amstat.org/publications/jse/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

In this paper, we present an interactive teaching approach to introduce the concept of optimal design of experiments to students. Our approach is based on the use of spreadsheets. One advantage of this approach is that no complex mathematical theory is needed nor that any design construction algorithm has to be discussed at the introductory stage. Another benefit is that the students build all necessary matrices for concrete examples starting from a sensible initial design. By modifying the initial design by trial and error, they can try to improve the properties of the parameter estimators interactively. For problems in which finding the optimal design is not evident, they can use optimization software which is readily available in the spreadsheet software.

1. Introduction

The optimal design approach is a powerful and flexible way to generate efficient experimental designs. It allows the researcher to compute a tailor-made design with any number of observations, any number of experimental variables and any type of variables (qualitative or quantitative variables) for fitting any type of model (e.g., linear or nonlinear in the unknown parameters). Moreover, it is easy to take into account any constraints that are imposed on the experimental variables so that the problem of designing mixture experiments can be tackled as well.

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

Physicians often use the so-called diving reflex to reduce abnormally rapid heartbeats in humans by submerging the patients' faces in cold water. Suppose that a research physician would like to conduct an experiment in order to investigate the effects of the water temperature on the pulse rates of six small children. One intuitive way to approach the problem is to select six temperatures and to assign each of the children in a random fashion to these temperatures. A reasonable set of temperatures for this problem might be 45, 50, 55, 60, 65 and 70oF. The matrix

(1)

where xi denotes the water temperature for the ith 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 Yi denotes the reduction in pulse for the ith 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 (FTF)-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 (FTF)-1. For the six selected temperatures, the determinant of (FTF)-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 FT.

  4. Calculate (FTF).

  5. Calculate (FTF)-1.

  6. Calculate the design criterion value, e.g. the D-criterion value, i.e. the determinant of (FTF)-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

In this section, we illustrate the optimal design of experiments for multiple linear regression by means of two examples. First, the case of a quadratic regression on a single explanatory variable will be considered. Next, a quadratic regression model in two variables will receive attention.

3.1 Quadratic regression on one variable

Consider an experiment involving quantum-well lasers to study the impact of waveguide aluminum (Al) mole fraction (x) on solar pumping threshold current (y, measured in Ampère/cm2). This type of experiment is used to study the viability of using semiconductor lasers for solar lighting in spaceborne applications. Typically, the data show that the relationship between x and y is not linear but quadratic. The model under investigation can then be written as

(8)

where Yi denotes the threshold current observed at the ith trial and xi 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 (FTF) 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

A nice example of an experiment involving two independent variables is described in Savova, Donev, Tepavicharova and Alexandrova (1989). In the experiment, nine runs were performed to investigate the impact of the amount of glycerine used (x1, measured in %) and the speed of reducing the temperature (x2, expressed in oF/minute) on the amount of surviving biological material (Y, measured in %) after freeze-drying. The amount of surviving biological material is essential for the preservation of the taste and the smell of a freeze-dried product. The amount of glycerine used in the experiment ranged from 10% to 30% and the speed of the temperature reductions lay between 10oF/minute and 30oF/minute. The model estimated was

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 32 factorial design, which is a design with observations at the factor level combinations (10%,10oF), (10%,20oF), (10%,30oF), (20%,10oF), (20%,20oF), (20%,30oF), (30%,10oF), (30%,20oF) and (30%,30oF). 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 spreadsheet approach can also be used to demonstrate the usefulness of optimal designs when constraints are imposed on the levels of the experimental variables. Consider, for example, the constrained mixture experiment for estimating the impact of three factors on the electric resistivity (Y) of a modified acrylonitrile powder described in Atkinson and Donev (1992). The components of the mixture under investigation were

x1 copper sulphate (CuSO4),
x2 sodium thiosulphate (Na2S2O3),
x3 glyoxal (CHO)2


Table 1. Components of the mixture

Starting design D-Optimal design
x1 x2 x3 x1 x2 x3
1 0.20 0.40 0.40 0.20 0.20 0.60
2 0.20 0.60 0.20 0.20 0.20 0.60
3 0.30 0.35 0.35 0.20 0.20 0.60
4 0.40 0.20 0.40 0.20 0.80 0.00
5 0.40 0.60 0.00 0.20 0.80 0.00
6 0.45 0.45 0.10 0.20 0.80 0.00
7 0.50 0.25 0.25 0.80 0.20 0.00
8 0.60 0.20 0.20 0.80 0.20 0.00
9 0.60 0.40 0.00 0.80 0.20 0.00


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 xi 1. In addition, the proportions in a mixture experiment have to add up to 100%, so that

x1 + x2 + x3 = 1.

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

Optimally designing an experiment for a model that is nonlinear in the unknown parameters is more complicated than for a linear model. This is because, for a nonlinear model, the variance-covariance matrix of the estimators and thus also the optimal design depends on the values of the unknown model parameters. Nevertheless, optimal designs can be computed if prior guesses of the unknown parameters are available. Consider, for example, the equation

(11)

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

(12)

where and represent the rates of conversion, and ti is the time at which Yi is observed. An example of two consecutive first order chemical reactions is the transition of O2 to H2O via H2O2. Another example of a process of consecutive first-order reactions is the nuclear decay of 92U238 which eventually leads to 82Pb206.

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 (FTF)-1 equals 1.024 if we would observe the concentration of substance B at time t1=1 in the first experimental run, at time t2=2 in the second run, and at times t3=3 and t4=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

Our spreadsheet approach for teaching optimal design has been used during a design of experiments course in the master of statistics program, during in-company courses on optimal design and during a doctoral course on design of experiments for economics students. As in the paper, simple but appealing examples were used in the courses to introduce the concept of optimal design and to show how an optimally designed experiment can be derived. An important feature of the teaching approach presented here is that it does not require students to study design construction algorithms as those of Fedorov (1972) and Meyer and Nachtsheim (1995), for example, nor does it require the use of specialized software like Design Expert, Genstat, JMP or SAS proc optex during the classes. Instead, the students can use spreadsheet software they are familiar with, for example Microsoft Excel. Another advantage of the approach presented here is its focus on applications instead of mathematical theory involving measures, continuous designs and the general equivalence theorem (see Kiefer and Wolfowitz (1959)).

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

We have used Microsoft Excel with standard Excel Solver to optimize the example problems in this paper. Standard Excel Solver is included in Excel. It is a product of Frontline Systems, Inc., a company that developed spreadsheet solver and optimization software for Excel, as well as for Lotus 1-2-3 and Quattro Pro. The applicability of this paper is therefore not limited to Excel users only.

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:



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 illustrate the use of Microsoft Excel for matrix computations, we start from the 6 x 2 extended design matrix for the example in Section 2. We assume that this matrix is in the cells A1:B6 of the worksheet.

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 FTF. 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, FTF 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

This article was finished while the first author was a member of the Department of Applied Economics of the Katholieke Universiteit Leuven as a post-doctoral fellow of the Fund for Scientific Research - Flanders (Belgium). The authors would like to thank the two anonymous referees for their valuable comments.


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


Volume 12 (2004) | Archive | Index | Data Archive | Information Service | Editorial Board | Guidelines for Authors | Guidelines for Data Contributors | Home Page | Contact JSE | ASA Publications