Coding, scaling, and the linear model#

Source worksheet: yint.org/w2 - week 2 of the applied DoE short course.

Module 1 read effects straight off corner averages of a 2x2. That was appropriate for two factors and four runs, but the same approach becomes cumbersome once you scale beyond it. Module 2 introduces the coded linear model

  • the workhorse that scales to any number of factors, copes with interactions, and gives you predictions outside the corners you ran.

The math is unchanged from Module 1; what we add is coded units (-1 and +1 instead of grams and minutes), a linear regression that recovers the same effects you computed by hand, and a sober look at extrapolation - what the model promises, and what it cannot promise.

Everything below uses Python with process_improve; the original worksheet’s R snippet at yint.org/2w1 is translated to its Python equivalent.

Tip

The first habit to build: in design space, every factor lives on the same scale, with -1 at its low extreme and +1 at its high extreme. A 30 mL setting becomes +1; 200 grams becomes -1. This is coding. It makes coefficients directly comparable and makes predictions outside the corner runs straightforward to write down. In process_improve you state the real-world range with c(..., lo=..., hi=...) and the package handles the coding.

Questions 1 and 2 - how many runs does a full factorial need?#

A full factorial runs every combination of every factor level. For k factors at L levels each, that is L^k experiments.

  1. Four factors, each at two levels: how many runs?

  2. Four factors, each at three levels: how many runs?

[1]:

def full_factorial_size(num_factors: int, levels: int) -> int:
    return levels ** num_factors


two_level = full_factorial_size(num_factors=4, levels=2)
three_level = full_factorial_size(num_factors=4, levels=3)
print(f"Full factorial, 4 factors x 2 levels: {two_level} runs")
print(f"Full factorial, 4 factors x 3 levels: {three_level} runs")

Full factorial, 4 factors x 2 levels: 16 runs
Full factorial, 4 factors x 3 levels: 81 runs

Solution

  1. 2**4 = 16 runs.

  2. 3**4 = 81 runs.

The rapid growth in run count motivates fractional factorials (Module 5): even modest screens become unaffordable at three or four levels.

Guidance

Two levels per factor is the default starting point because four factors fit in 16 runs, eight factors fit in 256 - but two levels already test the linear effect of every factor. A third level (often a center point at 0) is added when you suspect curvature; we revisit this in Module 4.

Questions 3 to 9 - greenhouse plant heights#

A greenhouse experiment varies two factors to maximize plant height:

  • A = soil amount [200 g or 400 g]

  • B = water amount [50 mL or 100 mL]

The runs come out of the lab notebook in experiment order, not in the textbook’s standard order:

Run

A (g)

B (mL)

Height (mm)

1

200

100

61

2

200

50

39

3

400

100

82

4

400

50

50

The worksheet asks for the rewrite in standard order, the main effects of A and B, the prediction equation, and three predicted heights at new conditions.

Check yourself

Before reading on, write the four rows in standard order (A = -, -, +, + and B = -, +, -, +) on paper and copy the four heights across. Standard order does not change the data, only the row order; it is the layout the model fitter expects and the easiest layout to spot patterns in.

[2]:

import plotly.graph_objects as go
import plotly.io as pio

from process_improve.experiments import c, gather, lm, main_effects_plot, predict
from process_improve.experiments.visualization import visualize_doe

pio.renderers.default = "notebook_connected"

# Build the design directly in standard order: A varies fastest, then B.

A = c(-1, +1, -1, +1, lo=200, hi=400, name="A", units="grams")
B = c(-1, -1, +1, +1, lo=50,  hi=100, name="B", units="mL")
height = c(39, 50, 61, 82, name="height")
plant = gather(A=A, B=B, height=height)
plant

[2]:
A B height
1 -1.0 -1.0 39.0
2 1.0 -1.0 50.0
3 -1.0 1.0 61.0
4 1.0 1.0 82.0 
[3]:

# Fit a main-effects model in coded units.

model = lm("height ~ A + B", plant)
params = model.get_parameters(drop_intercept=False)
print(params.to_string())

effect_A = 2 * params["A"]
effect_B = 2 * params["B"]
print(f"\nMain effect of soil amount (A):  {effect_A:+.1f} mm")
print(f"Main effect of water amount (B): {effect_B:+.1f} mm")

Intercept    58.0
A             8.0
B            13.5

Main effect of soil amount (A):  +16.0 mm
Main effect of water amount (B): +27.0 mm

[4]:

# Verify the prediction equation by predicting at three test points.

points = [
    ("300 g soil, 75 mL water", dict(A=0,  B=0)),
    ("200 g soil, 75 mL water", dict(A=-1, B=0)),
    ("500 g soil, 125 mL water (extrapolation!)", dict(A=2, B=2)),
]
for label, coded in points:
    y_hat = float(predict(model, **coded).iloc[0])
    print(f"{label:50s} -> {y_hat:6.1f} mm")

