5.17. Judging and comparing experimental designs

By this point we have several ways to build a design: full factorials, fractional factorials, central composite designs, and the more flexible optimal designs. A practical question follows almost immediately: when you are handed two or three candidate designs (perhaps a small screening design, a slightly larger one, and a classical response surface design), how do you decide which is better?

A good design has to do two separate jobs well, and it helps to keep them apart in your mind:

1. Separability. Can we tell the effects apart, so that the estimate of one term is not confounded with another? This is the language of aliasing and resolution.

2. Precision. Even when the effects are separable, how precisely can we estimate the coefficients and predict the response?

Both questions are answered by a single object built from the design: the information matrix \(\mathbf{M} = \mathbf{X}^T\mathbf{X}\), together with the optimality criteria that summarise it. Both were introduced in the earlier section on optimal designs, so here we take them as given. This subchapter shows how prediction variance is derived from \(\mathbf{M}\), how to read a fraction-of-design-space plot, and how the separability, bias, and power measures play out on one running comparison of two designs. The companion page then widens that comparison across the standard design families and closes with a checklist for choosing among them.

5.17.1. Prediction variance

We rarely stop at the coefficients; in response surface work we use the model to predict, and to chase an optimum, so the variance of a prediction is what we ultimately care about. The predicted response at any point \(\mathbf{x}\) is \(\widehat{y}(\mathbf{x}) = \mathbf{x}_m(\mathbf{x})^T \mathbf{b}\), where \(\mathbf{x}_m(\mathbf{x})\) is the model expansion of that point. Propagating the variance through this linear combination, using the rule \(\text{Var}(\mathbf{a}^T\mathbf{b}) = \mathbf{a}^T\,\text{Var}(\mathbf{b})\,\mathbf{a}\) with \(\mathbf{a} = \mathbf{x}_m(\mathbf{x})\),

\[\text{Var}\big(\widehat{y}(\mathbf{x})\big) = \mathbf{x}_m(\mathbf{x})^T \,\text{Var}(\mathbf{b})\, \mathbf{x}_m(\mathbf{x}) = \sigma^2\, \mathbf{x}_m(\mathbf{x})^T \mathbf{M}^{-1} \mathbf{x}_m(\mathbf{x})\]

Read this as a product of three things. The variance of a prediction at a new point depends on where you are predicting (the model expansion \(\mathbf{x}_m(\mathbf{x})\), which encodes how far that point sits from the centre of the design in coded space), on how the design is laid out (the inverse information matrix \(\mathbf{M}^{-1}\), the entire geometry of the runs distilled into one matrix), and on how noisy the system is (the baseline variance \(\sigma^2\)). Only the middle term is in the experimenter’s hands: the first is set by where you happen to want a prediction, and the last is a property of the process.

That is the entire derivation: it is just error propagation through the linear predictor, using \(\text{Var}(\mathbf{b}) = \sigma^2 \mathbf{M}^{-1}\).

For the three-run example from the previous section, multiplying \([1, x, x^2]\) through \(\mathbf{M}^{-1}\) and contracting gives a tidy polynomial:

\[\text{Var}\big(\widehat{y}(x)\big) = \sigma^2 \left(1 - 1.5\,x^2 + 1.5\,x^4\right)\]

At the three design points (\(x = -1, 0, +1\)) this equals \(\sigma^2\), as it must: a saturated design interpolates its own data. Between the points it dips to a minimum of \(0.625\,\sigma^2\) at \(x = \pm 0.707\), and beyond \(x = 1\) it climbs steeply (already \(5.2\,\sigma^2\) at \(x = 1.5\)): a quantitative warning against extrapolation.

This and every figure in this subchapter is reproducible with process_improve (pip install 'process-improve[all]', which installs every optional extra). Of those extras, expt (the pyDOE3 package) is needed for the Box-Behnken and central composite designs on the companion page, and ilp (the pulp integer-programming solver) for the twenty-five-run OMARS design there. The four-factor, thirteen-run OMARS design used on this page is small enough for generate_omars to find by exhaustive enumeration (every feasible design of that size is listed and scored), so it needs neither. Each block imports what it needs and reuses variables defined in the blocks before it, so paste them in order. The prediction variance of the three-run quadratic design is a closed form:

Python, 12 lines

import numpy as np
import plotly.graph_objects as go

# Single-factor quadratic design {-1, 0, +1}: the prediction variance is a closed form.
x = np.linspace(-1.6, 1.6, 321)
pred_var = 1 - 1.5 * x**2 + 1.5 * x**4              # Var(y_hat) / sigma^2

