Fit and check multiple regression and GLMs, design randomised, blocked and factorial experiments, and turn analyses into clear, reproducible reports for the researchers you work with.

Regression & Experimental Design

Write the question in one sentence, then name the response (what you measure, in what units) and the treatments (the factor and its 3 or 4 levels).

Name the experimental unit: the thing each treatment is applied to independently. Check you are not counting repeated measurements of one unit as replicates.

Identify one nuisance factor you can block on (day, person, batch of materials) and decide how many blocks you can afford. Every treatment must appear in every block.

Use a random number generator to randomise the run order separately within each block (for example `random.seed(...)` then `random.shuffle(...)` in Python, or a dice roll). Write down the seed or method.

Produce the run sheet: a table with one row per block and the randomised treatment order across the columns, plus a blank column for the response.

Write the analysis model you would use, for example `response ~ C(treatment) + C(block)`, and the error degrees of freedom: (treatments - 1) x (blocks - 1).

List two things that could go wrong during the run (a dropped measurement, a changed condition) and how you would record them.

## Why this matters

Module 1 analysed data that already existed. This module is about collecting data so that the analysis is easy and the conclusion is trustworthy. A well-designed experiment of 16 runs often tells you more than 1,000 rows of historical process data, because you chose the settings and controlled how they were assigned. R. A. Fisher put it bluntly in 1938: calling in a statistician after the experiment is done is often asking for a post-mortem. The three principles below, **randomisation, replication and blocking**, are what make an experiment defensible to a reviewer, a regulator or a sceptical plant manager.

## Vocabulary

- **Experimental unit:** the thing a treatment is applied to independently (a batch, a plot, a mouse, a wafer, a customer).
- **Factor:** an input you deliberately set (temperature, supplier, dose). Its settings are **levels**.
- **Treatment:** one combination of factor levels.
- **Response:** what you measure.
- **Nuisance factor:** something that affects the response but is not of interest (day, operator, raw-material lot, instrument).

## Randomisation

Assign treatments to units, and order the runs, **at random**. This protects you against things you did not think of. If you run all the high-temperature batches in the morning and the low ones in the afternoon, any morning-versus-afternoon difference (ambient humidity, a tired operator, a warming instrument) is mixed up with temperature, and no analysis can separate them.

Randomisation does not make groups identical; it makes the leftover differences **chance** differences, which is exactly what the error term and the p-value account for. In practice:

```python
import random

random.seed(2026)                    # record the seed in your lab notebook
runs = ["low"] * 4 + ["high"] * 4
random.shuffle(runs)
print(runs)    # e.g. ['high', 'low', 'low', 'high', ...]
```

Keep the seed and the resulting run sheet; it is part of the method.

## Replication

**Replication** means applying the same treatment to more than one independent experimental unit. It gives you an estimate of pure error, without which you cannot tell a real effect from noise, and more replicates give more precision (the standard error of a mean falls with `1/sqrt(n)`).

The trap is **pseudo-replication**: measuring the same unit several times and treating the measurements as independent. Three readings from one batch tell you about measurement error, not batch-to-batch variation. If treatments are applied to batches, the batch is the unit, and `n` is the number of batches.

## Blocking

A **block** is a group of units that are more alike than the units overall: runs on the same day, specimens cut from the same sheet, patients at the same site. **Blocking** means running every treatment within each block and comparing treatments inside blocks. Block-to-block differences then cancel out of the treatment comparisons and are removed from the error term, which makes the experiment more sensitive.

The rule of thumb, often attributed to George Box: **"Block what you can, randomise what you cannot."**

## Worked example: a randomised complete block design

A lab wants to compare four buffer recipes (A, B, C, D) on enzyme activity. Only four assays fit on one plate reader run per day, and day-to-day drift in the reader is known to be large. Design:

- **Blocks:** 4 days. **Treatments:** 4 buffers. Each buffer appears once per day.
- **Randomisation:** within each day, the order (and plate position) of the four buffers is shuffled separately.

```text
           run 1   run 2   run 3   run 4
Monday       C       A       D       B
Tuesday      B       D       A       C
Wednesday    A       C       B       D
Thursday     D       B       C       A
```

The analysis model is `activity ~ C(buffer) + C(day)`. The day term soaks up the drift, and the buffer comparison has `(4 − 1) × (4 − 1) = 9` error degrees of freedom. If you had ignored the days and randomised all 16 runs freely, the day-to-day variation would sit inside the error term and could hide a real difference between buffers.

