APTS Design of Experiments 2022

Optimal design via algorithmic search

library(apts.doe)
library(AlgDesign)
library(rsm)
library(dplyr)
library(plotly)

Point exchange algorithm

• Most common optimisation algorithm for optimal design
• Swaps between current design and a candidate list of possible design points
• Can be applied with most optimality criteria, including D-optimality
• Implemented in package AlgDesign

We will study three examples of using AlgDesign to find designs under linear models \[ Y = X\beta + \epsilon\,, \] with \(X\) the model matrix, \(\beta\) the \(p\)-vector of unknown model coefficients and \(\epsilon\) an \(n\)-vector of identically and independently normally distributed random errors, each having variance \(\sigma^2\).

A \(D\)-optimal design maximises the determinant of the information matrix \(X^T X\).

Example 1: one-factor quadratic regression

D-optimal design for the regression model

\[ Y_i = \beta_0 + \beta_1 x_i + \beta_2 x_i^2 + \epsilon_i \]

dopt1

$D
[1] 0.5088144

$A
[1] 3.139587

$Ge
[1] 0.941

$Dea
[1] 0.939

$design

$rows
[1]   1   2   3  50  51  52  99 100 101

The function returns

• D: the value of the determinant of normalised information matrix (divided by \(n\)) raised to the power \(1/p\)
• A: the value of the trace of the inverse of the information matrix, divided by \(p\)
• Ge and Dea: the \(G\)- and \(D\)-efficiencies relative to the theoretical optimal design (which may not be realisable)
• rows: the rows of the candidate list that correspond to the optimal design.

Example 2: a response surface design

A pharmaceutical industry experiment (Owen et al., 2001, Organic Process Research and Development) investigated four factors in \(n=30\) runs. We will find a \(V\)-optimal design for this experiment, assuming

\[ -2 \le x_1, x_2, x_3, x_4 \le 2\,. \]

Note that AlgDesign refers to \(V\)-optimality as \(I\)-optimality (\(I\) for integrated variance).

cand.list <- expand.grid(x1 = seq(-2, 2, length = 5), x2 = seq(-2, 2, length = 5), 
                         x3 = seq(-2, 2, length = 5), x4 = seq(-2, 2, length = 5))
dopt2 <- optFederov(~ quad(x1, x2, x3, x4), data = cand.list, nTrials = 30, 
                    criterion = "I")
dopt2

$D
[1] 3.952834

$A
[1] 0.7326728

$I
[1] 11.25339

$Ge
[1] 0.585

$Dea
[1] 0.492

$design

$rows
 [1]   1  13  25  55  63  71 101 103 115 121 250 255 273 303 311 313 315 361 373 396 476 503 511 515 551 573 575
[28] 605 613 621

It is interesting to compare the performance of this design under the various optimality criteria with the original design used in the study. This was a central composite design, or CCD, a common multi-factor design for fitting regression models. To do this we generate a CCD using the rsm package.

ccd_design

   run.order std.order x1.as.is x2.as.is x3.as.is x4.as.is
1          1         1       -1       -1       -1       -1
2          2         2        1       -1       -1       -1
3          3         3       -1        1       -1       -1
4          4         4        1        1       -1       -1
5          5         5       -1       -1        1       -1
6          6         6        1       -1        1       -1
7          7         7       -1        1        1       -1
8          8         8        1        1        1       -1
9          9         9       -1       -1       -1        1
10        10        10        1       -1       -1        1
11        11        11       -1        1       -1        1
12        12        12        1        1       -1        1
13        13        13       -1       -1        1        1
14        14        14        1       -1        1        1
15        15        15       -1        1        1        1
16        16        16        1        1        1        1
17        17        17        0        0        0        0
18        18        18        0        0        0        0
19        19        19        0        0        0        0
20        20        20        0        0        0        0
21        21        21        0        0        0        0
22        22        22        0        0        0        0
23         1         1       -2        0        0        0
24         2         2        2        0        0        0
25         3         3        0       -2        0        0
26         4         4        0        2        0        0
27         5         5        0        0       -2        0
28         6         6        0        0        2        0
29         7         7        0        0        0       -2
30         8         8        0        0        0        2

Data are stored in coded form using these coding formulas ...
x1 ~ x1.as.is
x2 ~ x2.as.is
x3 ~ x3.as.is
x4 ~ x4.as.is

Example 3: a multi-factor experiment with dependent factor settings

McNamara et al. (2005, Tech. Report) described an experiment to investigate the relationship between chemical yield and choice of solvent.

• 13 solvents available for the experiment
• characterised by values of three chemical descriptors: dielectric constant, dipole moment and water solubility
• Replacing a 13-level qualiative factor (solvent) with three quantitative variables
• combinations of the descriptors which are feasible for the experiment are determined by the list of available solvents.

solvents <- read.csv("solvents.csv", header = T)
solvents

\(D\)-optimal design to estimate a quadratic regression model with \(n=20\)

• the design problem is really picking the replication of each solvent
• optFederov requires the candidate list to have more rows than \(n\), so we use three replicates of the solvent list

dopt3 <- optFederov(~ quad(dc, dp, ws), data = rbind(solvents, solvents, solvents), 
                    nTrials = 20, criterion = "D") 
with(dopt3$design, {
  dopt3$design[order(solvent), ]
})

Notice that not all of the solvents have been selected in the design (so some have replication 0).

Compare this design to the \(D\)-optimal design that would be found if we assumed all possible combinations of descriptor values were possible (or a least a much wider range). That is, what design would be chosen if we were not restricted to our list of 13 solvents?

descrip.list <- expand.grid(dc = seq(-1, 1, len = 9), dp = seq(-1, 1, len = 9), 
                            ws = seq(-1, 1, len = 9))
dopt3b <- optFederov(~ quad(dc, dp, ws), data = descrip.list, nTrials = 20, 
                    criterion = "D")
with(dopt3b$design, {
  dopt3b$design[order(dc, dp, ws), ]  
})

This design has only one point in the interior of the design space (at the centre point), unlike the design constructed from the original candidate list. It has no replication.

dopt3b$design %>% count(dc, dp, ws)

round(100 * (dopt3$D/dopt3b$D) ^ (1 / 10))

[1] 82

Original design has low efficiency but is realisable (the solvents exist!), whereas we do not know if the combinations of the descriptors in the second design correspond to actual solvents.

plot_ly(dopt3$design %>% count(dc, dp, ws), x = ~dc, y = ~dp, z = ~ws, type = "scatter3d", mode = "markers", color = ~n)

plot_ly(dopt3b$design, x = ~dc, y = ~dp, z = ~ws, type = "scatter3d", mode = "markers")