Chapter 1: Fundamentals of DOE and Orthogonal Arrays

Understand the basic concepts of Design of Experiments (DOE) and learn efficient experimental design using orthogonal arrays. Master from one-way and two-way experiments to interpretation of main effect plots and interaction plots through Python practice.

1.1 What is Design of Experiments (DOE)

Definition and Purpose of DOE

Design of Experiments (DOE) is a statistical method to obtain maximum information with minimum number of experiments. It aims to efficiently evaluate the effects of multiple factors (variables) on products or processes and find optimal conditions.

Main purposes:

• Factor screening: Identify important factors among many candidates
• Optimization: Discover conditions that maximize process or product performance
• Robustness evaluation: Improve stability against disturbances
• Interaction detection: Understand mutual influences between factors

Differences from Conventional One-Variable-at-a-Time Experiments

Example: Evaluating 3 factors (temperature, pressure, catalyst amount) with 3 levels each in a chemical reaction

• OFAT: 3 × 3 × 3 = 27 experiments (evaluating each factor individually)
• DOE (Orthogonal array L9): 9 experiments evaluate 3 factors and their interactions

Three Principles of DOE

1. Replication
   • Conduct multiple experiments under the same conditions
   • Enables estimation of experimental error and statistical testing
2. Randomization
   • Randomize experimental order
   • Eliminate effects of time trends and systematic errors
3. Blocking
   • Group by known disturbance factors
   • Separate effects such as batch-to-batch variation

1.2 One-Way Experiment and ANOVA

Code Example 1: One-way ANOVA

Compare the performance of three types of catalysts (A, B, C) in a chemical reaction.

# Requirements:
# - Python 3.9+
# - matplotlib>=3.7.0
# - numpy>=1.24.0, <2.0.0
# - pandas>=2.0.0, <2.2.0
# - scipy>=1.11.0
# - seaborn>=0.12.0

"""
Example: Compare the performance of three types of catalysts (A, B, C

Purpose: Demonstrate data visualization techniques
Target: Intermediate
Execution time: 2-5 seconds
Dependencies: None
"""

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
import seaborn as sns

# Experimental data: Reaction yield (%) with 3 types of catalysts
# 5 experiments conducted for each catalyst
np.random.seed(42)

catalyst_A = [85.2, 84.8, 86.1, 85.5, 84.9]
catalyst_B = [88.3, 89.1, 88.7, 89.5, 88.2]
catalyst_C = [86.5, 87.2, 86.8, 87.0, 86.3]

# Organize into dataframe
data = pd.DataFrame({
    'Catalyst': ['A']*5 + ['B']*5 + ['C']*5,
    'Yield': catalyst_A + catalyst_B + catalyst_C
})

print("=== Experimental Data ===")
print(data.groupby('Catalyst')['Yield'].describe())

# One-way ANOVA
groups = [catalyst_A, catalyst_B, catalyst_C]
f_stat, p_value = stats.f_oneway(*groups)

print(f"\n=== One-way ANOVA ===")
print(f"F-statistic: {f_stat:.4f}")
print(f"p-value: {p_value:.6f}")

if p_value < 0.05:
    print("Conclusion: Significant difference among catalysts at 5% significance level")
else:
    print("Conclusion: No significant difference among catalysts")

# Visualization: Box plot
plt.figure(figsize=(10, 6))
sns.boxplot(x='Catalyst', y='Yield', data=data, palette='Set2')
plt.title('Comparison of Yield by Catalyst Type', fontsize=14, fontweight='bold')
plt.ylabel('Yield (%)', fontsize=12)
plt.xlabel('Catalyst', fontsize=12)
plt.grid(alpha=0.3)
plt.tight_layout()
plt.savefig('one_way_anova.png', dpi=300, bbox_inches='tight')
plt.show()

# Mean and 95% confidence interval for each catalyst
print("\n=== Mean and 95% Confidence Interval for Each Catalyst ===")
for catalyst in ['A', 'B', 'C']:
    subset = data[data['Catalyst'] == catalyst]['Yield']
    mean = subset.mean()
    ci = stats.t.interval(0.95, len(subset)-1, loc=mean, scale=stats.sem(subset))
    print(f"Catalyst {catalyst}: Mean={mean:.2f}%, 95%CI=[{ci[0]:.2f}, {ci[1]:.2f}]")

Example output:

=== Experimental Data ===
           count   mean       std   min    25%    50%    75%    max
Catalyst
A            5.0  85.30  0.487852  84.8  84.90  85.2  85.50  86.1
B            5.0  88.76  0.543139  88.2  88.30  88.7  89.10  89.5
C            5.0  86.76  0.382099  86.3  86.50  86.8  87.00  87.2

=== One-way ANOVA ===
F-statistic: 55.9821
p-value: 0.000002
Conclusion: Significant difference among catalysts at 5% significance level

=== Mean and 95% Confidence Interval for Each Catalyst ===
Catalyst A: Mean=85.30%, 95%CI=[84.70, 85.90]
Catalyst B: Mean=88.76%, 95%CI=[88.09, 89.43]
Catalyst C: Mean=86.76%, 95%CI=[86.29, 87.23]

Interpretation: With a large F-statistic and p-value less than 0.05, catalyst B has the highest yield, followed by C and A, with statistically significant differences.

Code Example 2: Two-way ANOVA without Replication

Evaluate the effects of temperature and pressure on chemical reaction yield.