Blocks are usually treated as fixed effects in a small design like this, or as random effects when there are many of them (module 3).

## Choosing the number of runs

Size the experiment before you run it. You need three things: the **smallest difference worth detecting** (from the scientist, not the statistician), an estimate of the **noise standard deviation** (from past data or a pilot), and the desired **power** (usually 80–90%). The power ideas from Probability & Statistical Inference apply directly; `statsmodels.stats.power` has calculators for t-tests and ANOVA. If the answer is "more runs than we can afford", it is better to know before spending the budget.

## Other designs you will meet

- **Completely randomised design (CRD):** no blocks; fine when units are homogeneous.
- **Latin square:** blocks in two directions at once (day *and* operator), each treatment once per row and column.
- **Split-plot:** some factors are hard to change (oven temperature), so they are randomised to large "whole plots" and easy factors to "subplots" inside them. Split-plots need a mixed-model analysis; analysing them as if fully randomised overstates precision for the hard-to-change factor.
- **Crossover:** each subject gets several treatments in random order, with washout periods.

## Common mistakes

- **Running treatments in a convenient order** (all of A, then all of B). It confounds treatment with time.
- **Calling repeated measurements replicates.**
- **Blocking but analysing as if unblocked**, throwing away the precision you paid for.
- **Changing the plan mid-experiment** without recording it. Deviations are normal; undocumented ones are fatal.

## Quick recap

- Randomise to protect against unknown nuisance variables.
- Replicate independent units to estimate error; beware pseudo-replication.
- Block known nuisance variables and compare treatments within blocks.
- Size the experiment with power before you start.

## Try it

The **Do** tab walks you through designing a small blocked experiment on paper, from choosing the unit and blocks to producing a randomised run sheet.


Apply the three principles of experimental design and lay out a randomised complete block experiment.

Randomisation, replication, blocking

Replication, Randomization & Blocking

Randomized Complete Block Design

## Why this matters

You have run your randomised experiment: three catalysts, five batches each. Are the catalysts really different, or are the differences in average yield just batch-to-batch noise? Running three separate t-tests (A vs B, A vs C, B vs C) inflates the chance of a false positive. **Analysis of variance (ANOVA)** answers the overall question with one test, and its table is the standard way results of designed experiments are reported in chemistry, biology, materials and engineering papers.

## The core idea

ANOVA compares two estimates of variability:

- **Between groups:** how far the group means are from the grand mean.
- **Within groups:** how much individual results scatter around their own group mean (pure noise, if the experiment was done properly).

If all treatment means are truly equal, both estimate the same noise variance and their ratio is near 1. If some treatments differ, the between-group estimate is inflated and the ratio grows. That ratio is the **F statistic**.

## The decomposition

With `k` groups, `n_i` observations in group `i`, `N` observations in total, group means `ȳ_i` and grand mean `ȳ`:

```text
SS_between = sum over groups  n_i * (ȳ_i - ȳ)^2          df = k - 1
SS_within  = sum over all obs (y_ij - ȳ_i)^2              df = N - k
SS_total   = SS_between + SS_within                       df = N - 1

MS = SS / df          F = MS_between / MS_within
```

Under the null hypothesis (all means equal) and the usual assumptions, `F` follows an F distribution with `(k − 1, N − k)` degrees of freedom.

## Worked example: three catalysts

```text
Catalyst A: 82 85 88 84 86   mean 85
Catalyst B: 90 92 89 93 91   mean 91
Catalyst C: 84 83 86 85 82   mean 84     grand mean 86.667
```

Between groups: `5 × [(85 − 86.667)² + (91 − 86.667)² + (84 − 86.667)²] = 5 × 28.667 = 143.33`.

Within groups: A deviations −3, 0, 3, −1, 1 give 20; B gives 10; C gives 10; total 40.

```text
Source     SS       df    MS       F
Catalyst   143.33    2    71.67   21.5
Error       40.00   12     3.33
Total      183.33   14
```

The 5% critical value of F(2, 12) is 3.89, and 21.5 is far above it (p ≈ 0.0001). We reject the hypothesis that all three catalysts give the same mean yield. A useful effect size is `η² = SS_between / SS_total = 0.78`: catalyst explains 78% of the variation in these runs.

## The same thing in Python