fig = go.Figure(go.Scatter(x=x, y=pred_var, mode="lines"))
fig.add_vrect(x0=-1.6, x1=-1, fillcolor="LightSalmon", opacity=0.2, line_width=0)
fig.add_vrect(x0=1, x1=1.6, fillcolor="LightSalmon", opacity=0.2, line_width=0)
fig.update_layout(xaxis_title="x (coded)", yaxis_title="Prediction variance / sigma^2")
fig.show()

Prediction variance for the three-run quadratic design. It equals \(\sigma^2\) at the three design points, dips to \(0.625\,\sigma^2\) midway between them, and climbs steeply once \(x\) leaves the design region \([-1, +1]\). The shaded bands are the extrapolation zone, where predictions become rapidly less certain.

To compare designs we strip out two nuisance factors. We divide by the unknown \(\sigma^2\) (a property of the process, not the design), and we multiply by the number of runs \(N\) (otherwise a design looks better merely for being larger: replicate any design and \(\mathbf{M}\) doubles, halving the variance). The result is the scaled prediction variance:

\[\text{SPV}(\mathbf{x}) = N\, \mathbf{x}_m(\mathbf{x})^T \mathbf{M}^{-1} \mathbf{x}_m(\mathbf{x}) = \frac{N\,\text{Var}\big(\widehat{y}(\mathbf{x})\big)}{\sigma^2}\]

The SPV depends only on the geometry of the design. The G-optimal value is its maximum over the region, and the I-optimal value is its average over the region (V-optimality is the same average taken over a chosen set of points rather than over the whole region).

5.17.2. A worked example: augmenting a small design

The previous section worked out the information matrix of the smallest design with curvature: a single factor at three levels, \(x = -1, 0, +1\), fitting the quadratic model \(y = b_0 + b_1 x + b_2 x^2 + e\). That single design is rarely the end of the story. Suppose the budget stretches to two more runs. Two natural options present themselves, and they pull in different directions:

• add two replicate centre points (two more runs at \(x = 0\)), or

• add a point at \(x = -1\) and one at \(x = +1\), reinforcing the extremes.

Which is better? It depends entirely on what we ask of the design, and the optimality criteria turn that vague question into an arithmetic one. The table below evaluates four designs on the same quadratic model: the base three-run design, the base design with all three runs repeated, the base plus two centre points, and the base plus two extreme points. The criteria are the raw, unnormalised summaries of \(\mathbf{M}\) from earlier (\(D = \det\mathbf{M}\), \(A = \text{trace}\,\mathbf{M}^{-1}\), \(E = \lambda_{\min}(\mathbf{M})\)), together with the maximum and average of the prediction variance \(d(x) = \mathbf{x}_m^T \mathbf{M}^{-1} \mathbf{x}_m\) over \(x \in [-1, +1]\) (these are the \(G\) and \(I\) quantities, in units of \(\sigma^2\)).

Four candidate designs for the single-factor quadratic model.

Design	\(N\)	\(\uparrow\ D=\det\mathbf{M}\)	\(\downarrow\ A=\text{trace}\,\mathbf{M}^{-1}\)	\(\uparrow\ E=\lambda_{\min}\)	\(\downarrow\ G\)	\(\downarrow\ I\)
Base \(\{-1, 0, +1\}\)	3	4	3.00	0.44	1.0	0.80
Base, all three runs repeated	6	32	1.50	0.88	0.5	0.40
Base + two centre points	5	12	1.67	1.00	1.0	0.44
Base + two points at \(\pm 1\)	5	16	2.50	0.47	1.0	0.67

It is worth seeing the matrices these summaries come from. Writing each run as \(\mathbf{x}_m = [\,1,\ x,\ x^2\,]\) and stacking the runs as the rows of \(\mathbf{X}\), the base three-run design (row 1) gives the matrix already worked out in the previous section,

\[\begin{split}\mathbf{X}_1 = \begin{bmatrix} 1 & -1 & 1 \\ 1 & 0 & 0 \\ 1 & 1 & 1 \end{bmatrix}, \qquad \mathbf{M}_1 = \mathbf{X}_1^T\mathbf{X}_1 = \begin{bmatrix} 3 & 0 & 2 \\ 0 & 2 & 0 \\ 2 & 0 & 2 \end{bmatrix}\end{split}\]

Repeating all three runs (row 2) stacks a second copy of every row, which simply doubles the information matrix,