# Requirements:
# - Python 3.9+
# - matplotlib>=3.7.0
# - numpy>=1.24.0, <2.0.0
# - pandas>=2.0.0, <2.2.0
# - scipy>=1.11.0
# - seaborn>=0.12.0

"""
Example: Evaluate the effects of temperature and pressure on chemical

Purpose: Demonstrate data visualization techniques
Target: Intermediate
Execution time: 2-5 seconds
Dependencies: None
"""

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
import seaborn as sns

# Experimental data: Temperature (3 levels) × Pressure (3 levels) = 9 experiments
# Temperature: 150°C, 175°C, 200°C
# Pressure: 1.0 MPa, 1.5 MPa, 2.0 MPa

np.random.seed(42)

# Yield data (%) (experimental order randomized)
data = pd.DataFrame({
    'Temperature': [150, 150, 150, 175, 175, 175, 200, 200, 200],
    'Pressure': [1.0, 1.5, 2.0, 1.0, 1.5, 2.0, 1.0, 1.5, 2.0],
    'Yield': [82.3, 85.1, 87.5, 86.2, 89.5, 91.3, 88.1, 90.8, 92.5]
})

print("=== Experimental Data ===")
print(data)

# Convert data to pivot table
pivot_data = data.pivot(index='Temperature', columns='Pressure', values='Yield')
print("\n=== Pivot Table (Yield %) ===")
print(pivot_data)

# Two-way ANOVA (without interaction)
# Using statsmodels
from statsmodels.formula.api import ols
from statsmodels.stats.anova import anova_lm

# Fit linear model
model = ols('Yield ~ C(Temperature) + C(Pressure)', data=data).fit()
anova_table = anova_lm(model, typ=2)

print("\n=== Two-way ANOVA ===")
print(anova_table)

# Visualize main effects
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Main effect of temperature
temp_means = data.groupby('Temperature')['Yield'].mean()
axes[0].plot(temp_means.index, temp_means.values, marker='o', linewidth=2, markersize=8, color='#11998e')
axes[0].set_xlabel('Temperature (°C)', fontsize=12)
axes[0].set_ylabel('Mean Yield (%)', fontsize=12)
axes[0].set_title('Main Effect of Temperature', fontsize=14, fontweight='bold')
axes[0].grid(alpha=0.3)

# Main effect of pressure
pressure_means = data.groupby('Pressure')['Yield'].mean()
axes[1].plot(pressure_means.index, pressure_means.values, marker='s', linewidth=2, markersize=8, color='#f59e0b')
axes[1].set_xlabel('Pressure (MPa)', fontsize=12)
axes[1].set_ylabel('Mean Yield (%)', fontsize=12)
axes[1].set_title('Main Effect of Pressure', fontsize=14, fontweight='bold')
axes[1].grid(alpha=0.3)

plt.tight_layout()
plt.savefig('two_way_anova_main_effects.png', dpi=300, bbox_inches='tight')
plt.show()

# Interaction plot
plt.figure(figsize=(10, 6))
for temp in [150, 175, 200]:
    subset = data[data['Temperature'] == temp]
    plt.plot(subset['Pressure'], subset['Yield'], marker='o', label=f'{temp}°C', linewidth=2, markersize=8)

plt.xlabel('Pressure (MPa)', fontsize=12)
plt.ylabel('Yield (%)', fontsize=12)
plt.title('Interaction Plot of Temperature and Pressure', fontsize=14, fontweight='bold')
plt.legend(title='Temperature', fontsize=10)
plt.grid(alpha=0.3)
plt.tight_layout()
plt.savefig('two_way_anova_interaction.png', dpi=300, bbox_inches='tight')
plt.show()

Example output:

=== Pivot Table (Yield %) ===
Pressure       1.0   1.5   2.0
Temperature
150           82.3  85.1  87.5
175           86.2  89.5  91.3
200           88.1  90.8  92.5

=== Two-way ANOVA ===
                    sum_sq   df         F    PR(>F)
C(Temperature)     57.3633  2.0  78.6275  0.000127
C(Pressure)        65.0433  2.0  89.1667  0.000094
Residual            1.4600  4.0       NaN       NaN

Conclusion: Both temperature and pressure significantly affect yield (p < 0.001)

Interpretation: Both temperature and pressure strongly influence yield, with higher temperature and pressure improving yield. The interaction plot shows nearly parallel lines, indicating minimal interaction.

1.3 Fundamentals of Orthogonal Arrays

What are Orthogonal Arrays

Orthogonal arrays are experimental design tables that efficiently arrange multiple factors. They are designed so that combinations of factor levels appear uniformly, enabling independent evaluation of factor effects with fewer experiments.

Major orthogonal arrays:

• L8 (2⁷): Up to 7 factors, 2 levels each, 8 experiments
• L9 (3⁴): Up to 4 factors, 3 levels each, 9 experiments
• L16 (2¹⁵): Up to 15 factors, 2 levels each, 16 experiments
• L27 (3¹³): Up to 13 factors, 3 levels each, 27 experiments

Code Example 3: Generation and Application of Orthogonal Array L8

Design an experiment with 3 factors (temperature, pressure, catalyst amount), 2 levels each, using orthogonal array L8.

# Requirements:
# - Python 3.9+
# - matplotlib>=3.7.0
# - numpy>=1.24.0, <2.0.0
# - pandas>=2.0.0, <2.2.0

