Design of Experiments (DOE)

Design of Experiments (DOE) is a systematic, rigorous methodology for planning controlled experiments that enable efficient estimation of factor effects and their interactions. Unlike observational studies, DOE allows researchers to establish causal relationships by systematically varying experimental factors while controlling extraneous sources of variation. The mathematical foundations of DOE draw from linear algebra, combinatorics, and the theory of orthogonal arrays, providing a framework for optimal information extraction with minimal experimental effort.

Fundamental Principles

Factorial Designs

Factorial designs are the cornerstone of DOE, simultaneously investigating multiple factors and their interactions.

Confounding and Resolution

When experimental resources are limited, fractional factorial designs confound (alias) some effects, assuming higher-order interactions are negligible.

Advanced Design Structures

Box-Behnken Designs

Box-Behnken designs are three-level designs that do not contain extreme points (corner points of the cube), making them suitable for second-order response surface modeling.

Central Composite Design (CCD)

Python Implementation

import numpy as np
import pandas as pd
from itertools import product
from scipy import stats
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Generate 2^k factorial design
def factorial_design(k, factors=None):
    """Generate a full 2^k factorial design."""
    if factors is None:
        factors = [f'X{i+1}' for i in range(k)]
    
    # Generate all combinations of -1 and +1
    levels = list(product([-1, 1], repeat=k))
    design = pd.DataFrame(levels, columns=factors)
    
    # Add interaction columns
    for i in range(k):
        for j in range(i+1, k):
            col_name = f'{factors[i]}*{factors[j]}'
            design[col_name] = design[factors[i]] * design[factors[j]]
    
    return design

# Effect estimation
def estimate_effects(design, response):
    """Estimate main effects and interactions."""
    effects = {}
    k = len([col for col in design.columns if '*' not in col])
    
    # Main effects
    for j in range(k):
        factor = design.columns[j]
        effects[factor] = np.mean(response[design[factor] == 1]) - \
                         np.mean(response[design[factor] == -1])
    
    # Two-factor interactions
    for i in range(k):
        for j in range(i+1, k):
            f1, f2 = design.columns[i], design.columns[j]
            interaction = design[f1] * design[f2]
            effects[f'{f1}*{f2}'] = np.mean(response[interaction == 1]) - \
                                   np.mean(response[interaction == -1])
    
    return effects

# Generate 2^3 factorial design
design_2k3 = factorial_design(3, ['Temperature', 'Pressure', 'Catalyst'])
print("2^3 Factorial Design:")
print(design_2k3.head(8))

# Simulated yield data
np.random.seed(42)
n_replicates = 3
design_expanded = pd.DataFrame(np.repeat(design_2k3.values, n_replicates, axis=0),
                               columns=design_2k3.columns)

# True effects
true_effects = {'Temperature': 8.0, 'Pressure': 5.0, 'Catalyst': 3.0,
                'Temperature*Pressure': -2.0, 'Temperature*Catalyst': 1.5,
                'Pressure*Catalyst': 0.8}

# Generate response with effects and noise
yield_data = []
for _, row in design_expanded.iterrows():
    y = 50  # Grand mean
    y += (row['Temperature'] * true_effects['Temperature'] / 2)
    y += (row['Pressure'] * true_effects['Pressure'] / 2)
    y += (row['Catalyst'] * true_effects['Catalyst'] / 2)
    y += (row['Temperature'] * row['Pressure'] * true_effects['Temperature*Pressure'] / 2)
    y += (row['Temperature'] * row['Catalyst'] * true_effects['Temperature*Catalyst'] / 2)
    y += (row['Pressure'] * row['Catalyst'] * true_effects['Pressure*Catalyst'] / 2)
    y += np.random.normal(0, 1.5)
    yield_data.append(y)

design_expanded['Yield'] = yield_data

# Estimate effects
effects = estimate_effects(design_2k3, 
                          [design_expanded['Yield'].iloc[i:i+n_replicates].mean() 
                           for i in range(0, len(design_expanded), n_replicates)])

print("\nEstimated Effects:")
for effect, value in effects.items():
    print(f"  {effect}: {value:.3f}")

# Fractional factorial 2^(3-1)
def fractional_factorial_2k1():
    """Generate 2^(3-1) fractional factorial design."""
    # Generator: C = AB
    design = pd.DataFrame({
        'A': [-1, -1, 1, 1],
        'B': [-1, 1, -1, 1],
        'C': [1, -1, -1, 1]  # C = A*B
    })
    return design

design_frac = fractional_factorial_2k1()
print("\n2^(3-1) Fractional Factorial:")
print(design_frac)
print("Alias structure: A aliased with BC, B aliased with AC, C aliased with AB")

# Box-Behnken design for 3 factors
def box_behnken_3factor():
    """Generate Box-Behnken design for 3 factors."""
    # Edge midpoints of the cube
    runs = [
        [-1, -1, 0], [-1, 1, 0], [1, -1, 0], [1, 1, 0],  # AB face
        [-1, 0, -1], [-1, 0, 1], [1, 0, -1], [1, 0, 1],  # AC face
        [0, -1, -1], [0, -1, 1], [0, 1, -1], [0, 1, 1],  # BC face
        [0, 0, 0], [0, 0, 0], [0, 0, 0]  # Center points
    ]
    return pd.DataFrame(runs, columns=['X1', 'X2', 'X3'])

bb_design = box_behnken_3factor()
print("\nBox-Behnken Design (3 factors):")
print(bb_design)

# Central Composite Design for 2 factors
def central_composite_2factor(alpha=None):
    """Generate CCD for 2 factors."""
    if alpha is None:
        alpha = np.sqrt(2)  # Rotatable
    
    factorial = [[-1, -1], [-1, 1], [1, -1], [1, 1]]
    axial = [[-alpha, 0], [alpha, 0], [0, -alpha], [0, alpha]]
    center = [[0, 0], [0, 0]]
    
    design = pd.DataFrame(factorial + axial + center, 
                          columns=['X1', 'X2'])
    return design

ccd_design = central_composite_2factor()
print("\nCentral Composite Design (2 factors, rotatable):")
print(ccd_design)

# Visualize CCD
fig, ax = plt.subplots(figsize=(8, 8))
factorial_pts = ccd_design.iloc[:4]
ax.scatter(factorial_pts['X1'], factorial_pts['X2'], s=100, 
          c='blue', label='Factorial', zorder=3)
ax.scatter(ccd_design.iloc[4:8]['X1'], ccd_design.iloc[4:8]['X2'], s=100, 
          c='red', label='Axial', zorder=3)
ax.scatter(ccd_design.iloc[8:]['X1'], ccd_design.iloc[8:]['X2'], s=100, 
          c='green', label='Center', zorder=3)
ax.axhline(0, color='gray', linestyle='--', alpha=0.5)
ax.axvline(0, color='gray', linestyle='--', alpha=0.5)
ax.set_xlabel('X1')
ax.set_ylabel('X2')
ax.set_title('Central Composite Design')
ax.legend()
ax.set_aspect('equal')
plt.grid(True, alpha=0.3)
plt.savefig('ccd_design.png', dpi=150)
plt.show()

Key Design Properties