\[\begin{split}\mathbf{X}_2 = \begin{bmatrix} 1 & -1 & 1 \\ 1 & 0 & 0 \\ 1 & 1 & 1 \\ 1 & -1 & 1 \\ 1 & 0 & 0 \\ 1 & 1 & 1 \end{bmatrix}, \qquad \mathbf{M}_2 = \begin{bmatrix} 6 & 0 & 4 \\ 0 & 4 & 0 \\ 4 & 0 & 4 \end{bmatrix} = 2\,\mathbf{M}_1\end{split}\]

Adding two centre points instead (row 3) repeats the centre run, feeding the intercept while leaving the boundary runs unchanged,

\[\begin{split}\mathbf{X}_3 = \begin{bmatrix} 1 & -1 & 1 \\ 1 & 0 & 0 \\ 1 & 1 & 1 \\ 1 & 0 & 0 \\ 1 & 0 & 0 \end{bmatrix}, \qquad \mathbf{M}_3 = \begin{bmatrix} 5 & 0 & 2 \\ 0 & 2 & 0 \\ 2 & 0 & 2 \end{bmatrix}\end{split}\]

while adding a point at each extreme (row 4) keeps the single centre run but doubles the two boundary runs,

\[\begin{split}\mathbf{X}_4 = \begin{bmatrix} 1 & -1 & 1 \\ 1 & 0 & 0 \\ 1 & 1 & 1 \\ 1 & -1 & 1 \\ 1 & 1 & 1 \end{bmatrix}, \qquad \mathbf{M}_4 = \begin{bmatrix} 5 & 0 & 4 \\ 0 & 4 & 0 \\ 4 & 0 & 4 \end{bmatrix}\end{split}\]

The \(D\), \(A\), and \(E\) columns of the table are nothing more than the determinant, the trace of the inverse, and the smallest eigenvalue of each of these \(\mathbf{M}\) matrices, and the prediction columns \(G\) and \(I\) follow from \(\mathbf{M}^{-1}\) exactly as in the previous section. Notice already how the two five-run designs compare: \(\mathbf{M}_3\) and \(\mathbf{M}_4\) share the same intercept information (\(M_{00} = 5\), since each adds two runs to the base), but the centre-point design holds the intercept-quadratic cross-term down at \(M_{02} = 2\), whereas doubling the boundary runs drives it up to \(M_{02} = 4\). That difference in entanglement, read straight off \(\mathbf{M}_3\) and \(\mathbf{M}_4\), is the seed of the comparison below.

Before reading the numbers, fix the direction each criterion should move towards. A design is better when \(D\) is larger (a bigger determinant is more joint information and a smaller confidence ellipsoid), when \(A\) is smaller (a smaller \(\text{trace}\,\mathbf{M}^{-1}\) is a lower average coefficient variance), when \(E\) is larger (a bigger smallest eigenvalue means the worst-estimated direction is better pinned down), and when \(G\) and \(I\) are smaller (lower worst-case and lower average prediction variance). The arrow at the head of each column marks that direction of improvement: \(\uparrow\) where larger is better and \(\downarrow\) where smaller is better. “Better” points a different way in each column, which is precisely why no single design is best in every column.

Now read the rows, and notice the run-count effect first. Going from the base design to the same design with every run repeated, every criterion improves: \(D\) jumps by a factor of \(2^p\) (here \(p = 3\), so \(4 \to 32\)), \(A\) halves, \(E\) doubles, and both \(G\) and \(I\) halve. Nothing about the design got better; we only ran more experiments. This is the warning to keep for later: the raw criteria scale with the number of runs \(N\), so they cannot be used to compare designs of different size. We undo that scaling below. For the same reason, the base design cannot be compared with either five-run design.

The fair comparison is between the two five-run designs, and here the criteria start to help make a decision. Adding the two extreme points maximises \(D\) (16 versus 12 when adding two replicate centre points): spreading runs to the boundary buys more joint information on the coefficient estimates, so it is the choice when the goal is to estimate the coefficients jointly.

Adding the two centre points instead reduces \(A\) (1.67 versus 2.50) because it lowers the average coefficient variance. The reason is visible in the information matrix. The weak spot of the base design is the intercept-quadratic pair: the two are correlated (\(M_{02} = 2\)) and the quadratic is the least precise coefficient (\(\text{Var}(b_2) = 1.5\,\sigma^2\)). A centre run expands to \(\mathbf{x}_m = [\,1,\ 0,\ 0\,]\), so it adds information to the intercept alone: it raises \(M_{00}\) from 3 to 5 while leaving the entangling cross-term \(M_{02}\) untouched. Pinning the intercept down this way partly de-correlates it from \(b_2\) and pulls \(\text{Var}(b_2)\) from \(1.5\,\sigma^2\) down to \(0.83\,\sigma^2\), and that drop in the largest coefficient variance is what lowers the average.