"""
Example: Design an experiment with 3 factors (temperature, pressure, 

Purpose: Demonstrate data visualization techniques
Target: Intermediate
Execution time: 2-5 seconds
Dependencies: None
"""

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# Definition of orthogonal array L8 (2^7: 7 factors, 2 levels each, 8 experiments)
# Using 3 factors (columns 1, 2, 3)
L8 = np.array([
    [1, 1, 1, 1, 1, 1, 1],
    [1, 1, 1, 2, 2, 2, 2],
    [1, 2, 2, 1, 1, 2, 2],
    [1, 2, 2, 2, 2, 1, 1],
    [2, 1, 2, 1, 2, 1, 2],
    [2, 1, 2, 2, 1, 2, 1],
    [2, 2, 1, 1, 2, 2, 1],
    [2, 2, 1, 2, 1, 1, 2]
])

# Factor definition (using columns 1, 2, 3)
# Factor A: Temperature (level 1 = 150°C, level 2 = 200°C)
# Factor B: Pressure (level 1 = 1.0 MPa, level 2 = 2.0 MPa)
# Factor C: Catalyst amount (level 1 = 0.5 g, level 2 = 1.0 g)

factors = L8[:, :3]  # Use first 3 columns

# Map to actual values
temperature_levels = {1: 150, 2: 200}
pressure_levels = {1: 1.0, 2: 2.0}
catalyst_levels = {1: 0.5, 2: 1.0}

# Create experimental design table
doe_table = pd.DataFrame({
    'Run': range(1, 9),
    'Temperature': [temperature_levels[x] for x in factors[:, 0]],
    'Pressure': [pressure_levels[x] for x in factors[:, 1]],
    'Catalyst': [catalyst_levels[x] for x in factors[:, 2]]
})

print("=== Experimental Design with Orthogonal Array L8 ===")
print(doe_table)

# Simulated experimental results (Yield %)
# True model: Yield = 70 + 10*(Temp-150)/50 + 5*(Press-1) + 3*(Cat-0.5)/0.5 + noise
np.random.seed(42)

yields = []
for _, row in doe_table.iterrows():
    temp = row['Temperature']
    press = row['Pressure']
    cat = row['Catalyst']

    # True effect (linear model)
    yield_true = (70 +
                  10 * (temp - 150) / 50 +
                  5 * (press - 1.0) +
                  3 * (cat - 0.5) / 0.5)

    # Add noise
    yield_obs = yield_true + np.random.normal(0, 1)
    yields.append(yield_obs)

doe_table['Yield'] = yields

print("\n=== Experimental Results ===")
print(doe_table)

# Verify orthogonal array properties (equal appearance of each factor level)
print("\n=== Verification of Orthogonal Array Properties ===")
print("Count of each factor level:")
for col in ['Temperature', 'Pressure', 'Catalyst']:
    print(f"\n{col}:")
    print(doe_table[col].value_counts().sort_index())

Example output:

=== Experimental Design with Orthogonal Array L8 ===
   Run  Temperature  Pressure  Catalyst
0    1          150       1.0       0.5
1    2          150       1.0       1.0
2    3          150       2.0       0.5
3    4          150       2.0       1.0
4    5          200       1.0       0.5
5    6          200       1.0       1.0
6    7          200       2.0       0.5
7    8          200       2.0       1.0

=== Experimental Results ===
   Run  Temperature  Pressure  Catalyst      Yield
0    1          150       1.0       0.5  70.494371
1    2          150       1.0       1.0  75.861468
2    3          150       2.0       0.5  75.646968
3    4          150       2.0       1.0  79.522869
4    5          200       1.0       0.5  79.647689
5    6          200       1.0       1.0  85.522232
6    7          200       2.0       0.5  84.233257
7    8          200       2.0       1.0  88.767995

=== Verification of Orthogonal Array Properties ===
Count of each factor level:

Temperature:
150    4
200    4

Pressure:
1.0    4
2.0    4

Catalyst:
0.5    4
1.0    4

Interpretation: Orthogonal array L8 evaluates 3 factors in 8 experiments. Each level of each factor appears uniformly 4 times, enabling independent evaluation of factor effects.

Code Example 4: Multi-factor Experiment with Orthogonal Array L16

Evaluate 5 factors (temperature, pressure, catalyst amount, reaction time, stirring speed) using orthogonal array L16.

# Requirements:
# - Python 3.9+
# - numpy>=1.24.0, <2.0.0
# - pandas>=2.0.0, <2.2.0

import numpy as np
import pandas as pd

# Generation of orthogonal array L16 (2^15: up to 15 factors, 2 levels each, 16 experiments)
# Using 5 factors
def generate_L16():
    """Generate orthogonal array L16"""
    L16 = []
    for i in range(16):
        row = []
        for j in range(15):
            # Generate 2 levels (1 or 2)
            val = ((i >> j) & 1) + 1
            row.append(val)
        L16.append(row)
    return np.array(L16)

L16 = generate_L16()

# Assign 5 factors to columns 1, 2, 4, 8, 15 (standard assignment)
factor_columns = [0, 1, 3, 7, 14]  # Python indices
factors = L16[:, factor_columns]

# Factor definition
factor_names = ['Temperature', 'Pressure', 'Catalyst', 'Time', 'Stirring']
levels = {
    'Temperature': {1: 150, 2: 200},
    'Pressure': {1: 1.0, 2: 2.0},
    'Catalyst': {1: 0.5, 2: 1.0},
    'Time': {1: 30, 2: 60},       # Reaction time (minutes)
    'Stirring': {1: 200, 2: 400}  # Stirring speed (rpm)
}

# Create experimental design table
doe_table = pd.DataFrame({'Run': range(1, 17)})

