Statistical Experimental Design: Designs for Response Surfaces

Figure 1: Initial runs are using a factorial design, while later runs use a CCD.

We have been fitting surfaces to find some configuration of factors that optimizes a response. Sometimes we found the optimum just from one experiment, but more often, we ran a series of experiments, gradually refining our understanding of the response surface.

For any of these experiments, what designs should we use?

In theory, any of the designs we have discussed in this course could be used, but we will prefer those which give us an accurate estimate of the response surface (especially around stationary points) using as few samples as possible.

Designs for First-Order Models

In first-order models, we are interested in linear effects associated with each factor. Reasonable choices in this setting are,

$2^{K}$ factorial designs
$2^{K - p}$ fractional factorial designs
$2^{K}$ and $2^{K - p}$ designs that have been augmented with center points. These center points allow estimation of measurement error $\sigma^2$ even when there is no replication.

Designs for Second-Order Models

For second-order models, we need to estimation nonlinear and interaction effects, which means we need more intensive sampling. The most common choice is the,

Central Composite Design (CCD): As discussed in Section 6.8, this is a full factorial $2^{K}$ design that has been supplemented by center and axial points.

though some alternatives are,

Box-Behnken design: The associated samples are all on edges of the cube (none at vertices or in the interior).
Equiradial design: This generalizes the simplex idea above. Instead of sampling at the corners of an equilateral triangle, we can sample at corners or regular polygons (+ center points)

In practice, you can use functions from the rsm package. This code produces a central composite design with three center points for both the factorial and axial runs. The package also has code for factorial, fractional factorial, and Box-Behnken designs.

library(rsm)

codings <- list(
  time_coded ~ (time - 35) / 5,
  temp_coded ~ (temp - 150) / 5
)
ccd_design <- ccd(~ time_coded + temp_coded, coding = codings, oneblock = TRUE)
head(ccd_design)

  run.order std.order     time     temp
1         1         3 35.00000 142.9289
2         2         1 30.00000 145.0000
3         3         6 35.00000 150.0000
4         4         1 27.92893 150.0000
5         5         8 35.00000 150.0000
6         6         2 42.07107 150.0000

Data are stored in coded form using these coding formulas ...
time_coded ~ (time - 35)/5
temp_coded ~ (temp - 150)/5

ggplot(decode.data(ccd_design)) +
  geom_point(
    aes(x = time, y = temp),
    position = position_jitter(w = 0.5, h = 0.5),
    size = 3
  ) +
  coord_fixed()

Figure 2: A generated CCD using the rsm package.

Using the same function, we can easily generate fractional or foldover designs, which are useful during the first-order phase of response surface modeling. For example, the code below generates a $2^{4 - 1}$ design.

codings_x = list(
  x1 ~ (x1_ - 10) / 3,
  x2 ~ (x2_ - 20) / 6,
  x3 ~ (x3_ + 1) / 4
)
ffactorial <- cube(
  ~ x1 + x2 + x3, 
  n0 = 0, 
  coding = codings_x,
  generators = x4 ~ x1 * x2 * x3
  )

And if we want to foldover this fractional factorial, we can use foldover.

foldover(ffactorial, randomize = FALSE)

mfactorial <- ffactorial %>%
  head() %>%
  melt(id.vars = "std.order") %>%
  filter(variable != "run.order")
ggplot(mfactorial) +
  geom_tile(
    aes(x = variable, y = as.factor(std.order), fill = as.factor(value))
  ) +
  coord_fixed() +
  labs(x = "factor", y = "run", fill = "level") +
  scale_fill_brewer(palette = "Set2")

$A $2^{4 - 1}$ fractional factorial design, made using the `rsm` package.$

Figure 3: A $2^{4 - 1}$ fractional factorial design, made using the rsm package.

mfactorial <- foldover(ffactorial, randomize=FALSE) %>%
  head() %>%
  melt(id.vars = "std.order") %>%
  filter(variable != "run.order")
ggplot(mfactorial) +
  geom_tile(
    aes(x = variable, y = as.factor(std.order), fill = as.factor(value))
  ) +
  coord_fixed() +
  labs(x = "factor", y = "run", fill = "level") +
  scale_fill_brewer(palette = "Set2")

$The foldover of the previous fractional factorial design.$

Figure 4: The foldover of the previous fractional factorial design.

Variance of Response Surfaces

At any point $x$, it’s possible to estimate the variance of the surface at that point, denoted $V\left(\hat{y}\left(x\right)\right)$.

For any response surface, we can also imagine an accompanying variance surface.
Away from factor configurations that we’ve sampled, we expect the variance to increase.

If there are only 2 factors, this surface can be directly visualized.

But in higher-dimensions, this is impossible.
In that case, use a variance dispersion graph to visualize the variance as a function of the distance from the center point of the design.
- Specifically, plot the maximum, minimum, and average variances along rings that emanate from 0
If the variances are constant along each ring, the design is called rotatable

ccd_design <- ccd(~ time_coded + temp_coded, coding = codings, oneblock = TRUE)
par(mfrow = c(1, 2))
varfcn(ccd_design, ~ SO(time_coded, temp_coded))
varfcn(ccd_design, ~ SO(time_coded, temp_coded), contour = TRUE, asp = 1)

Figure 5: The scaled prediction variance functions for a rotatable CCD.

ccd_design <- ccd(~ time_coded + temp_coded, coding = codings, alpha=2, oneblock = TRUE)
par(mfrow = c(1, 2))
varfcn(ccd_design, ~ SO(time_coded, temp_coded))
varfcn(ccd_design, ~ SO(time_coded, temp_coded), contour = TRUE, asp = 1)

Figure 6: The analogous functions for a nonrotatable CCD. Note that the axis and diagonal variance curves don’t overlap.