The centre points also maximise \(E\) (1.00 versus 0.47), the smallest eigenvalue of \(\mathbf{M}\), which gauges how well the worst-estimated direction in coefficient space is pinned down. That weakest direction is the very same intercept-quadratic combination, so the runs that de-correlate the pair are exactly the ones that strengthen it. The extreme points cannot help here: their rows \([\,1,\ -1,\ 1\,]\) and \([\,1,\ +1,\ 1\,]\) deepen the entanglement (the cross-term grows to \(M_{02} = 4\)), which is why \(E\) barely moves from the base value of 0.44.

Finally the centre points minimise \(I\) (0.44 versus 0.67), the average prediction variance over the region: concentrating runs in the interior predicts better across the bulk of the factor space, where the fitted model is most often used. The worst-case prediction variance \(G\) is a tie between the two five-run designs.

So the dilemma resolves based on your intentions: reinforce the extremes if you want the tightest joint estimate of the coefficients, add centre runs if you care about average precision and prediction. That is the \(D\)-for-estimation, \(I\)-for-prediction split we return to in the closing checklist.

Return for a moment to the dilemma table. The base three-run design has \(\text{SPV}(x) = 3\,(1 - 1.5\,x^2 + 1.5\,x^4)\), with a maximum of \(3\) at the design points and an average of \(2.4\) over the region. Now scale the design that simply repeated all three runs: it has \(N = 6\) and half the prediction variance, so its SPV is \(6 \times \tfrac{1}{2}(1 - 1.5\,x^2 + 1.5\,x^4)\), the identical curve. Scaling by \(N\) has exactly cancelled the artificial gain from pure replication: on the SPV scale the two designs are correctly seen as one and the same. This is the fix promised earlier, and it is why the fraction-of-design-space plot below is built from the SPV rather than from the raw prediction variance.

That is as far as a single factor can take us. The one-factor example has supplied the working vocabulary: prediction variance, the \(D\), \(A\), \(E\), \(G\), and \(I\) criteria, and the per-run (SPV) scaling that stops sheer replication from masquerading as quality. The rest of this subchapter puts that vocabulary to work on a single running comparison of two realistic multi-factor designs, introduced in the next section.

5.17.3. A running comparison: a DSD and an OMARS design

The running example for the rest of the subchapter is a pair of four-factor designs, both fitted on the main-effects-plus-quadratic model: a nine-run definitive screening design and a thirteen-run orthogonal minimally aliased response surface (OMARS) design that spends its four extra runs to buy two estimable two-factor interactions. OMARS designs are a recent generalization of the definitive screening design: they keep the main effects orthogonal to every second-order term while trading a handful of runs for interaction estimability, and the definitive screening designs are themselves a special case within that family.

As each new metric is introduced (the FDS reading, then separability, the variance inflation factor, the alias bias, and power), it is read off these same two designs, and a single summary table collects every value at the end.

Both designs come straight from process_improve: generate_design builds the nine-run definitive screening design, and generate_omars builds the thirteen-run OMARS member directly (the DSD is the minimal member of the same foldover family), each confirmed with the library’s is_omars verifier. The fds helper defined here wraps evaluate_design, which integrates the prediction variance over the design region and returns both the FDS curve and the average and worst-case values; the next section uses it to draw the plot.

Python, 27 lines

import numpy as np
import pandas as pd
import plotly.graph_objects as go
from process_improve.experiments import (Factor, evaluate_design, generate_design,
                                         generate_omars, is_omars)

def fds(design, model, *, n_samples, seed=1):
    """Region prediction-variance summary from ``evaluate_design``: the FDS
    curve (fraction of the design space vs prediction variance, scaled and
    unscaled) together with the average and worst-case values. The 2**k cube
    vertices, where the worst case usually sits, are included by default."""
    cols = [chr(ord("A") + i) for i in range(np.shape(design)[1])]
    df = pd.DataFrame(np.asarray(design, float), columns=cols)
    return evaluate_design(df, model=model, metric="fds", n_samples=n_samples,
                           random_seed=seed, fds_resolution=200)["fds"]

# Both four-factor designs come from process_improve: the 9-run DSD and the
# precision-optimal (A-optimal) 13-run OMARS member of the same foldover family,
# the latter carrying two estimable two-factor interactions. is_omars confirms each.
factors4 = [Factor(name=c, low=-1, high=1) for c in "ABCD"]
dsd4 = np.asarray(generate_design(factors4, design_type="dsd").design[list("ABCD")],
                  float)
