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

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.

(1) |

where *x _{i}* denotes the water temperature for the

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

(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**^{T}**F**)^{-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**^{T}**F**)^{-1}. For the six selected temperatures, the
determinant of (**F**^{T}**F**)^{-1} amounts to 0.000381. This
can be verified by using a spreadsheet, constructed following the
steps below:

- 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. - 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. - Transpose the extended design matrix to obtain
**F**^{T}. - Calculate (
**F**^{T}**F**). - Calculate (
**F**^{T}**F**)^{-1}. - Calculate the design criterion value, e.g. the D-criterion value, i.e.
the determinant of (
**F**^{T}**F**)^{-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.

(8) |

where *Y _{i}* denotes the threshold current observed at the

. | (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**^{T}**F**) 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).

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^{2} factorial
design, which is a design with observations at the factor level
combinations (10%,10^{o}F), (10%,20^{o}F),
(10%,30^{o}F), (20%,10^{o}F),
(20%,20^{o}F), (20%,30^{o}F),
(30%,10^{o}F), (30%,20^{o}F) and
(30%,30^{o}F). 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).)

x_{1} | copper sulphate (CuSO_{4}), |

x_{2} | sodium thiosulphate (Na_{2}S_{2}O_{3}), |

x_{3} | glyoxal (CHO)_{2} |

Starting design | D-Optimal design | |||||

x_{1} | x_{2} | x_{3} |
x_{1} | x_{2} | x_{3} | |

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 *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 | x_{1} | 0.8, |

0.2 | x_{2} | 0.8, |

0.0 | x_{3} | 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.

(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_{2} to H_{2}O via
H_{2}O_{2}. Another example of a process of
consecutive first-order reactions is the nuclear decay of
_{92}U^{238} which eventually leads to _{82}Pb^{206}.

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**^{T}**F**)^{-1} equals 1.024 if we would observe the
concentration of substance B at time *t*_{1}=1 in the first
experimental run, at time *t*_{2}=2 in the second run, and at times
*t*_{3}=3 and *t*_{4}=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.

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.

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.

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**^{T}**F**. 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**^{T}**F** 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.

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