Blog
EDM 1 - 2 - 102 :: Part - 3 :: Design of Experiments and Research Method Concepts
Filed under MathematicsScienceEconomicsBooksDesign
J. Application of Algebras
There is a rich connection between experimental design and algebraic structures, particularly combinatorial matrix theory and group theory. In the theory of experimental designs, one useful concept is the association algebra (or design algebra) which arises from treating treatments or experimental units as vertices in a structure and defining relationships among them.
For block designs, one can define an association scheme based on how often pairs of treatments occur together in blocks. Specifically, in a block design with treatments, define a relation on the set of treatments such that if treatment and occur together in exactly blocks. For example, in a BIBD, two treatments are either in blocks together or not at all, so there are typically two such relations (ignoring for not co-occurring, and for co-occurring). These relations can be represented by (0,1)-matrices , where if and 0 otherwise.
The span of these matrices (including and often the all-ones matrix for completeness) can form a commutative algebra known as the Bose–Mesner algebra of an association scheme. The matrices have properties such as:
and
with constants that count how many treatment pairs in relation and a given treatment in relation imply a relation . For a BIBD, this algebra is relatively simple because there are just two relations (either two treatments co-occur or not), and it leads to the properties of the information matrix being proportional to structure mentioned.
Why is this algebra useful? One reason is it provides a systematic way to obtain the eigenvalues of the information matrix or C-matrix of the design. The matrix for a balanced design can often be expressed as a linear combination of these matrices. Then the eigen-decomposition of (which determines variances of estimates) can be obtained from the known eigenvectors of the association matrices (which typically have a few eigenvalues with known multiplicities from the association scheme theory). For example, in a BIBD, the treatment information matrix (adjusted for one constraint) has two eigenvalues: one of multiplicity 1 (for the mean direction) and one of multiplicity (for the contrasts), resulting directly from the two eigenvalues of the matrices in the association scheme.
Another algebraic approach is through the use of Hadamard matrices and orthogonal arrays for constructing designs. An orthogonal array can be seen as a set of vectors of length over an alphabet of sise with certain orthonormality properties. The rows of an orthogonal array generate an algebra (under element-wise multiplication mod if we treat levels as elements of a group) that helps in analyzing confounding. For two-level designs, the design generators can be considered as forming a vector space over (the field with 2 elements). The aliasing is determined by the structure of this vector space. Algebraic coding theory methods apply to analyze these designs: the design is like a linear code where the minimum distance is related to resolution.
In summary, algebraic methods provide powerful tools for:
-
Characterising designs: using parameters and association schemes, one can derive combinatorial constraints (like Fisher’s inequality, BIBD existence conditions, etc.) and optimality criteria in elegant ways.
-
Analyzing designs: obtaining eigenvalues of matrices, which then yield formulas for variances of estimates and efficiency factors.
-
Constructing designs: e.g., using finite field arithmetic to create orthogonal arrays, or using group theory to construct PBIBD (partially balanced incomplete block designs) where treatments are arranged in association classes (like PBIBD with two associate classes corresponds to a -class association scheme).
One concrete example: Partially Balanced Incomplete Block Designs (PBIBD). These are incomplete block designs where not all pairs of treatments occur together equally often, but there are a few classes of co-occurrences. For instance, in a PBIBD with 2 associate classes, any two treatments are either first associates (occur together times) or second associates (occur together times), and these relationships form an association scheme with two classes. The matrices and (plus identity) form an algebra. The analysis of a PBIBD yields two different eigenvalues for treatment comparisons, meaning the variances of estimates might not all be equal but can be described succinctly by two values. The ANOVA of a PBIBD can be performed by decomposing treatment sums of squares according to these eigencomponents (sometimes called contrasts corresponding to eigenvectors of and ). PBIBDs are less simple than BIBDs but allow more flexibility in parameters; algebraic techniques are crucial in their construction (often via finite geometries or group-divisible designs) and analysis.
Algebraic methods in design theory connect to topics like group theory (e.g., using cyclic groups to generate cyclic designs), field theory (finite fields to generate factorial designs and orthogonal arrays), and linear algebra (spectral decomposition of design matrices). They provide a unifying framework to understand seemingly disparate designs.
K. Design for Two-Way Elimination of Heterogeneity
When nuisance variability exists in two different directions or criteria, a two-way elimination of heterogeneity can be performed by using a row-column design. The classic example of this is a Latin Square Design. A Latin square controls variation in two orthogonal directions (often called “rows” and “columns”) in addition to the treatments.
A Latin square of order is an arrangement of treatments in a grid such that each treatment appears exactly once in each row and exactly once in each column. Thus, treatments, rows, and columns. The design requires experimental units arranged in a square (or in rows and columns). One way to visualise this: think of a field that is divided into rows and columns (thus forming plots). We have treatment levels to apply, and we want to ensure each treatment occurs once per row and once per column.
Statistically, the model for a Latin square design is:
where indexes the treatment, the row, and the column. Here is the treatment effect, the row effect, and the column effect. The Latin square property ensures treatments and rows and columns are orthogonal in the design matrix: each treatment appears once in each row and once in each column. Typically we impose constraints for identifiability.
The analysis of a Latin square partitions total variance into:
- Rows ( df),
- Columns ( df),
- Treatments ( df),
- Error ( df).
Why error df ? Because total df is , and we remove df for row, column, and treatment. Simplify: . For this is valid (the smallest Latin square is , but that yields df for error, not useful for analysis, so must be at least 3 for an actual experiment).
Latin squares are particularly useful when there are two nuisance factors (e.g., in an agricultural trial, varieties (treatments) might be tested where the field has a fertility gradient in two perpendicular directions; rows and columns account for those gradients). Another application is in taste testing or controlled experiments where rows could be time slots and columns could be locations, etc., and one wants each treatment to appear exactly once in each time and each location.
Graeco-Latin Squares: A further extension is a Graeco-Latin square, which involves two different sets of treatments (say Latin letters and Greek letters). A Graeco-Latin square of order is a pair of Latin squares of order , one using Latin symbols and one using Greek symbols, arranged such that when superimposed, each Latin letter occurs exactly once with each Greek letter. Graeco-Latin squares allow studying two treatment factors while controlling for two blocking factors (rows and columns). They require treatments in set A and in set B (two treatment factors each of sise ), and use rows and columns as blocks. It’s like adding a second set of treatments on top of a Latin square structure. These exist for many orders (notably all except have Graeco-Latin squares).
ANOVA of Latin Square: The F-test for treatments uses MS$_{\text{trt}} / $ MS with degrees of freedom. Row and column effects are usually not of intrinsic interest but are accounted for; one can also check if they are significant to verify block effectiveness.
Latin square designs assume no interaction between row and treatment, column and treatment, or row and column (which is why each appears once in each combination). If there are potential interactions (e.g., if a treatment reacts differently to different rows), the Latin square does not allow estimating that interaction separately; such interaction would inflate error variance or bias treatment comparisons if systematic. Therefore, Latin squares are primarily for situations where the blocking factors do not interact with treatments.
Extensions one can extend the concept to higher dimensions:
-
Youden Squares: A Youden square is a rectangular design (not necessarily square) that is essentially a Latin square with one dimension incomplete. It allows one missing blocking level in one direction.
-
Latin Rectangle: If the number of treatments is less than , one can arrange treatments in a rectangle such that each treatment appears at most once per row and at most once per column. However, if , a full Latin square cannot be formed, but partial arrangements exist.
-
Graeco-Latin Cube: A hypercube for more blocking criteria can be conceptualised but becomes impractical beyond two blocking factors.
Latin square designs are widely used in agricultural and industrial experiments when two orthogonal blocking factors are present. They achieve a very efficient design: only replications for treatments yet controlling two nuisance variables.
L. Balanced Array and Orthogonal Array
The concept of balanced arrays and orthogonal arrays is a generalisation of the idea of fractional factorial designs and other orthogonal designs to possibly more levels and factors.
A balanced array refers to an arrangement of experimental runs where certain combinations of factor levels occur equally often. An orthogonal array is a specific kind of balanced array with a strict combinatorial definition:
An is an array filled with symbols (levels) such that in any selection of columns, each of the possible level combinations appears in exactly runs. Here is the number of runs (rows), is the number of factors (columns), is the number of levels each factor can take (assumed equal for all factors in a standard orthogonal array), and is the strength (the array is orthogonal up to any columns).
Orthogonal arrays are a general design tool. Many designs can be seen as orthogonal arrays:
-
A full factorial of levels in factors is denoted by .
-
A fractional factorial that is orthogonal for all main effects and some interactions is an , where equals the maximum interaction order for which orthogonality holds.
-
A Latin square of order is an : runs, factors (row, column, treatment), each with levels, and strength ; meaning any two columns give all pairs exactly once.
-
A Graeco-Latin square of order is an with factors: row, column, Latin treatment, and Greek treatment. Any pair of factors yields all level combinations exactly once, satisfying the Graeco-Latin property.
-
A two-factor fractional factorial design of resolution III is an : strength 2 means that any pair of factors is balanced (i.e., main effects are not aliased with each other).
-
An array of strength 3, i.e., , corresponds to a resolution IV fractional factorial design. Such a design ensures that any three columns; such as one main effect and one two-factor interaction; produce all combinations equally often.
Orthogonal arrays ensure that projections onto a small number of factors give full coverage of level combinations, which implies that main effects (with ) and some interactions (up to factor interactions) can be estimated without confounding.
Taguchi popularised the use of orthogonal arrays in robust parameter design. For example, Taguchi’s design is an : 8 runs, 7 factors, each at 2 levels, strength 2 (it’s a specific half-fraction of a that is resolution III).
Orthogonal arrays can often be constructed from mathematical objects. For instance: A Hadamard matrix of order (when is multiple of 4 beyond 4) gives an , since a Hadamard matrix of order yields runs for factors (excluding one factor which would be the product of all others, akin to one generator) at 2 levels with balance.
Finite field constructions: For example, an orthogonal array can be constructed using a Latin square of order . Another example: using a -ary linear code of length and dimension (over ) yields an if the code has no short blocks (this is related to error-correcting codes).
The transpose of the incidence matrix of a symmetric BIBD yields an if one treats blocks as runs and each treatment as a factor with two levels (0/1 indicating absence/presence in block). However, usually orthogonal arrays are thought of with equal levels for each factor. A BIBD translates to an orthogonal array with factors having 2 levels (treatment absent or present), but each factor (treatment) appears times in runs and each pair of factors co-occurs times, giving strength 2 with some parameter constraints (regular fractional).
One particularly important family is the Plackett-Burman designs, which are for multiple of 4 (other than 4 itself yields trivial). These are resolution III designs (main effects nonconfounded, but main effects aliased with some two-factor interactions) often used for screening. They are not full factorial fractions but rather non-regular fractions constructed via Hadamard matrices.
Orthogonal arrays and balanced arrays also play a role in response surface methodology (RSM) for selecting sample points that fit higher-order models efficiently (though for RSM, designs like central composite designs or Box-Behnken designs are commonly used, which are not exactly orthogonal arrays but have other properties like rotatability).
In conclusion, balanced and orthogonal arrays provide a unifying framework for understanding and constructing experimental designs that achieve balance and orthogonality in factor level combinations. They facilitate highly economical designs (fractional factorials, Latin/Graeco-Latin squares, etc.) which still allow clear estimation of certain effects. The use of orthogonal arrays is especially prominent in industrial and engineering experimentation, where they are used as standard test plans for robust product design (to study many factors in few runs).
M. Response Surface
So far we have largely discussed designs for determining the presence or absence of factor effects or comparing treatment means. In many experiments, especially in industrial and chemical processes, the goal is to optimise a response that depends on quantitative factors. Response surface methodology (RSM) deals with experiments where the factors are continuous variables (like temperature, pressure, time, concentrations, etc.) and the objective is to model the response as a smooth function of these variables, then find the settings that optimise (maximise or minimise) the response or achieve a target.
We consider quantitative factors . The response is a function of these factors:
Usually is unknown and we approximate it by a polynomial model based on experimental data. For instance, a first-order model (a linear model) or a second-order model (including quadratic terms and possibly interactions) are common approximations.
In RSM, experiments are often done in stages:
-
Factor Screening: If is large, first use fractional factorial or Plackett-Burman designs at two levels to identify the subset of factors that significantly influence .
-
Steepest Ascent/Descent: Using a first-order model, identify a direction in the factor space where the response seems to improve and move sequentially in that direction (adjust factors step by step) until no further improvement is observed. This finds a region near the optimum.
-
Optimisation Experiment: Once near the optimum, run a second-order design to fit a full quadratic model for . Then analytically find the stationary point of that quadratic model and determine if it is a maximum, minimum, or saddle point. This stationary point is an estimate of the optimum conditions.
Common designs in RSM: Central Composite Design (CCD): This is a popular second-order design. For factors, a CCD consists of:
-
A factorial or fractional factorial portion (often or runs at the cube points, usually at coded levels and for each factor).
-
An axial (star) portion: runs where one factor is set to an extreme level (often coded or ) and all other factors at the center (0). These points allow estimation of pure quadratic terms.
-
Center points: several runs at the center to estimate pure error and check for curvature.
The axial distance can be chosen to ensure certain properties: for example, yields rotatability, meaning the variance of predicted depends only on the distance from the design center, not the direction. Another common choice is when augmenting a factorial, though that doesn’t give rotatability unless or design is scaled.
Box-Behnken Design: Another second-order design for factors that does not require extreme factorial corners. It is constructed by combining -level factorial designs for subsets of factors. For example, for , a Box-Behnken has midpoints of edges of the cube and center (totals 13 runs for 3 factors). For general , Box-Behnken requires factors be 3 levels each but uses fewer runs than a full . It doesn’t include corners of the cube, only middle of edges and center, which is safer if extreme conditions are hard to run. Box-Behnken designs are nearly rotatable and have good efficiency for quadratic model estimation.
3-level factorial or fractional factorial: A full factorial provides a full second-order fit (since a quadratic in each factor can be captured by 3 points per factor), but grows quickly. Often a fractional is used, analogous to 2-level fractions but for 3 levels (though ensuring estimability of all needed quadratic terms requires careful selection of fraction).
Key concepts in RSM:
-
Rotatability: A design is rotatable if the variance of the predicted response at a point depends only on the distance from the design center, not on the direction. Rotatable designs allow uniform precision in all directions around the center, useful when the direction of the optimum is unknown. Many CCDs and Box-Behnken designs are or can be made approximately rotatable.
-
Scaled Units and Coding: Factors are usually coded to dimensionless units (like , , ) to make all factors comparable and to use design tables. It also improves numerical stability in regression.
-
Lack of Fit and Pure Error: By including replicate center points (and possibly other replicates), one can test for lack-of-fit of the fitted polynomial model. Lack-of-fit is detected if the residual error is significantly larger than the pure error (from replicated points).
-
Response Optimisation and Canonical Analysis: Once a quadratic model is fit, it can be written in matrix form:
where is a symmetric matrix of quadratic coefficients. Setting the derivative to zero,
yields the stationary point
By eigen-decomposing , one can classify the stationary point: if all eigenvalues of are negative, it’s a local maximum; if all are positive, a local minimum; mixed signs indicate a saddle point. The magnitudes indicate curvature in each principal direction. This analysis is known as canonical analysis of the response surface.
- Sequential Nature: RSM encourages sequential designs: starting with a fractional factorial (resolution III or IV) for screening, followed by a CCD or similar augmentation, rather than executing a large design in one go. This approach is adaptive and cost-effective.
For example, consider an experiment with two factors (temperature and pressure ) that influence yield . We might first run a factorial plus center to see if a linear model suffices. If curvature is detected, we augment with axial points to make it a CCD. Suppose after fitting a quadratic model, we find:
with coded such that 0 is current operating condition, +1 and -1 correspond to high/low settings. We can solve for optimum by setting derivatives zero: , . Solving: from first, . From second, . Solve these simultaneously: from second, . Substitute into first: . That gives . Multiply by 3: . So . Rearr: , , . Then . In coded units that’s beyond the original range (+1), meaning the optimum is outside initial factor range in . If this is so, one might need to extend the region or accept boundary optimum. Also check classification by second derivatives matrix: in this model (half of interaction coefficient in each off-diagonal). Eigenvalues might be negative (likely yes, because diag -2, -4 and not too large off diag, so looks negative definite, implying a local maximum). So it suggests maximum yield at (which is slightly above our high temp) and (just slightly above center pressure).
The actual optimum yield predicted can be computed by plugging these in the model. If the region had to be extended, the experimenter would need to run a new set of experiments covering higher temperature to verify and fine-tune optimum.
In practice, process optimisation via RSM has led to significant improvements in manufacturing, chemistry, etc., by finding better operating conditions.
One must also consider validation after finding an optimum via RSM, it’s common to run confirmatory experiments at the predicted optimum to verify the predicted gain in response.
RSM underlines the philosophy of DOE as an iterative process: use designed experiments to model the system and then exploit the model to improve the system.