omars4 = np.asarray(
    generate_omars(factors4, n_runs=13, model="main_quadratic",
                   selection_criterion="a_optimal").design[list("ABCD")], float)
assert is_omars(dsd4) and is_omars(omars4)
model4 = " + ".join(list("ABCD") + [f"I({c}**2)" for c in "ABCD"])

5.17.4. The fraction-of-design-space (FDS) plot

A single design has a different SPV at every point in the factor region, so quoting one number hides a lot. The fraction-of-design-space plot (FDS plot) shows the whole distribution. It is worth being precise about how it is built, because the recipe is far less mysterious than the finished plot can look.

Start from the key fact that \(\text{SPV}(\mathbf{x}) = N\,\mathbf{x}_m(\mathbf{x})^T \mathbf{M}^{-1} \mathbf{x}_m(\mathbf{x})\) is a closed-form function of position. A “point” here is not an experimental run: it is any location \(\mathbf{x}\) in the factor region, and once the design has fixed \(\mathbf{M}\) we can evaluate the SPV there by simply plugging \(\mathbf{x}\) into the formula, with no data and no experiment required. The whole plot is a property of the design’s geometry alone.

So we sample the region. Scatter a large number of locations, tens of thousands of them (the figure below uses 80,000), uniformly at random across the coded factor region: for \(k\) factors that means drawing each coordinate independently and uniformly on \([-1, +1]\), so the points fill the cube evenly. One point is easy to miss: uniform random sampling almost never lands exactly on a corner of the cube, yet the prediction variance is usually largest at a vertex, where the G-optimal worst case sits. Add the \(2^k\) extreme vertices back to the sampled set explicitly, so the worst-case (right-hand) tail of the curve is represented and the region is covered fairly. At every one of those locations we evaluate the SPV formula. We now hold tens of thousands of SPV values, one per sampled location, a dense picture of how precisely the fitted model would predict across the entire region.

Finally we turn that cloud of numbers into a curve. Sort the SPV values from smallest to largest. The \(i\)-th value in the sorted list is plotted at horizontal position \(f = i / (\text{number of points})\), and its height is the SPV itself. The horizontal axis is therefore just the running fraction (percentiles), and the curve is the empirical cumulative distribution of SPV over the region: with enough points it converges to the true distribution.

The horizontal axis is the part that confuses people on first sight: it is not a factor axis. A point at horizontal position \(f\) means “a fraction \(f\) of the design region is predicted at least this well (this SPV or lower).” So \(f = 0.5\) is the median SPV, and \(f = 1.0\) is the single worst point in the region. Because the values are sorted, the curve only rises: the familiar gently-rising-then-steepening shape. What to read off it:

• the left end (\(f \to 0\)): the best-predicted regions, usually near the centre;

• the median (\(f = 0.5\)): the typical precision you will experience;

• the right end (\(f \to 1\)): the worst SPV in the region, which is the G-optimal value, usually at a corner;

• the overall height: lower is more precise everywhere, and the area under the curve is essentially the V-/I-optimal (average) value;

• the flatness: a flat curve means uniform precision across the whole region, which is the response-surface ideal: you predict equally well everywhere. A curve that ramps up sharply at the right predicts well in the middle but poorly at the edges.

In short, you are looking for low and flat, and when comparing two designs, for which curve sits underneath and which is flatter.

Applying that recipe to the running example reuses dsd4, omars4, model4 and the fds helper from the previous section:

Python, 8 lines

fig = go.Figure()
for design, label in [(dsd4, "DSD [n=9]"), (omars4, "OMARS [n=13]")]:
    curve = fds(design, model4, n_samples=80_000)["curve"]
    fig.add_trace(go.Scatter(x=curve["fraction"],
                             y=curve["scaled_prediction_variance"], name=label))
fig.update_layout(xaxis_title="Fraction of design space",
                  yaxis_title="Scaled prediction variance, SPV")
fig.show()

FDS plot for a nine-run definitive screening design and a thirteen-run OMARS design, both in four factors on the main-effects-plus-quadratic model. The curves cross near a fraction of 0.73: the larger design predicts better on average but worse at the extreme corners.

The thirteen-run OMARS curve sits below the nine-run DSD curve for roughly the first three quarters of the region: it has lower best-case, median, and average prediction variance. But the two curves cross near \(f \approx 0.73\), and the thirteen-run OMARS curve then rises well above: its worst-case prediction variance is noticeably higher. This crossing is the practical tension between V- (average) and G- (worst-case) optimality, and it has a physical cause: the larger design here places fewer runs out near the edge of the region, so prediction there behaves like mild extrapolation. The reading is concrete: if you care about prediction on average across the space (typical optimization work), prefer the larger design; if you must predict reliably even in the worst spot, the flatter nine-run DSD curve is the safer choice.

