A balanced incomplete block design needs \(\lambda = r(k-1)/(v-1)\) to be a whole number, and even then may need far more blocks than an experimenter can afford. Fisher's inequality \(b \ge v\) alone rules out many useful sizes. The relaxation is to allow two values of \(\lambda\): pairs of treatments are classified as first or second associates, and a pair appears together \(\lambda_1\) or \(\lambda_2\) times according to its class.
The association scheme must itself be regular, in three senses:
The third condition is the substantial one: it is what makes the design's normal equations solvable in closed form.
The second relation counts the other treatments; the third counts, for a fixed treatment, the \(r(k-1)\) units it shares a block with, sorted by associate class — the same argument as relation (ii) for a BIBD, now split in two. A BIBD is the special case \(\lambda_1 = \lambda_2\).
The scheme. Take the six unordered pairs from \(\{1,2,3,4\}\) as the six treatments:
\[ 12,\quad 13,\quad 14,\quad 23,\quad 24,\quad 34. \]Two treatments are first associates if they share a symbol and second associates if they do not. This is the triangular association scheme.
The design. Form one block for each symbol \(x\), containing the three pairs not involving \(x\):
| Block | Symbol omitted | Treatments |
|---|---|---|
| 1 | 1 | 23, 24, 34 |
| 2 | 2 | 13, 14, 34 |
| 3 | 3 | 12, 14, 24 |
| 4 | 4 | 12, 13, 23 |
Step 1 — read off the parameters. \(v = 6\), \(b = 4\), \(k = 3\); each pair avoids two of the four symbols, so \(r = 2\).
Step 2 — the first relation. \(vr = 6 \times 2 = 12\) and \(bk = 4 \times 3 = 12\). \(\checkmark\)
Step 3 — count the associates of \(12\). Sharing a symbol: \(13, 14, 23, 24\), so \(n_1 = 4\). Not sharing one: \(34\) alone, so \(n_2 = 1\). Then \(n_1 + n_2 = 5 = v - 1\). \(\checkmark\)
Step 4 — the two lambdas. \(12\) and \(13\) appear together only in block 4, so \(\lambda_1 = 1\). \(12\) and \(34\) never appear together — block 1 holds \(34\) but not \(12\), block 3 holds \(12\) but not \(34\) — so \(\lambda_2 = 0\).
Step 5 — the third relation.
\[ n_1\lambda_1 + n_2\lambda_2 = 4(1) + 1(0) = 4 = r(k-1) = 2 \times 2. \checkmark \]Step 6 — the \(p\) matrices. Take the first-associate pair \((12, 13)\) and classify the remaining four treatments against each:
\[ \mathbf{P}^{1} = \begin{pmatrix} p^{1}_{11} & p^{1}_{12}\\ p^{1}_{21} & p^{1}_{22}\end{pmatrix} = \begin{pmatrix}2 & 1\\ 1 & 0\end{pmatrix}, \qquad \mathbf{P}^{2} = \begin{pmatrix}4 & 0\\ 0 & 0\end{pmatrix}. \]Check the row sums against \(\sum_w p^{m}_{uw} = n_u - \delta_{um}\): for \(\mathbf{P}^{1}\), row 1 gives \(2 + 1 = 3 = n_1 - 1\) and row 2 gives \(1 + 0 = 1 = n_2\); for \(\mathbf{P}^{2}\), \(4 + 0 = 4 = n_1\) and \(0 + 0 = 0 = n_2 - 1\). \(\checkmark\)
Interpretation. Would a BIBD have done instead? With \(v = 6\), \(k = 3\), \(r = 2\) it would need \(\lambda = r(k-1)/(v-1) = 4/5\), which is not a whole number: no BIBD exists with these parameters. The partially balanced design does, at the cost that not every pair of treatments is compared with the same precision — first associates are compared more precisely than second.
The adjusted totals are formed exactly as for a BIBD, \(Q_i = T_i - \frac{1}{k}\sum_{j \ni i}B_j\), and the normal equations are \(\mathbf{C}\hat{\boldsymbol\tau} = \mathbf{Q}\) with
\[ \mathbf{C} = \frac{r(k-1)}{k}\,\mathbf{I} - \frac{\lambda_1}{k}\,\mathbf{B}_1 - \frac{\lambda_2}{k}\,\mathbf{B}_2, \]where \(\mathbf{B}_m\) is the association matrix with \((i,i')\) entry 1 when \(i\) and \(i'\) are \(m\)th associates and 0 otherwise. Written entrywise this is simply
\[ c_{ii} = \frac{r(k-1)}{k}, \qquad c_{ii'} = -\frac{\lambda_m}{k} \;\text{ when } i, i' \text{ are } m\text{th associates}. \]What the regularity of the scheme buys is that \(\mathbf{C}\) has only three distinct entries, so its generalised inverse also has only three, and the solution can be written in closed form:
\[ \hat\tau_i = \frac{k}{r(k-1)}\left(Q_i + c_1 S_1(Q_i) + c_2 S_2(Q_i)\right), \]where \(S_m(Q_i)\) sums \(Q\) over the \(m\)th associates of \(i\) and the constants \(c_1, c_2\) are fixed by the \(p^{m}_{uw}\). Two treatments are then compared with one of two standard errors, according to whether they are first or second associates — which is the practical meaning of “partially” balanced, and the reason a BIBD is preferred when one exists.
A simple lattice takes \(v = s^{2}\) treatments, writes them in an \(s \times s\) array, and uses two replicates: in the first, the blocks are the rows; in the second, the columns. It is the \(t = 0\) case of the construction at the end of Unit 3 — no Latin squares used at all.
At \(s = 3\):
\[ \begin{array}{ccc} 1 & 2 & 3\\ 4 & 5 & 6\\ 7 & 8 & 9 \end{array} \]| Replicate | Blocks |
|---|---|
| I (rows) | \(\{1,2,3\}\), \(\{4,5,6\}\), \(\{7,8,9\}\) |
| II (columns) | \(\{1,4,7\}\), \(\{2,5,8\}\), \(\{3,6,9\}\) |
The parameters. \(v = 9\), \(b = 6\), \(r = 2\), \(k = 3\), and \(vr = 18 = bk\). \(\checkmark\) Treatment 1 shares a block with \(2, 3\) (its row) and with \(4, 7\) (its column), so \(n_1 = 4\) with \(\lambda_1 = 1\); the remaining four treatments \(5, 6, 8, 9\) never share a block with it, so \(n_2 = 4\) with \(\lambda_2 = 0\). Then
\[ n_1\lambda_1 + n_2\lambda_2 = 4 = r(k-1) = 2 \times 2. \checkmark \]So a simple lattice is a PBIBD(2) with a square association scheme. A BIBD with \(v = 9\), \(k = 3\), \(r = 2\) would need \(\lambda = 2 \times 2 / 8 = 0.5\), which is impossible — the lattice exists precisely where the balanced design cannot.
Resolvability is the practical attraction: each replicate is a complete set of the treatments, so the experiment can be stopped after any whole number of replicates and still be analysable. A triple lattice adds a third replicate built from a Latin square, and a balanced lattice uses all \(s+1\) systems and is then a genuine BIBD with \(\lambda = 1\).
A Latin square controls two nuisance factors at once, but needs as many rows and columns as treatments. When only \(k < v\) rows are available — three positions on a machine and four treatments, three days and four operators — the square cannot be completed. A Youden square is a \(k \times b\) arrangement in which
It is an incomplete Latin square, and the two-way control is real: rows are eliminated exactly as in a Latin square, columns as the blocks of a BIBD.
At \(v = b = 4\), \(k = r = 3\), \(\lambda = 2\):
| Col 1 | Col 2 | Col 3 | Col 4 | |
|---|---|---|---|---|
| Row 1 | A | B | C | D |
| Row 2 | B | C | D | A |
| Row 3 | C | D | A | B |
The four columns are \(\{A,B,C\}, \{B,C,D\}, \{C,D,A\}, \{D,A,B\}\), a BIBD with
\[ vr = 4 \times 3 = 12 = bk = 4 \times 3, \qquad \lambda(v-1) = 2 \times 3 = 6 = r(k-1) = 3 \times 2. \checkmark \]Every treatment appears once in each of the three rows, so the row totals are equal and the rows are orthogonal to the treatments.
The analysis is the BIBD intra-block analysis of Unit 3 with one extra line for rows, taken out before anything else because it is orthogonal to the rest:
| Source | df, general | df here |
|---|---|---|
| Rows | \(k-1\) | 2 |
| Columns (unadjusted) | \(b-1\) | 3 |
| Treatments (adjusted) | \(v-1\) | 3 |
| Error | by subtraction | 3 |
| Total | \(bk-1\) | 11 |
A Youden square exists whenever a symmetric BIBD does, because deleting any one row of a symmetric BIBD's incidence structure and arranging the rest suitably always produces one.
Until now the factor levels have been labels: three varieties, four machines. In response surface methodology they are quantities — a temperature, a concentration, a time — and the object is not to compare the levels tried but to find the setting that optimises the response. The response is treated as a smooth function
\[ y = f(x_1, \ldots, x_k) + e, \]approximated locally by a polynomial. The coded variables are always centred and scaled so that the design points are at \(0, \pm 1\) and perhaps \(\pm\alpha\), which makes the algebra and the interpretation independent of the original units.
The strategy is two-phase. Far from the optimum a plane fits well, and the experiment should move; near the optimum the plane is flat and curvature dominates, and the experiment should map. Those are two different designs.
Fit
\[ y = \beta_0 + \sum_{i=1}^{k}\beta_i x_i + e \]from a \(2^{k}\) factorial (or a resolution III fraction of one) with \(n_0\) runs added at the centre. The factorial points give the \(\beta_i\) with the usual contrast formula; the centre points do two jobs that the factorial points cannot do at all:
on one degree of freedom, against the pure error mean square.
The path of steepest ascent leaves the centre in the direction of the gradient \((\beta_1, \ldots, \beta_k)\), so steps along it are proportional to the fitted coefficients. Choose a convenient step in the variable with the largest coefficient and scale the rest; run single trials along the path until the response stops improving; then start a new first-order experiment there. The path is computed on the coded scale and then translated back, which is why the coding matters.
Given. A \(2^{2}\) factorial with five centre runs, nine observations:
| \(x_1\) | \(x_2\) | \(y\) |
|---|---|---|
| −1 | −1 | 39.3 |
| +1 | −1 | 40.9 |
| −1 | +1 | 40.0 |
| +1 | +1 | 41.5 |
| 0 | 0 | 40.3, 40.5, 40.7, 40.2, 40.6 |
Step 1 — the coefficients. Each is a contrast of the four factorial points divided by 4:
\[ b_1 = \frac{(40.9 + 41.5) - (39.3 + 40.0)}{4} = \frac{82.4 - 79.3}{4} = \frac{3.1}{4} = 0.775000, \] \[ b_2 = \frac{(40.0 + 41.5) - (39.3 + 40.9)}{4} = \frac{81.5 - 80.2}{4} = \frac{1.3}{4} = 0.325000, \] \[ b_{12} = \frac{(39.3 + 41.5) - (40.9 + 40.0)}{4} = \frac{-0.1}{4} = -0.025000. \]The interaction is negligible, which supports the plane.
Step 2 — the intercept. The mean of all nine observations is
\[ b_0 = \frac{364.0}{9} = 40.444444. \]Step 3 — the curvature check.
\[ \bar y_F = \frac{161.7}{4} = 40.425000, \qquad \bar y_C = \frac{202.3}{5} = 40.460000, \qquad \bar y_F - \bar y_C = -0.035000, \] \[ SS_{\text{curvature}} = \frac{4 \times 5 \times (-0.035)^{2}}{9} = \frac{20 \times 0.001225}{9} = 0.002722. \]Step 4 — the pure error. The five centre runs have mean \(40.46\), and
\[ SS_{\text{pure error}} = 0.172000 \text{ on } 4 \text{ df}, \qquad MS_{\text{pure error}} = 0.043000. \]So \(F = 0.002722/0.043 = 0.063\) — no evidence of curvature whatever. The plane is adequate and the experiment should move.
Step 5 — the path of steepest ascent. The gradient is \((0.775, 0.325)\). Taking a step of \(1\) in \(x_1\), the matching step in \(x_2\) is
\[ \frac{b_2}{b_1} = \frac{0.325}{0.775} = 0.419355. \]| Step | \(x_1\) | \(x_2\) |
|---|---|---|
| 1 | 1.0000 | 0.4194 |
| 2 | 2.0000 | 0.8387 |
| 3 | 3.0000 | 1.2581 |
| 4 | 4.0000 | 1.6774 |
| 5 | 5.0000 | 2.0968 |
Interpretation. A single observation is taken at each point of the path until the response falls. The path is a search, not an experiment: its runs are unreplicated and no test is performed on them. When the response stops improving, a fresh \(2^{2}\) with centre points is run at the best point found, and the curvature test is repeated. When it becomes significant, the plane has done its work and the second-order design of section 5 is called for.
with \(1 + 2k + \binom{k}{2}\) parameters — six at \(k = 2\), ten at \(k = 3\). A \(2^{k}\) factorial cannot fit it: with only two levels per factor there is no way to estimate \(\beta_{ii}\), since \(x_i^{2} = 1\) at every point. At least three levels of each factor are needed.
Rotatability. Before the experiment, the location of the optimum is unknown, so there is no direction in which precision should be preferred. A design is rotatable if the variance of the predicted response, \(\operatorname{Var}(\hat y(\mathbf{x}))\), depends on \(\mathbf{x}\) only through the distance \(\rho = \sqrt{\sum_i x_i^{2}}\) from the centre. Rotating the design about the centre then changes nothing, which is exactly the neutrality wanted.
The condition, in terms of the design moments, is
\[ \sum_u x_{iu} = \sum_u x_{iu}x_{ju} = \sum_u x_{iu}^{3} = \sum_u x_{iu}^{2}x_{ju} = 0, \qquad \sum_u x_{iu}^{4} = 3\sum_u x_{iu}^{2}x_{ju}^{2} \;\; (i \ne j). \]The first group is symmetry, and any balanced design has it. The second is the real condition, and it is a constraint on how far out the extra points are placed.
A central composite design adds two things to a \(2^{k}\) factorial:
Total runs \(N = F + 2k + n_0\) where \(F = 2^{k}\) (or a fraction of it). The design is built in two stages, so a first-order experiment can be augmented into a second-order one without discarding anything — which is its main practical attraction.
Choosing \(\alpha\) for rotatability. With \(F\) factorial points at \(\pm 1\), \(2k\) axial points at \(\pm\alpha\) and \(n_0\) centre points,
\[ \sum_u x_{iu}^{2} = F + 2\alpha^{2}, \qquad \sum_u x_{iu}^{4} = F + 2\alpha^{4}, \qquad \sum_u x_{iu}^{2}x_{ju}^{2} = F, \]because the axial points have one non-zero coordinate each and so contribute nothing to the cross moment. The rotatability condition \(\sum x_i^{4} = 3\sum x_i^{2}x_j^{2}\) becomes
\[ F + 2\alpha^{4} = 3F \;\Longrightarrow\; \alpha^{4} = F \;\Longrightarrow\; \boxed{\;\alpha = F^{1/4}\;} \]— a remarkably simple answer, and one that depends only on the number of factorial points, not on \(n_0\).
| \(k\) | \(F = 2^{k}\) | \(\alpha = F^{1/4}\) | \(N\) with \(n_0\) centre runs |
|---|---|---|---|
| 2 | 4 | 1.414214 | \(8 + n_0\) |
| 3 | 8 | 1.681793 | \(14 + n_0\) |
| 4 | 16 | 2.000000 | \(24 + n_0\) |
Note what \(n_0\) does not do: it does not affect rotatability at all. What it affects is the shape of the variance function, and choosing it so that the variance is nearly constant over the region of interest gives a uniform precision design — \(n_0 = 5\) at \(k = 2\), \(n_0 = 6\) at \(k = 3\).
Given. The \(k = 2\) central composite design of Fig 4.1: \(F = 4\), \(\alpha = \sqrt 2\), \(n_0 = 5\), so \(N = 13\).
Step 1 — the moments.
\[ \sum_u x_{1u}^{2} = 4 + 2(\sqrt 2)^{2} = 4 + 4 = 8, \qquad \sum_u x_{1u}^{4} = 4 + 2(\sqrt 2)^{4} = 4 + 8 = 12, \] \[ \sum_u x_{1u}^{2}x_{2u}^{2} = 4 \quad\text{(only the factorial points contribute)}. \]Step 2 — test rotatability. \(3 \times 4 = 12 = \sum x_1^{4}\). \(\checkmark\) The design is rotatable.
Step 3 — the variance function. Writing \(\mathbf{x}_m' = (1, x_1, x_2, x_1^{2}, x_2^{2}, x_1x_2)\),
\[ \frac{\operatorname{Var}\left(\hat y(\mathbf{x})\right)}{\sigma^{2}} = \mathbf{x}_m'\left(\mathbf{X}'\mathbf{X}\right)^{-1}\mathbf{x}_m. \]Evaluating it at six different directions at each of several distances gives:
| \(\rho\) | along \(x_1\) | at 30° | at 45° | at 60° | along \(x_2\) | spread |
|---|---|---|---|---|---|---|
| 0.000000 | 0.200000 | 0.200000 | 0.200000 | 0.200000 | 0.200000 | 0 |
| 0.500000 | 0.190234 | 0.190234 | 0.190234 | 0.190234 | 0.190234 | 0 |
| 1.000000 | 0.268750 | 0.268750 | 0.268750 | 0.268750 | 0.268750 | 0 |
| 1.414214 | 0.625000 | 0.625000 | 0.625000 | 0.625000 | 0.625000 | 0 |
| 2.000000 | 2.200000 | 2.200000 | 2.200000 | 2.200000 | 2.200000 | 0 |
The six figures in each row are identical to twelve decimal places, which is what rotatability means in practice.
Step 4 — read the shape. The variance falls from \(0.200000\) at the centre to \(0.190234\) at \(\rho = 0.5\), then rises to \(0.268750\) at the edge of the factorial and \(0.625000\) at the axial points. Prediction is most precise just inside the factorial region, not at the very centre, and it deteriorates sharply beyond the design — at \(\rho = 2\) the variance is eleven times its value at \(\rho = 0.5\).
Interpretation. The flat interior is the uniform-precision property that \(n_0 = 5\) was chosen for. The sharp rise outside is the reason a response surface must never be used to extrapolate: the fitted second-order model will happily report an optimum outside the design region, and its standard error there makes the report worthless.
Once the second-order model is fitted, write it in matrix form
\[ \hat y = b_0 + \mathbf{x}'\mathbf{b} + \mathbf{x}'\mathbf{B}\mathbf{x}, \qquad \mathbf{B} = \begin{pmatrix} b_{11} & b_{12}/2\\ b_{12}/2 & b_{22}\end{pmatrix}, \]the off-diagonal halving being because \(b_{12}x_1x_2\) is split between two entries. Setting the gradient \(\mathbf{b} + 2\mathbf{Bx}\) to zero gives the stationary point
\[ \mathbf{x}_s = -\tfrac12\,\mathbf{B}^{-1}\mathbf{b}, \qquad \hat y_s = b_0 + \tfrac12\,\mathbf{x}_s'\mathbf{b}. \]Its nature is settled by the eigenvalues of \(\mathbf{B}\), by the canonical form of Linear Algebra, Unit 3: all negative gives a maximum, all positive a minimum, mixed signs a saddle, and one near zero a ridge along which the response barely changes. A ridge is often the most useful outcome of all, because it means the process can be run anywhere along a line without loss — and the cheapest point on that line can then be chosen.