for i, fname in enumerate(factor_names):
    doe_table[fname] = [levels[fname][x] for x in factors[:, i]]

print("=== 5-Factor Experimental Design with Orthogonal Array L16 ===")
print(doe_table.to_string(index=False))

# Simulated experimental results
np.random.seed(42)

yields = []
for _, row in doe_table.iterrows():
    temp = row['Temperature']
    press = row['Pressure']
    cat = row['Catalyst']
    time = row['Time']
    stir = row['Stirring']

    # True model (main effects only)
    yield_true = (60 +
                  8 * (temp - 150) / 50 +
                  4 * (press - 1.0) +
                  3 * (cat - 0.5) / 0.5 +
                  2 * (time - 30) / 30 +
                  1 * (stir - 200) / 200)

    # Add noise
    yield_obs = yield_true + np.random.normal(0, 1.5)
    yields.append(yield_obs)

doe_table['Yield'] = yields

print("\n=== Experimental Results (First 5 experiments) ===")
print(doe_table.head())

print(f"\nTotal experiments: {len(doe_table)}")
print(f"Number of factors evaluated: {len(factor_names)}")
print(f"Efficiency: Can be conducted at 50% of full factorial design (2^5=32 runs)")

Example output:

=== 5-Factor Experimental Design with Orthogonal Array L16 ===
 Run  Temperature  Pressure  Catalyst  Time  Stirring
   1          150       1.0       0.5    30       200
   2          200       1.0       0.5    30       200
   3          150       2.0       0.5    30       200
   4          200       2.0       0.5    30       200
   5          150       1.0       1.0    30       200
...(omitted)

=== Experimental Results (First 5 experiments) ===
   Run  Temperature  Pressure  Catalyst  Time  Stirring      Yield
0    1          150       1.0       0.5    30       200  59.494371
1    2          200       1.0       0.5    30       200  67.861468
2    3          150       2.0       0.5    30       200  63.646968
3    4          200       2.0       0.5    30       200  71.522869
4    5          150       1.0       1.0    30       200  65.647689

Total experiments: 16
Number of factors evaluated: 5
Efficiency: Can be conducted at 50% of full factorial design (2^5=32 runs)

Interpretation: Using orthogonal array L16, evaluation of 5 factors can be conducted in 16 experiments (half of full factorial design).

1.4 Main Effect Plots and Interaction Plots

Code Example 5: Visualization of Interactions

Analyze the interaction between temperature and pressure in detail.

# Requirements:
# - Python 3.9+
# - matplotlib>=3.7.0
# - numpy>=1.24.0, <2.0.0
# - pandas>=2.0.0, <2.2.0

"""
Example: Analyze the interaction between temperature and pressure in 

Purpose: Demonstrate data visualization techniques
Target: Intermediate
Execution time: 2-5 seconds
Dependencies: None
"""

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# Experimental data: Model with temperature×pressure interaction
np.random.seed(42)

temperatures = [150, 175, 200]
pressures = [1.0, 1.5, 2.0]

data = []
for temp in temperatures:
    for press in pressures:
        # Model with interaction term
        # Yield = β0 + β1*Temp + β2*Press + β3*Temp*Press + ε
        yield_true = (50 +
                      0.15 * temp +
                      20 * press +
                      0.05 * temp * press)  # Interaction term

        yield_obs = yield_true + np.random.normal(0, 1.5)
        data.append({
            'Temperature': temp,
            'Pressure': press,
            'Yield': yield_obs
        })

df = pd.DataFrame(data)

print("=== Experimental Data ===")
print(df)

# Interaction plots
plt.figure(figsize=(12, 5))

# Left: Fix temperature and observe pressure effect
plt.subplot(1, 2, 1)
for temp in temperatures:
    subset = df[df['Temperature'] == temp]
    plt.plot(subset['Pressure'], subset['Yield'],
             marker='o', linewidth=2, markersize=8, label=f'{temp}°C')

plt.xlabel('Pressure (MPa)', fontsize=12)
plt.ylabel('Yield (%)', fontsize=12)
plt.title('Pressure Effect by Temperature (with interaction)', fontsize=14, fontweight='bold')
plt.legend(title='Temperature', fontsize=10)
plt.grid(alpha=0.3)

# Right: Fix pressure and observe temperature effect
plt.subplot(1, 2, 2)
for press in pressures:
    subset = df[df['Pressure'] == press]
    plt.plot(subset['Temperature'], subset['Yield'],
             marker='s', linewidth=2, markersize=8, label=f'{press} MPa')

plt.xlabel('Temperature (°C)', fontsize=12)
plt.ylabel('Yield (%)', fontsize=12)
plt.title('Temperature Effect by Pressure (with interaction)', fontsize=14, fontweight='bold')
plt.legend(title='Pressure', fontsize=10)
plt.grid(alpha=0.3)

plt.tight_layout()
plt.savefig('interaction_plot.png', dpi=300, bbox_inches='tight')
plt.show()

# Quantitative evaluation of interaction
print("\n=== Evaluation of Interaction ===")
print("Parallel lines → No interaction")
print("Crossing lines → Strong interaction")
print("\nThis data: Changing slopes → Interaction exists between temperature and pressure")

Interpretation: When lines are not parallel and slopes change, interaction exists. In this case, a synergistic effect is observed in the high temperature × high pressure combination.

Code Example 6: Creating Main Effect Plots

Create main effect plots from orthogonal array L8 results.

# Requirements:
# - Python 3.9+
# - matplotlib>=3.7.0
# - numpy>=1.24.0, <2.0.0
# - pandas>=2.0.0, <2.2.0