```python
from scipy import stats
import statsmodels.formula.api as smf
from statsmodels.stats.anova import anova_lm

a = [82, 85, 88, 84, 86]; b = [90, 92, 89, 93, 91]; c = [84, 83, 86, 85, 82]
print(stats.f_oneway(a, b, c))          # F = 21.5, p ≈ 1e-4

# The regression view: ANOVA is a linear model with a categorical predictor
model = smf.ols("yield_pct ~ C(catalyst)", data=df).fit()
print(anova_lm(model, typ=2))
```

The second form matters. ANOVA *is* multiple regression with dummy variables, so everything from module 1 applies: residual diagnostics, adding a block term (`+ C(day)`), and adding covariates (ANCOVA).

## After a significant F: which groups differ?

The F test says "not all equal". To say which pairs differ, use a **multiple-comparison procedure** that controls the family-wise error rate:

```python
from statsmodels.stats.multicomp import pairwise_tukeyhsd
print(pairwise_tukeyhsd(df["yield_pct"], df["catalyst"], alpha=0.05))
```

**Tukey's HSD** compares all pairs. Here it finds B higher than both A and C, while A and C are not clearly different. If you only care about comparing each treatment with a control, **Dunnett's test** is more powerful. If you planned a few specific contrasts in advance, test those directly.

## Checking the assumptions

ANOVA assumes independent observations, roughly normal residuals and **equal variances** across groups. Check:

- Residuals vs fitted and a Q-Q plot, as in the diagnostics lesson.
- Group standard deviations: a ratio of largest to smallest SD below about 2 is usually fine, especially with equal group sizes.
- If variances clearly differ, use **Welch's ANOVA** (`pingouin.welch_anova` or `statsmodels.stats.oneway.anova_oneway(..., use_var="unequal")`).
- For badly skewed data with small groups, the rank-based **Kruskal–Wallis** test is an alternative, though it tests a slightly different hypothesis.

Balanced designs (equal `n` per group) make ANOVA much more robust to unequal variances, one more reason to plan the experiment rather than take what comes.

## Common mistakes

- **Running all pairwise t-tests without correction.** With 5 groups there are 10 pairs, and the chance of at least one false positive at 5% each is about 40%.
- **Concluding "no difference" from a non-significant F.** It may mean the experiment was too small. Report confidence intervals for the differences.
- **Ignoring the block** you designed in. If the experiment was blocked, the model must include the block.
- **Treating repeated measurements as replicates**, which shrinks the error term and inflates F.

## Quick recap

- ANOVA splits total variation into between-group and within-group parts.
- `F = MS_between / MS_within`; compare with F(k − 1, N − k).
- A significant F needs a follow-up (Tukey, Dunnett, planned contrasts).
- ANOVA is a linear model, so diagnose residuals and include design terms such as blocks.

## Try it

In the **Do** tab you will compute the ANOVA table for the catalyst data from scratch with the standard library and confirm the F statistic of 21.5.


Compute and interpret a one-way ANOVA table and follow it with the right multiple comparisons.

ANOVA

ANOVA: Crash Course Statistics #33

One-Way ANOVA [Analysis of Variance] simply explained

## Why this matters

The intuitive way to study several factors is **one factor at a time (OFAT)**: hold everything fixed, vary temperature, pick the best; then vary pressure, pick the best; and so on. It is slow, and worse, it can miss the best settings completely when factors interact. **Factorial designs** vary all factors together in a planned pattern. With the same number of runs they estimate every main effect more precisely than OFAT and they reveal interactions. The **2^k factorial**, with every factor at two levels, is the backbone of industrial experimentation, process development and screening studies.

## The 2^k design

Each of `k` factors gets a **low** and a **high** level, coded −1 and +1. Run every combination: `2^k` runs. For three factors, temperature (A), pressure (B) and catalyst type (C), that is 8 runs. In **standard (Yates) order**, A alternates fastest:

```text
run    A    B    C    yield
 1    -1   -1   -1     60
 2    +1   -1   -1     72
 3    -1   +1   -1     54
 4    +1   +1   -1     68
 5    -1   -1   +1     52
 6    +1   -1   +1     83
 7    -1   +1   +1     45
 8    +1   +1   +1     80
```

Standard order is for writing the design down. You **run** it in a randomised order, as in lesson 4.

## Estimating effects

The **main effect** of a factor is the average response at its high level minus the average at its low level:

```text
A effect = mean(runs 2,4,6,8) - mean(runs 1,3,5,7)
         = (72 + 68 + 83 + 80)/4 - (60 + 54 + 52 + 45)/4 = 75.75 - 52.75 = 23.0
```