300 g soil, 75 mL water                            ->   58.0 mm
200 g soil, 75 mL water                            ->   50.0 mm
500 g soil, 125 mL water (extrapolation!)          ->  101.0 mm

[5]:

fig = main_effects_plot(model, factor_labels={"A": "Soil [g]", "B": "Water [mL]"})
fig.update_layout(width=620, height=380, title="Main-effects plot - greenhouse")
fig

Solution

Standard order and effects. Sorted (A, B) = (-,-)=39, (+,-)=50, (-,+)=61, (+,+)=82.

  • Main effect of soil: (50 + 82)/2 - (39 + 61)/2 = 66 - 50 = +16 mm.

  • Main effect of water: (61 + 82)/2 - (39 + 50)/2 = 71.5 - 44.5 = +27 mm.

Coefficients. The intercept is the grand mean, b_0 = (39 + 50 + 61 + 82)/4 = 58. Each coefficient is half its main effect, so b_A = +8 and b_B = +13.5.

Prediction equation. In coded units, y = 58 + 8 * x_A + 13.5 * x_B.

Three predictions.

  • 300 g soil and 75 mL water sit at the center of the design (x_A = 0, x_B = 0); the model returns the intercept, y = 58 mm.

  • 200 g soil and 75 mL water moves A to -1 but keeps B at 0; predicted height y = 58 - 8 = 50 mm.

  • 500 g soil and 125 mL water sits at x_A = +2, x_B = +2 - one whole step outside the design on both axes. The model extrapolates to y = 58 + 16 + 27 = 101 mm; in real life plants do not grow proportionally forever.

Questions 10 to 12 - extrapolation, one-factor-at-a-time, and interactions#

The worksheet now steps back from the maths:

  1. Can you extrapolate outside the design region? How far? And what are contour lines useful for?

  2. A colleague wants to run their next study by changing one factor at a time against a baseline. Why is that a bad idea?

  3. Give an everyday example of an interaction where increasing factor B:

    • amplifies the effect of factor A;

    • cancels or reverses the effect of factor A.

Solution

Q10 - Extrapolation. A linear model lets you compute a prediction at any (x_A, x_B); it does not promise that prediction is right. Inside the four corners you ran, the model interpolates between observed responses and is usually trustworthy. Just outside, predictions are still useful as a direction: they tell you which way to push for the next experiment. Far outside, expect the world to push back - biological limits, melting points, saturated reactions.

Contour plots are the most useful summary view: they show you how to trade off one factor against another while holding the response constant. A flat contour is a robust operating window; crowded contours are where the process is sensitive.

Q11 - One factor at a time (OFAT). Three failures of OFAT compared with a factorial:

  1. It misses interactions completely. With OFAT you tweak A while holding B at its baseline, then tweak B while holding A at its baseline; you never observe the A = +, B = + corner where the interaction lives.

  2. It uses runs less efficiently. A 2x2 spends every run estimating both main effects; OFAT spends each run estimating just one.

  3. The “winning” combination depends on which factor you chose to vary first. That conclusion is fragile, and it is hard to defend when a colleague asks “what if we had started with B?”.

Q12 - Everyday interactions.

  • Amplification. Coffee and sleep deprivation: the same dose of caffeine wakes you up much more after a short night than after a long one. A and B amplify each other.

  • Cancellation. Salt and soup volume: a teaspoon of salt is perfect in one mug of soup and unnoticeable in a stock-pot. The effect of “adding salt” depends on (and is cancelled by) the dilution from extra water.

Question 13 - a baked food with a negative interaction#

A product is baked for time T using a certain amount of fat F. The outcome is crispiness. Both factors have plausible main effects, but the data hide a negative interaction that turns up only when you fit the model.

Run

T (min)

F (g)

Crispiness

1

80

20

37

2

100

20

57

3

80

30

49

4

100

30

53

The worksheet asks for the cube plot with contours, an interaction plot in both orientations, a verification of the prediction equation, and two extrapolated crispness predictions.

[6]:

T = c(-1, +1, -1, +1, lo=80, hi=100, name="T", units="min")
F = c(-1, -1, +1, +1, lo=20, hi=30,  name="F", units="g")
crisp = c(37, 57, 49, 53, name="crisp")
bake = gather(T=T, F=F, crisp=crisp)
bake

[6]:
T F crisp
1 -1.0 -1.0 37.0
2 1.0 -1.0 57.0
3 -1.0 1.0 49.0
4 1.0 1.0 53.0 
[7]:

square = visualize_doe(
    plot_type="square_plot",
    design_data=bake.to_dict(orient="records"),
    response_column="crisp",
    factors_to_plot=["T", "F"],
    factor_labels={"T": "Time [min]", "F": "Fat [g]"},
    backend="plotly",
)
fig = go.Figure(square["plotly"])
fig.update_layout(width=520, height=440, title="Baked-food 2x2: crispiness")
fig

[8]:

model = lm("crisp ~ T + F + T:F", bake)
params = model.get_parameters(drop_intercept=False)
print(params.to_string())
print()
print("Verification - the equation is:")
print(f"  crisp = {params['Intercept']:.0f} "
      f"+ {params['T']:+.0f} * x_T "
      f"+ {params['F']:+.0f} * x_F "
      f"+ {params['T:F']:+.0f} * x_T * x_F")

Intercept    49.0
T             6.0
F             2.0
T:F          -4.0

Verification - the equation is:
  crisp = 49 + +6 * x_T + +2 * x_F + -4 * x_T * x_F

[9]:

# Interaction plot: lines with very different slopes signal interaction.
# visualize_doe(plot_type="interaction") draws both orientations from
# the same call (Plotly subplots) so you can read it from either factor's
# point of view.

interaction = visualize_doe(
    plot_type="interaction",
    design_data=bake.to_dict(orient="records"),
    response_column="crisp",
    factors_to_plot=["T", "F"],
    factor_labels={"T": "Time [min]", "F": "Fat [g]"},
    backend="plotly",
)
fig = go.Figure(interaction["plotly"])
fig.update_layout(width=720, height=380, title="Interaction T x F (both orientations)")
fig

[10]:

# Predictions at the two extrapolation points the worksheet asks for.
# In coded units: T-center = 90 min, half-range = 10 min;
#                 F-center = 25 g,  half-range = 5 g.

points = [
    ("110 min, 20 g fat", dict(T=2, F=-1)),
    ("110 min, 15 g fat", dict(T=2, F=-2)),
]
for label, coded in points:
    y_hat = float(predict(model, **coded).iloc[0])
    print(f"{label:25s} -> crispness {y_hat:.1f}")

110 min, 20 g fat         -> crispness 67.0
110 min, 15 g fat         -> crispness 73.0

Guidance

An interaction plot in two orientations is the most reliable way to read a negative interaction. Plot mean crispiness against T for the two levels of F: at low fat the line is steep and rises fast; at high fat the line is much flatter. Same plot with the axes swapped: at short time the line rises with fat; at long time it almost flattens. Two lines with very different slopes - or crossing - mean an interaction is present.

Solution

Prediction equation. crisp = 49 + 6 * x_T + 2 * x_F - 4 * x_T * x_F, which is exactly what the code recovers. The negative T:F coefficient is the interaction: fat helps at short bake times, but penalises crispiness at long bake times.

Advice to the manager. The factor with the biggest pull is time; the interaction means more fat hurts once you push time to the high level. The maximum-crispiness corner is T = +1, F = -1 (100 min, 20 g) at 57 - confirmed by the data.

Physical sense. At long bake times, extra fat seals moisture in and softens the crust; at short bake times the crust has not set yet, so additional fat makes the result richer and crisper. The interaction is real, and matches the physical reasoning above.

Two predictions.

  • 110 min, 20 g fat: x_T = +2, x_F = -1 gives 49 + 12 - 2 + 8 = 67.

  • 110 min, 15 g fat: x_T = +2, x_F = -2 gives 49 + 12 - 4 + 16 = 73.

Both are well above any of the four corners, but both are extrapolations - the model assumes the linear trend keeps going, which a real oven will not.

Question 14 - the code from the source worksheet#

The week-2 worksheet links to an R snippet at yint.org/2w1 that fits the baked-food interaction model. The same job in process_improve is the block above:

from process_improve.experiments import c, gather, lm

T = c(-1, +1, -1, +1)
F = c(-1, -1, +1, +1)
y = c(37, 57, 49, 53)
bake = gather(T=T, F=F, y=y)
model = lm("y ~ T + F + T:F", bake)
print(model.get_parameters(drop_intercept=False))

No R required. The library exposes lm() with the same formula syntax (patsy under the hood) and a Model.get_parameters() helper that returns the coefficients as a pandas Series.

Wrap-up#

You have now moved from “four numbers and a manual calculation” to a coded linear model that predicts anywhere in the design space, and beyond it with appropriate caution. Two notable habits:

  • Coding factors to [-1, +1] makes coefficients directly comparable and turns the prediction equation into something you can read off in your head.

  • An interaction term costs nothing extra to fit in a 2x2 and tells you when “more is better” stops being true.

Next: Module 3 scales the same toolkit to three factors (the cube plot is back) and replicates the worksheet from week 3, where a fifth experiment at the center is added to gain a degree of freedom and start estimating curvature.