"""
Example: Create main effect plots from orthogonal array L8 results.

Purpose: Demonstrate data visualization techniques
Target: Intermediate
Execution time: 2-5 seconds
Dependencies: None
"""

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# Using L8 experimental data (doe_table and yields from Code Example 3)
# Regenerating here
np.random.seed(42)

# Orthogonal array L8
L8 = np.array([
    [1, 1, 1],
    [1, 1, 2],
    [1, 2, 1],
    [1, 2, 2],
    [2, 1, 1],
    [2, 1, 2],
    [2, 2, 1],
    [2, 2, 2]
])

temperature_levels = {1: 150, 2: 200}
pressure_levels = {1: 1.0, 2: 2.0}
catalyst_levels = {1: 0.5, 2: 1.0}

doe_table = pd.DataFrame({
    'Run': range(1, 9),
    'Temp_code': L8[:, 0],
    'Press_code': L8[:, 1],
    'Cat_code': L8[:, 2],
    'Temperature': [temperature_levels[x] for x in L8[:, 0]],
    'Pressure': [pressure_levels[x] for x in L8[:, 1]],
    'Catalyst': [catalyst_levels[x] for x in L8[:, 2]]
})

# Simulated yield
yields = []
for _, row in doe_table.iterrows():
    yield_true = (70 +
                  10 * (row['Temperature'] - 150) / 50 +
                  5 * (row['Pressure'] - 1.0) +
                  3 * (row['Catalyst'] - 0.5) / 0.5)
    yield_obs = yield_true + np.random.normal(0, 1)
    yields.append(yield_obs)

doe_table['Yield'] = yields

# Calculate main effects (average yield at each level of each factor)
main_effects = {}

# Main effect of temperature
temp_level1 = doe_table[doe_table['Temp_code'] == 1]['Yield'].mean()
temp_level2 = doe_table[doe_table['Temp_code'] == 2]['Yield'].mean()
main_effects['Temperature'] = {150: temp_level1, 200: temp_level2}

# Main effect of pressure
press_level1 = doe_table[doe_table['Press_code'] == 1]['Yield'].mean()
press_level2 = doe_table[doe_table['Press_code'] == 2]['Yield'].mean()
main_effects['Pressure'] = {1.0: press_level1, 2.0: press_level2}

# Main effect of catalyst amount
cat_level1 = doe_table[doe_table['Cat_code'] == 1]['Yield'].mean()
cat_level2 = doe_table[doe_table['Cat_code'] == 2]['Yield'].mean()
main_effects['Catalyst'] = {0.5: cat_level1, 1.0: cat_level2}

# Create main effect plots
fig, axes = plt.subplots(1, 3, figsize=(16, 5))

# Main effect of temperature
axes[0].plot([150, 200], [temp_level1, temp_level2],
             marker='o', linewidth=2.5, markersize=10, color='#11998e')
axes[0].axhline(y=doe_table['Yield'].mean(), color='red', linestyle='--',
                linewidth=1.5, label='Overall mean')
axes[0].set_xlabel('Temperature (°C)', fontsize=12)
axes[0].set_ylabel('Mean Yield (%)', fontsize=12)
axes[0].set_title('Main Effect of Temperature', fontsize=14, fontweight='bold')
axes[0].legend()
axes[0].grid(alpha=0.3)

# Main effect of pressure
axes[1].plot([1.0, 2.0], [press_level1, press_level2],
             marker='s', linewidth=2.5, markersize=10, color='#f59e0b')
axes[1].axhline(y=doe_table['Yield'].mean(), color='red', linestyle='--',
                linewidth=1.5, label='Overall mean')
axes[1].set_xlabel('Pressure (MPa)', fontsize=12)
axes[1].set_ylabel('Mean Yield (%)', fontsize=12)
axes[1].set_title('Main Effect of Pressure', fontsize=14, fontweight='bold')
axes[1].legend()
axes[1].grid(alpha=0.3)

# Main effect of catalyst amount
axes[2].plot([0.5, 1.0], [cat_level1, cat_level2],
             marker='^', linewidth=2.5, markersize=10, color='#7b2cbf')
axes[2].axhline(y=doe_table['Yield'].mean(), color='red', linestyle='--',
                linewidth=1.5, label='Overall mean')
axes[2].set_xlabel('Catalyst Amount (g)', fontsize=12)
axes[2].set_ylabel('Mean Yield (%)', fontsize=12)
axes[2].set_title('Main Effect of Catalyst Amount', fontsize=14, fontweight='bold')
axes[2].legend()
axes[2].grid(alpha=0.3)

plt.tight_layout()
plt.savefig('main_effects_plot.png', dpi=300, bbox_inches='tight')
plt.show()

# Calculate magnitude of main effects (effect = high level mean - low level mean)
print("=== Magnitude of Main Effects ===")
print(f"Temperature effect: {temp_level2 - temp_level1:.2f} %")
print(f"Pressure effect: {press_level2 - press_level1:.2f} %")
print(f"Catalyst amount effect: {cat_level2 - cat_level1:.2f} %")

print("\n=== Factor Importance Ranking ===")
effects = {
    'Temperature': abs(temp_level2 - temp_level1),
    'Pressure': abs(press_level2 - press_level1),
    'Catalyst amount': abs(cat_level2 - cat_level1)
}
sorted_effects = sorted(effects.items(), key=lambda x: x[1], reverse=True)
for i, (factor, effect) in enumerate(sorted_effects, 1):
    print(f"{i}. {factor} (effect: {effect:.2f} %)")