One detail of method is worth stating, because it is easy to get wrong. The worst-case figure \(G\) is a maximum over the whole design region, and that maximum can sit exactly at an extreme corner (a vertex of the \([-1, 1]\) cube), where random interior sampling rarely lands. evaluate_design therefore adds the cube vertices to its interior sample by default (its include_vertices argument). Including them lifts the nine-run DSD’s \(G\) from \(8.98\) to \(9.00\), a maximum that turns out to sit precisely at a corner, and leaves the thirteen-run OMARS value at \(12.50\) because its worst case does not sit at a vertex. The shift is tiny, so it changes no conclusion here, but including the extreme points is the correct procedure, and the omnibus comparison relies on it.

5.17.5. Separability is not the same as precision

The optimality criteria and the FDS plot all speak to precision. They are silent on the other question, separability, which is governed by the off-diagonal structure of \(\mathbf{M}\): how correlated the effects are with one another. The correlation between two model-term columns, \(\mathbf{c}_a\) and \(\mathbf{c}_b\), is just the cosine of the angle between them once each has been centred:

\[r_{ab} = \frac{\widetilde{\mathbf{c}}_a^{\,T} \widetilde{\mathbf{c}}_b} {\|\widetilde{\mathbf{c}}_a\| \, \|\widetilde{\mathbf{c}}_b\|}\]

where \(\widetilde{\mathbf{c}}\) is the column after regressing out the intercept and the main-effect columns, that is, the residual

\[\widetilde{\mathbf{c}} = \mathbf{c} - \mathbf{X}_0\,(\mathbf{X}_0^T \mathbf{X}_0)^{-1}\,\mathbf{X}_0^T\,\mathbf{c}, \qquad \mathbf{X}_0 = [\,\mathbf{1} \;\; \mathbf{x}_1 \;\; \cdots \;\; \mathbf{x}_k\,]\]

with \(\mathbf{X}_0\) holding the intercept column \(\mathbf{1}\) and the \(k\) main-effect columns. Two distinct pieces are removed here. Subtracting the intercept removes the column’s mean (this is the centring), which matters because the quadratic columns \(x_i^2\) have a positive mean that would otherwise inflate every correlation. Subtracting the main-effect columns then removes any linear trend the column shares with the factors. What is left is the genuine entanglement between the second-order terms. For a design whose main effects are already orthogonal to every second-order term (the definitive screening designs and their generalizations being the prime example), the main-effect part is already zero, so there only the centring does any work; in general both pieces are removed.

We report the absolute value of \(r\) because its sign is an artefact of how the factor levels happen to be coded: flip the direction of one factor and the sign of every term containing it flips with it. The magnitude is the coding-invariant quantity, and it is what governs separability. A value of \(|r| = 0\) means the two effects are orthogonal and can be estimated independently; \(|r| = 1\) means they are the same column and cannot be told apart at all. Two summaries are useful: the maximum \(|r|\) over all pairs of second-order effects (the quadratics and two-factor interactions, which is where a screening design’s entanglement lives) is the single tightest confounding anywhere in the design (the worst case you would have to defend), while the mean \(|r|\) is the overall level of entanglement. In the comparison table, the definitive screening design has a worst-pair value of \(0.707\) (a structural hallmark of DSDs, and a confounding between two of its second-order effects), which the thirteen-run OMARS design improves to \(0.570\); the mean values, \(0.322\) against \(0.307\), are much closer, telling us the extra runs help most with the worst pair rather than with the average.

It is essential to treat separability and precision as two axes, because a design can be excellent on one and poor on the other. A one-factor-at-a-time design, for instance, has almost no correlation between its effects (good separability) and yet very poor precision (its information is spread thinly, giving a low determinant and high prediction variance). Ranking designs on any single number (including a correlation summary) will eventually recommend something you would never want to run. Look at both axes.

5.17.6. Variance inflation factors

The correlation \(r\) above is a pairwise measure. Its multivariate cousin, which asks how much a coefficient suffers from its entanglement with all the other terms at once, is the variance inflation factor (VIF). For term \(j\),

\[\text{VIF}_j = \frac{1}{1 - R_j^2}\]