Every run is used for every effect. That is where the efficiency comes from.

For an **interaction**, multiply the coded columns. The AC column is `A × C`, and the AC effect is the mean response where AC = +1 minus the mean where AC = −1:

```text
AC signs:  +1 -1 +1 -1 -1 +1 -1 +1
AC effect = (60 + 54 + 83 + 80)/4 - (72 + 68 + 52 + 45)/4 = 69.25 - 59.25 = 10.0
```

The full set for this example:

```text
A = 23.0   B = -5.0   C = 1.5   AB = 1.5   AC = 10.0   BC = 0.0   ABC = 0.5
```

## Interpreting an interaction

AC = 10 means the effect of temperature **depends on** the catalyst. Split by catalyst:

```text
catalyst C = -1:  temperature effect = mean(72, 68) - mean(60, 54) = 70 - 57 = 13
catalyst C = +1:  temperature effect = mean(83, 80) - mean(52, 45) = 81.5 - 48.5 = 33
```

The AC effect is half the difference: (33 − 13)/2 = 10. In an **interaction plot** (temperature on the x-axis, one line per catalyst) the lines are not parallel. The main effect of C is small (1.5), but saying "catalyst does not matter" would be wrong: catalyst +1 is worse at low temperature and much better at high temperature. When an interaction is large, interpret the factors together, not separately.

OFAT can easily miss this. If the team had first compared the two catalysts at the low temperature and low pressure, catalyst +1 would have looked worse (52 against 60) and been dropped, and the best combination in the design (high temperature with catalyst +1, 83) would never have been run.

## Which effects are real?

With a single replicate there is no pure-error estimate, so you cannot run an ordinary ANOVA. Options:

1. **Replicate** the design (16 runs for two replicates) and use ANOVA or regression with the usual error term.
2. **Add centre points** (runs at the middle of every factor): they give a pure-error estimate and a test for curvature.
3. Use a **normal probability (Daniel) plot** or **half-normal plot** of the effects. Most effects in a real process are near zero (effect sparsity), so they fall on a straight line; real effects stand out. **Lenth's method** turns that into a formal test.
4. Assume high-order interactions (ABC) are negligible and pool them into error.

In regression form, a 2^k analysis is simply `y ~ A * B * C` with coded −1/+1 predictors. Each regression coefficient is **half** the effect (because the coded range from −1 to +1 is 2 units):

```python
import statsmodels.formula.api as smf
model = smf.ols("yld ~ A * B * C", data=df).fit()   # A, B, C coded -1/+1
print(2 * model.params)                              # effects
```

With 8 runs and 8 parameters this model fits perfectly, which is why you need replicates, centre points or a pooled error term for inference.

## Fractional factorials

With 7 factors a full factorial needs 128 runs. A **fractional factorial** runs a planned subset, for example `2^(7−4) = 8` runs, at the cost of **aliasing**: some effects cannot be separated from others. The **resolution** tells you how bad the aliasing is. In a resolution III design, main effects are aliased with two-factor interactions; in resolution IV, main effects are clear of two-factor interactions; in resolution V, two-factor interactions are clear of each other. Screening experiments often start with a resolution III or IV fraction, then "fold over" or augment the design to untangle what matters.

## Common mistakes

- **Running in standard order** instead of randomising.
- **Choosing levels too close together**, so effects are lost in noise. Make the range as wide as is safe and practical.
- **Interpreting main effects when a large interaction is present.**
- **Ignoring curvature.** A two-level design only sees straight lines. Centre points tell you whether the optimum may lie inside the region, which leads to the response-surface methods in module 3.

## Quick recap

- A 2^k factorial runs every combination of k two-level factors.
- Effect = mean at +1 minus mean at −1, using contrast columns; interactions use products of columns.
- Interactions mean "it depends": interpret factors together.
- Without replication, judge effects with a normal plot, centre points or pooled interactions.

## Try it

In the **Do** tab you will build contrast columns for a 2³ design and compute all seven effects for the yield data above.


Estimate main effects and interactions from a 2^k factorial experiment and interpret them.

Factorial designs

The 2^2 Factorial Design, Part 1

Full Factorial Design (DoE - Design of Experiments) Simply explained

Plan experiments with randomisation, replication and blocking, analyse them with one-way ANOVA and follow-up comparisons, and estimate main effects and interactions in 2^k factorial designs.