Example output:

=== Magnitude of Main Effects ===
Temperature effect: 10.15 %
Pressure effect: 5.08 %
Catalyst amount effect: 3.04 %

=== Factor Importance Ranking ===
1. Temperature (effect: 10.15 %)
2. Pressure (effect: 5.08 %)
3. Catalyst amount (effect: 3.04 %)

Interpretation: Main effect plots show that temperature has the largest influence, followed by pressure and catalyst amount. Optimal conditions are temperature 200°C, pressure 2.0 MPa, catalyst amount 1.0 g.

Code Example 7: Creating Interaction Plots

Analyze the temperature×pressure interaction in detail.

# Requirements:
# - Python 3.9+
# - matplotlib>=3.7.0
# - numpy>=1.24.0, <2.0.0
# - pandas>=2.0.0, <2.2.0

"""
Example: Analyze the temperature×pressure interaction in detail.

Purpose: Demonstrate data visualization techniques
Target: Intermediate
Execution time: 2-5 seconds
Dependencies: None
"""

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# Full factorial design for 2 factors (to evaluate interaction)
np.random.seed(42)

data = []
for temp in [150, 200]:
    for press in [1.0, 2.0]:
        # Model with interaction term
        # Yield = β0 + β1*Temp + β2*Press + β3*Temp*Press + ε
        # Intentionally adding interaction term here
        yield_true = (50 +
                      0.10 * temp +
                      15 * press +
                      0.03 * temp * press)  # Interaction term

        # 3 replicates
        for rep in range(3):
            yield_obs = yield_true + np.random.normal(0, 1.0)
            data.append({
                'Temperature': temp,
                'Pressure': press,
                'Replicate': rep + 1,
                'Yield': yield_obs
            })

df = pd.DataFrame(data)

# Calculate mean for each condition
avg_df = df.groupby(['Temperature', 'Pressure'])['Yield'].mean().reset_index()

print("=== Mean Yield for Each Condition ===")
print(avg_df)

# Create interaction plot
plt.figure(figsize=(10, 6))

for temp in [150, 200]:
    subset = avg_df[avg_df['Temperature'] == temp]
    plt.plot(subset['Pressure'], subset['Yield'],
             marker='o', linewidth=2.5, markersize=10, label=f'{temp}°C')

plt.xlabel('Pressure (MPa)', fontsize=12)
plt.ylabel('Mean Yield (%)', fontsize=12)
plt.title('Interaction Plot of Temperature×Pressure', fontsize=14, fontweight='bold')
plt.legend(title='Temperature', fontsize=11)
plt.xticks([1.0, 2.0])
plt.grid(alpha=0.3)
plt.tight_layout()
plt.savefig('interaction_plot_detailed.png', dpi=300, bbox_inches='tight')
plt.show()

# Quantitative evaluation of interaction
# Interaction = (high temp, high press - high temp, low press) - (low temp, high press - low temp, low press)
y_150_10 = avg_df[(avg_df['Temperature'] == 150) & (avg_df['Pressure'] == 1.0)]['Yield'].values[0]
y_150_20 = avg_df[(avg_df['Temperature'] == 150) & (avg_df['Pressure'] == 2.0)]['Yield'].values[0]
y_200_10 = avg_df[(avg_df['Temperature'] == 200) & (avg_df['Pressure'] == 1.0)]['Yield'].values[0]
y_200_20 = avg_df[(avg_df['Temperature'] == 200) & (avg_df['Pressure'] == 2.0)]['Yield'].values[0]

interaction = (y_200_20 - y_200_10) - (y_150_20 - y_150_10)

print(f"\n=== Magnitude of Interaction ===")
print(f"Temperature×Pressure interaction: {interaction:.2f} %")

if abs(interaction) > 2:
    print("Judgment: Interaction exists (|interaction| > 2%)")
else:
    print("Judgment: Interaction is negligible (|interaction| ≤ 2%)")

print("\n=== Interpretation ===")
if interaction > 0:
    print("Synergistic effect observed in high temperature×high pressure combination")
else:
    print("Negative interaction (effect decreases with combination of high levels)")

Example output:

=== Mean Yield for Each Condition ===
   Temperature  Pressure      Yield
0          150       1.0  80.329885
1          150       2.0  95.246314
2          200       1.0  85.294928
3          200       2.0 106.022349

=== Magnitude of Interaction ===
Temperature×Pressure interaction: 5.81 %

Judgment: Interaction exists (|interaction| > 2%)

=== Interpretation ===
Synergistic effect observed in high temperature×high pressure combination

Interpretation: When interaction plot lines are not parallel (crossing or different slopes), interaction exists. In this case, combining high temperature and high pressure produces a synergistic effect.

1.5 Case Study: Optimization of Chemical Reaction Yield

Code Example 8: 3-Factor Experiment with Orthogonal Array L8 and Optimal Condition Search

Optimize temperature, catalyst concentration, and reaction time in an esterification reaction using orthogonal array L8.

# Requirements:
# - Python 3.9+
# - matplotlib>=3.7.0
# - numpy>=1.24.0, <2.0.0
# - pandas>=2.0.0, <2.2.0
# - seaborn>=0.12.0

"""
Example: Optimize temperature, catalyst concentration, and reaction t

Purpose: Demonstrate data visualization techniques
Target: Advanced
Execution time: 2-5 seconds
Dependencies: None
"""

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