where \(R_j^2\) is the coefficient of determination from regressing column \(j\) of the model matrix on all of the other columns. Equivalently, \(\text{VIF}_j\) is the \(j\)-th diagonal element of the inverse of the correlation matrix of the model terms. It is the factor by which the variance of \(b_j\) is inflated relative to a perfectly orthogonal design:

\[\text{Var}(b_j) = \text{VIF}_j \cdot \frac{\sigma^2}{S_{jj}}\]

with \(S_{jj}\) the corrected sum of squares of column \(j\). In an orthogonal design every \(R_j^2 = 0\), so every \(\text{VIF}_j = 1\): the ideal. A value of 4 means the standard error of that coefficient is doubled (\(\sqrt{4}\)) by the correlation; a common rule of thumb raises a flag past 5, and a serious one past 10.

The link to \(r\) is direct: if a term were correlated with just one other term at level \(r\), its VIF would be \(1/(1 - r^2)\). The VIF generalizes this to the joint effect of every other term, and is computed on whichever model you actually intend to fit. The quadratic columns \(x_i^2\) have a nonzero mean, but because each column is centered when the factor is formed (the corrected sum of squares \(S_{jj}\) above reflects the same centering), that mean does not by itself inflate the VIF; only genuine correlation with the other fitted terms does.

In the comparison table below, both summaries are reported for the main-effects-and-quadratic model. The definitive screening design shows \(\text{VIF} = 1.0\) throughout: on that model its terms are mutually orthogonal. (This is not in tension with the worst-pair \(|r| = 0.707\) quoted above: that correlation involves the two-factor interactions, which this model leaves out. Restricted to the main effects and quadratics, the DSD really is orthogonal.) The thirteen-run OMARS design carries a maximum VIF of \(1.18\) and a mean of \(1.08\), which inflates the worst standard error by only \(\sqrt{1.18} \approx 1.09\), about nine percent. That is a mild and entirely acceptable price for the residual degrees of freedom and the interaction estimates that the extra runs provide.

5.17.7. Bias from the terms left out: the alias matrix

The variance inflation factor measures entanglement among the terms we fit. A separate question sits underneath the whole comparison: the model we have chosen leaves the two-factor interactions out, and if any of them is not in fact zero, leaving it out pushes its effect onto the coefficients we do estimate. How much, and onto which ones, is read from the alias matrix.

Split the model terms into the ones we keep and the ones we drop. For this four-factor running example, let \(\mathbf{X}_1\) hold the nine fitted columns (the intercept, the four linear terms, and the four pure quadratics) and \(\mathbf{X}_2\) hold the six two-factor interaction columns \(x_i x_j\) left out. If the true response contains those interactions with coefficients \(\boldsymbol{\beta}_2\), the least-squares estimates of the fitted coefficients are biased by a fixed, design-dependent amount:

\[E[\mathbf{b}_1] = \boldsymbol{\beta}_1 + \mathbf{A}\,\boldsymbol{\beta}_2, \qquad \mathbf{A} = (\mathbf{X}_1^T\mathbf{X}_1)^{-1}\mathbf{X}_1^T\mathbf{X}_2\]

Each entry of the alias matrix \(\mathbf{A}\) is the amount by which one omitted interaction shifts one fitted coefficient: an entry of zero leaves that coefficient untouched, an entry of one adds a full unit of the interaction to it. This is the same kind of entanglement the VIF describes, turned outward. The VIF measures overlap among the terms in the model; the alias matrix measures overlap between those terms and the terms the model omits.

The two screening-oriented designs are built to control the rows of \(\mathbf{A}\) that matter most. A definitive screening design and an OMARS design keep every main effect orthogonal to every second-order term, so the main-effect rows of \(\mathbf{A}\) are exactly zero: interactions that are present but omitted do not bias the estimated main effects. That is what minimally aliased names. The price is carried by the quadratics, whose rows are not zero. The omnibus comparison reports the largest absolute entry of \(\mathbf{A}\) for each design, so this bias sits in the same table as the variance it trades against.

5.17.8. Statistical power

Everything so far describes how precisely a design estimates. Power asks the question the experimenter actually cares about: if an effect is really there, how likely are we to detect it? We are testing \(H_0\!: \beta_j = 0\) against the alternative that \(\beta_j\) equals some effect size \(\delta\) we consider practically important. Under that alternative the usual \(t\)- or \(F\)-test statistic is no longer central; it follows a non-central distribution whose non-centrality parameter is

\[\lambda = \frac{\delta^2}{\text{Var}(b_j)} = \frac{\delta^2}{\sigma^2\, c_{jj}}, \qquad c_{jj} = \left[(\mathbf{X}^T\mathbf{X})^{-1}\right]_{jj}\]

The power is the probability that the statistic clears its critical value under this non-central distribution, \(\text{power} = P\!\left(F_{1,\nu} > F_{\text{crit}} \mid \lambda\right)\), where \(\nu\) is the residual degrees of freedom. Anything that improves precision (a larger \(1/c_{jj}\), i.e. more information about that term) raises \(\lambda\) and therefore the power, for a fixed effect size and significance level.

Two consequences deserve to be stated plainly. First, power requires residual degrees of freedom: the non-central distribution needs an estimate of \(\sigma^2\), and a saturated design with \(\nu = 0\) supplies none. This is exactly why the definitive screening design’s power entries in the comparison table below are marked “n/a”: with nine runs and nine terms in the main-effects-and-quadratic model it has nothing left over to estimate the noise, so no test can be run at all. The four extra runs of the thirteen-run OMARS design buy \(\nu = 4\), and with them the ability to test.

Second, power is always quoted for a stated effect size and \(\alpha\). The table’s values assume an effect of one noise standard deviation (\(\delta = \sigma\)) at \(\alpha = 0.05\). Here \(\delta\) is the coefficient itself, so for a linear main effect \(\delta = \sigma\) means the response shifts by \(2\sigma\) across a factor’s \(-1\) to \(+1\) range, a signal-to-noise ratio of 2. That is the default several design packages (for example Stat-Ease Design-Expert) use for power calculations, so the entries here are directly comparable to them; a larger assumed effect would raise every entry in the power rows. Read that way, the thirteen-run OMARS design has roughly a \(0.46\) chance of flagging a true one-sigma main effect as significant (its main effects are not all estimated with equal precision, so the exact figure varies from factor to factor, and the table reports the weakest), and about \(0.25\) for a quadratic of the same size. The gap is expected: quadratic effects are estimated with larger variance (we saw this in the worked example, where \(\text{Var}(b_2)\) was the largest of the three), so they are intrinsically harder to detect. A screening study that must catch curvature will need either more runs or a larger assumed effect size.

5.17.9. Putting the metrics side by side

Collecting every measure we have defined into one place, here are the two designs from the FDS plot, a nine-run definitive screening design and the thirteen-run OMARS design, evaluated on the four-factor main-effects-and-quadratic model.

Comparing the nine-run DSD with the thirteen-run OMARS design.

metric (arrow shows the preferred direction)	DSD [n=9], 4 factors	OMARS [n=13], 4 factors
\(\uparrow\) D-efficiency	42.8 %	39.0 %
\(\downarrow\ A\), summed coefficient variance	3.67	2.52
\(\downarrow\ I\), average SPV	6.59	6.19
\(\downarrow\ G\), maximum SPV	9.00	12.50
\(\downarrow\) maximum \(|r|\)	0.707	0.570
\(\downarrow\) mean \(|r|\)	0.322	0.307
\(\downarrow\) maximum VIF	1.00	1.18
\(\downarrow\) mean VIF	1.00	1.08
\(\uparrow\) residual degrees of freedom	0	4
\(\uparrow\) power, main effect at \(\delta = \sigma\)	n/a	0.46
\(\uparrow\) power, quadratic at \(\delta = \sigma\)	n/a	0.25
\(\uparrow\) two-factor interactions estimable	0	2

No design wins every row, which is the entire point. The thirteen-run OMARS design is better on average coefficient variance, average prediction, both correlation summaries, and is the only one of the two that can estimate interactions or test anything at all; the nine-run DSD holds a higher D-efficiency and a lower worst-case prediction variance \(G\). Read the D-efficiency and the per-run figures with care, though: as the dilemma table showed, these quantities shift with the number of runs, so a head-to-head on \(D\) across a nine-run DSD and a thirteen-run OMARS design is not a like-for-like comparison. The reading that carries weight leans on the quantities that carry real meaning here, separability (\(|r|\), VIF), prediction (\(I\), \(G\)), and the ability to test at all (residual degrees of freedom), and lets the purpose of the study break the ties.

Read this for the method, not for a verdict on the design types. The rows rank these two particular designs on one model; they do not say that an OMARS design is in general better or worse than a definitive screening design. The thirteen-run design here is one member of a large family: for a given factor count the OMARS catalogue holds many designs of different run sizes and aliasing trade-offs, so a different member would sit differently on every row. What transfers to the next study is the procedure, reading precision, separability, and power off the information matrix, not the ranking of any one design.

That is the running comparison in full. The next page widens it from two designs to a shortlist of six, comparing the standard design families for five factors on the same model, and closes with a checklist for choosing among them.