# Case study: Optimization of esterification reaction yield
# Factor A: Reaction temperature (60°C vs 80°C)
# Factor B: Catalyst concentration (0.1 M vs 0.5 M)
# Factor C: Reaction time (2 hours vs 4 hours)

np.random.seed(42)

# Orthogonal array L8
L8 = np.array([
    [1, 1, 1],
    [1, 1, 2],
    [1, 2, 1],
    [1, 2, 2],
    [2, 1, 1],
    [2, 1, 2],
    [2, 2, 1],
    [2, 2, 2]
])

# Factor definition
factor_levels = {
    'Temperature': {1: 60, 2: 80},
    'Catalyst': {1: 0.1, 2: 0.5},
    'Time': {1: 2, 2: 4}
}

# Create experimental design table
doe_table = pd.DataFrame({
    'Run': range(1, 9),
    'Temp_code': L8[:, 0],
    'Cat_code': L8[:, 1],
    'Time_code': L8[:, 2],
    'Temperature': [factor_levels['Temperature'][x] for x in L8[:, 0]],
    'Catalyst': [factor_levels['Catalyst'][x] for x in L8[:, 1]],
    'Time': [factor_levels['Time'][x] for x in L8[:, 2]]
})

# Simulated yield (realistic model)
# True model: Yield = f(Temp, Cat, Time) + interaction + noise
yields = []
for _, row in doe_table.iterrows():
    temp = row['Temperature']
    cat = row['Catalyst']
    time = row['Time']

    # Main effects
    yield_base = 50
    temp_effect = 15 * (temp - 60) / 20
    cat_effect = 20 * (cat - 0.1) / 0.4
    time_effect = 10 * (time - 2) / 2

    # Interaction (temperature×catalyst concentration)
    interaction = 5 * ((temp - 60) / 20) * ((cat - 0.1) / 0.4)

    yield_true = yield_base + temp_effect + cat_effect + time_effect + interaction

    # Add noise
    yield_obs = yield_true + np.random.normal(0, 2)
    yields.append(yield_obs)

doe_table['Yield'] = yields

print("=== Esterification Reaction Experimental Design and Results ===")
print(doe_table[['Run', 'Temperature', 'Catalyst', 'Time', 'Yield']])

# Main effects analysis
print("\n=== Main Effects Analysis ===")

# Main effect of temperature
temp_low = doe_table[doe_table['Temp_code'] == 1]['Yield'].mean()
temp_high = doe_table[doe_table['Temp_code'] == 2]['Yield'].mean()
temp_effect = temp_high - temp_low
print(f"Temperature effect: {temp_effect:.2f}% (low level: {temp_low:.2f}%, high level: {temp_high:.2f}%)")

# Main effect of catalyst concentration
cat_low = doe_table[doe_table['Cat_code'] == 1]['Yield'].mean()
cat_high = doe_table[doe_table['Cat_code'] == 2]['Yield'].mean()
cat_effect = cat_high - cat_low
print(f"Catalyst concentration effect: {cat_effect:.2f}% (low level: {cat_low:.2f}%, high level: {cat_high:.2f}%)")

# Main effect of reaction time
time_low = doe_table[doe_table['Time_code'] == 1]['Yield'].mean()
time_high = doe_table[doe_table['Time_code'] == 2]['Yield'].mean()
time_effect = time_high - time_low
print(f"Reaction time effect: {time_effect:.2f}% (low level: {time_low:.2f}%, high level: {time_high:.2f}%)")

# Factor importance ranking
effects = {
    'Catalyst concentration': abs(cat_effect),
    'Temperature': abs(temp_effect),
    'Reaction time': abs(time_effect)
}
sorted_effects = sorted(effects.items(), key=lambda x: x[1], reverse=True)

print("\n=== Factor Importance Ranking ===")
for i, (factor, effect) in enumerate(sorted_effects, 1):
    print(f"{i}. {factor} (effect: {effect:.2f}%)")

# Determine optimal conditions
print("\n=== Optimal Conditions ===")
print("Conditions to maximize yield:")
print(f"  Temperature: {80 if temp_effect > 0 else 60}°C")
print(f"  Catalyst concentration: {0.5 if cat_effect > 0 else 0.1} M")
print(f"  Reaction time: {4 if time_effect > 0 else 2} hours")

# Predicted yield (optimal conditions)
predicted_yield_max = doe_table['Yield'].mean() + abs(temp_effect)/2 + abs(cat_effect)/2 + abs(time_effect)/2
print(f"  Predicted yield: {predicted_yield_max:.1f}%")

# Visualize main effects
fig, axes = plt.subplots(1, 3, figsize=(16, 5))

# Main effect of temperature
axes[0].plot([60, 80], [temp_low, temp_high],
             marker='o', linewidth=2.5, markersize=10, color='#11998e')
axes[0].axhline(y=doe_table['Yield'].mean(), color='red', linestyle='--',
                linewidth=1.5, label='Overall mean', alpha=0.7)
axes[0].set_xlabel('Reaction Temperature (°C)', fontsize=12)
axes[0].set_ylabel('Mean Yield (%)', fontsize=12)
axes[0].set_title('Main Effect of Temperature', fontsize=14, fontweight='bold')
axes[0].legend()
axes[0].grid(alpha=0.3)

# Main effect of catalyst concentration
axes[1].plot([0.1, 0.5], [cat_low, cat_high],
             marker='s', linewidth=2.5, markersize=10, color='#f59e0b')
axes[1].axhline(y=doe_table['Yield'].mean(), color='red', linestyle='--',
                linewidth=1.5, label='Overall mean', alpha=0.7)
axes[1].set_xlabel('Catalyst Concentration (M)', fontsize=12)
axes[1].set_ylabel('Mean Yield (%)', fontsize=12)
axes[1].set_title('Main Effect of Catalyst Concentration', fontsize=14, fontweight='bold')
axes[1].legend()
axes[1].grid(alpha=0.3)

# Main effect of reaction time
axes[2].plot([2, 4], [time_low, time_high],
             marker='^', linewidth=2.5, markersize=10, color='#7b2cbf')
axes[2].axhline(y=doe_table['Yield'].mean(), color='red', linestyle='--',
                linewidth=1.5, label='Overall mean', alpha=0.7)
axes[2].set_xlabel('Reaction Time (hours)', fontsize=12)
axes[2].set_ylabel('Mean Yield (%)', fontsize=12)
axes[2].set_title('Main Effect of Reaction Time', fontsize=14, fontweight='bold')
axes[2].legend()
axes[2].grid(alpha=0.3)

plt.tight_layout()
plt.savefig('esterification_main_effects.png', dpi=300, bbox_inches='tight')
plt.show()

# Visualize results with heatmap
pivot_temp_cat = doe_table.pivot_table(values='Yield',
                                        index='Temperature',
                                        columns='Catalyst',
                                        aggfunc='mean')

plt.figure(figsize=(8, 6))
sns.heatmap(pivot_temp_cat, annot=True, fmt='.1f', cmap='RdYlGn',
            cbar_kws={'label': 'Yield (%)'}, linewidths=2, linecolor='white')
plt.title('Yield Map of Temperature×Catalyst Concentration', fontsize=14, fontweight='bold')
plt.ylabel('Temperature (°C)', fontsize=12)
plt.xlabel('Catalyst Concentration (M)', fontsize=12)
plt.tight_layout()
plt.savefig('esterification_heatmap.png', dpi=300, bbox_inches='tight')
plt.show()

print("\n=== Case Study Summary ===")
print("✅ Evaluated 3 factors with 8 experiments using orthogonal array L8")
print("✅ Identified catalyst concentration as the most important factor")
print("✅ Optimal conditions: 80°C, 0.5 M, 4 hours")
print(f"✅ Predicted yield: {predicted_yield_max:.1f}% (confirmation experiment recommended)")

Example output:

=== Esterification Reaction Experimental Design and Results ===
   Run  Temperature  Catalyst  Time      Yield
0    1           60       0.1     2  50.987420
1    2           60       0.1     4  60.723869
2    3           60       0.5     2  70.294968
3    4           60       0.5     4  80.046738
4    5           80       0.1     2  65.296378
5    6           80       0.1     4  75.044464
6    7           80       0.5     2  90.466514
7    8           80       0.5     4 100.535989

=== Main Effects Analysis ===
Temperature effect: 15.09% (low level: 65.51%, high level: 80.59%)
Catalyst concentration effect: 20.25% (low level: 63.01%, high level: 83.26%)
Reaction time effect: 9.99% (low level: 68.14%, high level: 78.13%)

=== Factor Importance Ranking ===
1. Catalyst concentration (effect: 20.25%)
2. Temperature (effect: 15.09%)
3. Reaction time (effect: 9.99%)

=== Optimal Conditions ===
Conditions to maximize yield:
  Temperature: 80°C
  Catalyst concentration: 0.5 M
  Reaction time: 4 hours
  Predicted yield: 95.7%

=== Case Study Summary ===
✅ Evaluated 3 factors with 8 experiments using orthogonal array L8
✅ Identified catalyst concentration as the most important factor
✅ Optimal conditions: 80°C, 0.5 M, 4 hours
✅ Predicted yield: 95.7% (confirmation experiment recommended)

Interpretation: Using orthogonal array L8, we efficiently evaluated the effects of 3 factors in only 8 experiments and determined optimal conditions. Conventional one-variable-at-a-time experiments would require at least 3×2×2×2=24 runs, but DOE reduced experimental runs by 67%.

1.6 Chapter Summary

What We Learned

1. Fundamentals of DOE — A statistical method to obtain maximum information with minimum experiments, its differences from conventional OFAT (one-variable-at-a-time) experiments, and the three principles of DOE: replication, randomization, and blocking.

2. One-way and Two-way Experiments — One-way ANOVA for comparison between levels of a single factor, two-way ANOVA for main effects and interactions of two factors, and determination of statistical significance by F-test.

3. Application of Orthogonal Arrays — Efficient experimental design with orthogonal arrays L8 and L16, the property of equal appearance of each factor level, and significant reduction in experimental runs (50-75%).

4. Main Effect Plots and Interaction Plots — Main effects for visualization of individual factor effects, interactions for detection of mutual effects between factors, and methods to determine optimal conditions.

5. Application to Chemical Processes — Optimization of esterification reaction yield, identification of important factors and search for optimal conditions, and calculation of predicted yield with necessity of confirmation experiments.

Important Points

DOE can reduce experimental runs by 50-75% while enabling evaluation of interactions. Orthogonal arrays are optimal for factor screening to identify important factors. Main effect plots enable visual understanding of factor influences. When interaction plot lines are not parallel, interaction exists. Optimal conditions are combinations of levels that maximize main effects.

To the Next Chapter

In Chapter 2, we will learn about Factorial Designs and ANOVA in detail, covering full factorial design, fractional factorial design, details of ANOVA and F-test, multiple comparison tests (Tukey HSD), decomposition of variance components and calculation of contribution ratios, and a case study exploring factors affecting catalyst activity.