Design of Experiments Worksheet Template (Statistics, R and Python)

Design of Experiments Worksheet: [Your experiment or course]

A design of experiments (DoE) worksheet for students and analysts who must plan an experiment before running it: choose a design, write the randomisation plan, fix the analysis model, and check the plan against the classic principles. It speeds up learning by putting each design next to its model, its code and a worked example, then fading the help until you plan one yourself. The worked examples use fictional data; replace them with your own. Source basis: Fisher (1935), Box, Hunter and Hunter (2005) and the NIST/SEMATECH e-Handbook of Statistical Methods, chapter 5 (public domain), written in original words.

Document control

This table identifies the experiment plan and the software used to create the randomisation and the analysis.

When to use this worksheet

Use this worksheet when you control which experimental unit receives which treatment. That control is what separates an experiment from an observational study, and it is the reason a well-designed experiment supports a causal conclusion: if the treatment was assigned by a random mechanism, then nothing about the units can explain a systematic difference between the groups except the treatment itself (plus chance, which the analysis quantifies).

Alternatives worth knowing: split-plot designs (when one factor is hard to change), response surface designs such as central composite designs (when you want to find an optimum after screening), and crossover designs (when each person receives every treatment in sequence).

Key terms and notation

An experimental unit is the smallest thing to which a treatment is independently applied: a batch of dough, a loaf, a tray. A factor is a variable you set (proof time), its levels are the values you choose (40 and 60 minutes), and a treatment is one combination of levels. The response is what you measure (loaf height in mm). A block is a group of units expected to be alike (one day of baking), and a run is one execution of a treatment on a unit.

The four principles

Fisher set out the logic of experiments in The Design of Experiments (1935). Four principles carry almost all of it.

1. Randomisation. Assign treatments to units, and order the runs, by a random mechanism with a recorded seed. Randomisation protects against bias from things you did not think of (the oven warming up during the morning, a tired baker in the afternoon), and it justifies the probability statements in the analysis.
2. Replication. Apply each treatment to more than one independent unit. Replication lets you estimate the error variance $\sigma^2$, and more replicates make estimates more precise. Measuring the same loaf three times is repeated measurement, not replication: it tells you about the ruler, not about the process.
3. Blocking. When a nuisance source is known (day, flour batch, operator), group the runs so that every treatment appears within each block, then remove the block-to-block variation in the analysis. A rule of thumb popularised by Box, Hunter and Hunter: block the nuisance sources you can identify, and randomise to protect against the rest.
4. Control. Hold constant everything that is not under study (same mixer, same scale, same measurement method), and include a control or reference treatment where a comparison needs one. Local control also means measuring the response the same way, ideally by someone who does not know which treatment a unit received (blinding).

Which design?

The flow below turns the four principles into a design choice. Start with how many factors you want to study, then ask whether a nuisance source must be blocked.

Completely randomised design

In a completely randomised design (CRD), every unit has the same chance of receiving any treatment. It is the simplest design and the most flexible (unequal group sizes are fine), but it puts all unit-to-unit variation into the error term. Use it when the units are genuinely alike.

The model is

Example: Brightwell Bakery compares three yeast strains (Y1, Y2, Y3) on loaf volume in mL, with 12 dough batches assigned at random, four per strain. The analysis is a one-way ANOVA (see academic/stats-one-way-anova).

crd <- data.frame(yeast = factor(rep(c("Y1", "Y2", "Y3"), each = 4)),
                  volume = c(1510, 1495, 1530, 1505, 1560, 1545, 1575, 1550,
                             1520, 1540, 1515, 1535))
summary(aov(volume ~ yeast, data = crd))

            Df Sum Sq Mean Sq F value  Pr(>F)
yeast        2   4617  2308.3   12.98 0.00223 **
Residuals    9   1600   177.8

The strains differ, F(2, 9) = 12.98, p = .002. The Python version with anova_lm(smf.ols("volume ~ C(yeast)", data=crd).fit()) gives the same table. A follow-up such as Tukey HSD tells you which strains differ.

Randomised complete block design

In a randomised complete block design (RCBD), units are grouped into blocks that are alike, every treatment appears once in every block, and the order within each block is randomised separately. The block effect is then removed from the error, so the comparison of treatments is made within blocks, where conditions are similar.

Blocks are not tested as a hypothesis of interest: they are in the model because the design put them there. The usual analysis has no treatment-by-block interaction term, which is an assumption (additivity) that you can check with a plot of the response against block for each treatment: the lines should be roughly parallel. Worked example B uses an RCBD to compare three flours over four baking days.

Latin square design

A Latin square blocks two nuisance sources at once. With $t$ treatments, arrange a $t \times t$ grid whose rows are levels of the first nuisance source and whose columns are levels of the second; each treatment letter appears exactly once in each row and once in each column. It is efficient (only $t^2$ runs) but assumes no interactions among rows, columns and treatments, and it gives few error degrees of freedom when $t$ is small.

Example: four recipes (A, B, C, D) are baked in four ovens (rows) across four shifts (columns), and a trained taster scores crumb softness from 0 to 100. To randomise, start from any standard square, then shuffle the rows, shuffle the columns and assign recipes to letters at random.

ls_dat <- data.frame(
  oven = factor(rep(1:4, each = 4)), shift = factor(rep(1:4, 4)),
  recipe = factor(c("A","B","C","D", "B","C","D","A", "C","D","A","B", "D","A","B","C")),
  score = c(70, 71.5, 76, 69, 74.5, 75, 71.5, 73, 74, 65.5, 70, 72.5, 69, 68.5, 75.5, 77))
summary(aov(score ~ oven + shift + recipe, data = ls_dat))

            Df Sum Sq Mean Sq F value   Pr(>F)
oven         3  19.55    6.52   40.35 0.000226 ***
shift        3  23.42    7.81   48.35 0.000135 ***
recipe       3 110.80   36.93  228.74 1.42e-06 ***
Residuals    6   0.97    0.16

Recipe means are A 70.375, B 73.500, C 75.500, D 68.750. Only 6 error df remain out of 15, which is why Latin squares smaller than 4 by 4 are usually replicated.

Full factorial 2^k designs

A full factorial runs every combination of the factor levels. With $k$ factors at two levels each, that is $2^k$ runs per replicate. Factorials beat "one factor at a time" experiments on two counts: each run informs every main effect (so the same number of runs gives more precise estimates), and they reveal interactions, where the effect of one factor depends on the level of another.

In plain words, the main effect of A is "the average response at high A minus the average response at low A". Write the design in standard order (A alternates fastest, then B, then C), give each run a sign in the A column, and the effect is a signed sum (a contrast) divided by half the number of runs:

The sign column of an interaction such as AB is the product of the A and B columns, so the same formula estimates interactions. With an error standard deviation $s$ estimated from replicates or centre points,

That factor of one half trips up many learners: software reports coefficients, textbooks often report effects. Worked example A estimates all seven effects of a $2^3$ design by contrasts.

Fractional factorial designs for screening

With five or more factors a full factorial gets expensive ($2^6 = 64$ runs). A fractional factorial $2^{k-p}$ runs a carefully chosen $1/2^p$ of the combinations. The price is aliasing: some effects share the same sign column, so the data can only estimate their sum. You choose the fraction through a generator, for example $D = ABC$ for a $2^{4-1}$ design. Multiplying both sides by $D$ (any column times itself is a column of $+1$, written $I$) gives the defining relation

To find what any effect is aliased with, multiply it by the defining relation and cancel squared letters: $A \cdot ABCD = BCD$, $AB \cdot ABCD = CD$. So in this design A is aliased with the three-factor interaction BCD (usually negligible), but AB is aliased with CD (both may matter).

The resolution is the length of the shortest word in the defining relation and summarises the damage:

The $2^{4-1}$ design with $I = ABCD$ is resolution IV. Screening rests on effect sparsity (few factors matter), hierarchy (main effects tend to be larger than interactions) and heredity (an interaction is more likely when its parent main effects are active). When an alias set is ambiguous, add a second fraction chosen to break it. For a resolution III design, the full fold-over (the same runs with every sign reversed) frees all main effects from two-factor interactions. For a resolution IV design the full fold-over does not help (every word in the defining relation has even length, so it is unchanged); instead reverse the sign of one factor, which here is the same as running the complementary half ($I = -ABCD$) and together gives the full $2^4$. Box and Hunter (1961) set out the theory of these $2^{k-p}$ designs.

Centre points and curvature

A two-level design assumes the response changes linearly between the low and high levels. Running a few replicates at the centre (every coded factor at 0, for example 50 minutes and 210 °C) checks that assumption without adding factor levels to every run. Centre points give a pure-error estimate $s$ even when the factorial part is unreplicated, and they test curvature by comparing the centre mean with the factorial mean:

In example A, the centre mean is 97.25 mm and the factorial mean is 95.5 mm, so the curvature estimate is 1.75 mm, t(3) = 2.99, p = .058. That is weak evidence of curvature. A significant result would say that the optimum may lie inside the region and that a response surface design is the next step. Spread centre points through the run order (at the start, middle and end), so they also detect drift.

Randomisation plan generated in R and Python

The plan below randomises the 12 runs of example A (8 factorial runs plus 4 centre points) and, separately, the flour order within each day for example B. Record the seed, the software and the version in Document control so anyone can regenerate the same order.

set.seed(2026)
design <- expand.grid(A = c(-1, 1), B = c(-1, 1), C = c(-1, 1))
design <- rbind(design, data.frame(A = 0, B = 0, C = 0)[rep(1, 4), ])
design$std_order <- 1:12
design$run_order <- sample(12)
print(design[order(design$run_order), ], row.names = FALSE)

set.seed(2026)
for (d in c("Mon", "Tue", "Wed", "Thu")) cat(d, ":", sample(c("F1", "F2", "F3")), "\n")

  A  B  C std_order run_order
  1 -1 -1         2         1
  1  1  1         8         2
  0  0  0        10         3
  1  1 -1         4         4
  1 -1  1         6         5
 -1  1 -1         3         6
  0  0  0         9         7
 -1 -1  1         5         8
 -1 -1 -1         1         9
  0  0  0        12        10
  0  0  0        11        11
 -1  1  1         7        12
Mon : F1 F3 F2
Tue : F2 F1 F3
Wed : F3 F2 F1
Thu : F1 F3 F2

import itertools
import numpy as np
import pandas as pd

rng = np.random.default_rng(2026)
levels = [(a, b, c) for c, b, a in itertools.product([-1, 1], repeat=3)]
design = pd.DataFrame(levels + [(0, 0, 0)] * 4, columns=["A", "B", "C"])
design["std_order"] = range(1, 13)
design["run_order"] = rng.permutation(12) + 1
print(design.sort_values("run_order").to_string(index=False))

 A  B  C  std_order  run_order
 1 -1  1          6          1
 1 -1 -1          2          2
-1  1  1          7          3
 1  1  1          8          4
 0  0  0         11          5
 0  0  0          9          6
-1 -1 -1          1          7
 0  0  0         10          8
 1  1 -1          4          9
-1  1 -1          3         10
 0  0  0         12         11
-1 -1  1          5         12

R and Python give different orders for the same seed. That is by design: R uses the Mersenne Twister with its own sampling method, and numpy's default_rng uses the PCG64 generator. Both orders are valid randomisations. What matters is that you pick one tool, record it with the seed and version, and follow that printed order. Do not rerun with new seeds until you like the order: that quietly undoes the randomisation.

Analysis model for each design

The design fixes the analysis. The table pairs each design with its model and the code that fits it.

For balanced designs, R's aov and statsmodels' anova_lm (type I sums of squares by default) agree. If runs are lost and the design becomes unbalanced, the order of terms matters for type I sums of squares; see academic/stats-factorial-anova for type II and III.

Worked example A: 2^3 factorial at Brightwell Bakery (full solution)

Question. Which of three settings raise loaf height: A proof time (40 vs 60 min), B oven temperature (200 vs 220 °C), C sugar (2% vs 4%)? The bakery ran one replicate of the $2^3$ design plus four centre points, in the R run order above. Responses in standard order:

Step 1, effect of A by hand. Contrast = $-92 + 96 - 93 + 103 - 91 + 95 - 92 + 102 = 28$, so $\text{Effect}_A = 28 / 4 = 7$ mm. Equivalently, the mean at high A is 99 mm and at low A is 92 mm.

Step 2, interaction AB by hand. The AB sign column is $(+,-,-,+,+,-,-,+)$. Contrast $= 92 - 96 - 93 + 103 + 91 - 95 - 92 + 102 = 12$, so $\text{Effect}_{AB} = 3$ mm.

Step 3, all effects by software.

fac <- expand.grid(A = c(-1, 1), B = c(-1, 1), C = c(-1, 1))
fac$height <- c(92, 96, 93, 103, 91, 95, 92, 102)
X <- model.matrix(~ A * B * C, data = fac)
round(colSums(X * fac$height) / 4, 2)[-1]   # contrast / (N/2)

    A     B     C   A:B   A:C   B:C A:B:C
    7     4    -1     3     0     0     0

import itertools
import numpy as np
import pandas as pd

levels = [(a, b, c) for c, b, a in itertools.product([-1, 1], repeat=3)]
fac = pd.DataFrame(levels, columns=["A", "B", "C"])
fac["height"] = [92, 96, 93, 103, 91, 95, 92, 102]
for name in ["A", "B", "C", "A:B", "A:C", "B:C", "A:B:C"]:
    col = np.prod([fac[f] for f in name.split(":")], axis=0)
    print(name, (col * fac["height"]).sum() / 4)

A 7.0
B 4.0
C -1.0
A:B 3.0
A:C 0.0
B:C 0.0
A:B:C 0.0

Step 4, error and significance. The centre points give $s = 0.957$ mm with 3 df, so $SE(\text{Effect}) = 2(0.957)/\sqrt{8} = 0.677$ mm. Then A: t(3) = 10.34, p = .002; B: t(3) = 5.91, p = .010; AB: t(3) = 4.43, p = .021; C: t(3) = -1.48, p = .236. If you instead fit lm(height ~ A + B + C + A:B) to all 12 runs, the coefficients are half the effects (A 3.500, B 2.000, C -0.500, A:B 1.500) with SE 0.442, because the residual now pools pure error with the dropped interactions and curvature. Both routes lead to the same conclusion.

Step 5, picture. The interaction plot shows the cell means for A and B (averaged over C). Non-parallel lines are the AB interaction: longer proofing helps much more in the hotter oven.

Step 6, interpretation and report. Proof time and oven temperature both raise height, and they interact; sugar has no detectable effect in this range. The best setting tested is 60 minutes at 220 °C (cell mean 102.5 mm). Curvature is borderline (p = .058), so a follow-up around the high corner is sensible.

Report sentence: "In an unreplicated $2^3$ factorial with four centre points, longer proofing increased loaf height by 7.0 mm (SE 0.68, t(3) = 10.34, p = .002) and the higher oven temperature by 4.0 mm (t(3) = 5.91, p = .010), with a proof by temperature interaction of 3.0 mm (t(3) = 4.43, p = .021); sugar had no detectable effect (-1.0 mm, p = .236)."

Worked example B: blocking the flour trial (half solved)

Question. Do three flours (F1, F2, F3) differ in dough rise after 60 minutes (mm)? Humidity changes from day to day, so the bakery used an RCBD with four days as blocks and the randomised order printed in the plan above.

rcbd <- data.frame(day = factor(rep(c("Mon", "Tue", "Wed", "Thu"), each = 3)),
                   flour = factor(rep(c("F1", "F2", "F3"), 4)),
                   rise = c(42, 45, 47, 38, 42, 44, 44, 46, 50, 40, 43, 45))
summary(aov(rise ~ day + flour, data = rcbd))   # blocked: the correct model
summary(aov(rise ~ flour, data = rcbd))         # ignores the blocks

            Df Sum Sq Mean Sq F value   Pr(>F)
day          3  49.00  16.333      49  0.00013 ***
flour        2  60.67  30.333      91 3.25e-05 ***
Residuals    6   2.00   0.333

            Df Sum Sq Mean Sq F value Pr(>F)
flour        2  60.67  30.333   5.353 0.0294 *
Residuals    9  51.00   5.667

1. Flour means: F1 = 41.0, F2 = 44.0, F3 = [ ] mm.
2. The model is $y_{ij} = \mu + \tau_i + \beta_j + \varepsilon_{ij}$, where $\beta_j$ is the effect of [ ].
3. Blocking moved [ ] of the 51.00 units of sum of squares out of the error, so the error mean square falls from 5.667 to [ ].
4. The F statistic for flour rises from 5.353 to [ ], because $F = MS_{\text{flour}} / MS_{\text{error}}$ = 30.333 / [ ].
5. Report: "Flour affected dough rise, F(2, [ ]) = [ ], p < .001, after blocking on baking day."
6. Why is it wrong to report the unblocked analysis even though it is also significant? [ ]

Exercise C: screening with a half fraction

The bakery adds a fourth factor, D dough hydration (65% vs 70%), and can only afford 8 runs. It uses the $2^{4-1}$ design with generator $D = ABC$. Responses (loaf height, mm) in standard order of A, B, C:

1. Write the defining relation and the resolution.
2. List the alias of each main effect and the three alias pairs of two-factor interactions.
3. Compute the effects of A, B, C and D, and the contrast estimating AB + CD.
4. Which factors would you carry forward, and what follow-up would separate AB from CD?

My experiment plan

Fill this plan for your own experiment. Each line answers one decision in the order you should make it.

Common mistakes and misconceptions

• Pseudo-replication. Cutting one loaf into four slices and treating them as four replicates. The unit that received the treatment is the loaf, so $n = 1$. The error is underestimated and p-values are far too small.
• One factor at a time. Changing proof time with temperature fixed, then temperature with proof fixed. It needs more runs for the same precision and cannot see the AB interaction in example A.
• Blocking in the design but not in the analysis. If the runs were blocked by day, the model must contain day. Dropping it inflates the error (example B) and misstates the df.
• Convenience order instead of a random order. Running all low-temperature runs first because the oven is cold confounds temperature with warm-up drift. If a factor is truly hard to change, plan a split-plot design openly.
• Confusing effects with coefficients. On the $\pm 1$ scale, a regression coefficient is half the effect.
• Reading a fractional factorial as if it were full. In a resolution IV design, "AB is significant" really means "AB + CD is significant".
• Treating a non-significant interaction as proof of additivity. With few runs, power to detect interactions is low. Absence of evidence is not evidence of absence.
• Re-drawing the randomisation until it looks nice. That is no longer random. Fix the seed first.

Self-test questions

1. A teacher at Lakemoor University gives two teaching methods to two classes and measures each of the 60 students. How many experimental units are there per method, and why?
2. Name the four principles of experimental design and give one bakery example of each.
3. In a $2^3$ design with responses 10, 14, 11, 17, 10, 15, 12, 18 in standard order, what is the effect of A?
4. In a $2^{5-2}$ design with generators $D = AB$ and $E = AC$, what is the resolution?
5. When would you choose a Latin square over an RCBD?
6. Your R and Python scripts give different run orders for seed 2026. Is one of them wrong?

My notes and next steps

• Confidence (1 to 5) with: principles [ ], CRD/RCBD [ ], Latin square [ ], factorial effects [ ], aliasing [ ].
• Open questions: [ ]
• Mistakes to log in academic/stats-error-log: [ ]
• Next review date (spacing): YYYY-MM-DD. Next template: academic/stats-factorial-anova or academic/stats-analysis-plan.

Sources and further reading

Structure aligned with the NIST/SEMATECH e-Handbook of Statistical Methods, chapter 5 (Process Improvement), and the design principles in Fisher (1935) and Box, Hunter and Hunter (2005), written in original words with fictional data.

• Fisher, R. A. (1935). The Design of Experiments. Oliver and Boyd.
• Box, G. E. P., Hunter, J. S., and Hunter, W. G. (2005). Statistics for Experimenters: Design, Innovation, and Discovery (2nd ed.). Wiley.
• Box, G. E. P., and Hunter, J. S. (1961). The 2^(k-p) fractional factorial designs, Part I. Technometrics, 3(3), 311–351; Part II, Technometrics, 3(4), 449–458.
• Yates, F. (1937). The Design and Analysis of Factorial Experiments (Technical Communication No. 35). Imperial Bureau of Soil Science.
• NIST/SEMATECH e-Handbook of Statistical Methods, section 5 Process Improvement, https://www.itl.nist.gov/div898/handbook/ (accessed 2026-10-25). Public domain.
• Montgomery, D. C. (2020). Design and Analysis of Experiments (10th ed.). Wiley.
• Further reading (by name only): de Smith, M. J. (2024). Statistical Analysis Handbook, 2024 edition. StatsRef, https://www.statsref.com/, sections on experimental design.
• Dunlosky, J., Rawson, K. A., Marsh, E. J., Nathan, M. J., and Willingham, D. T. (2013). Improving students' learning with effective learning techniques. Psychological Science in the Public Interest, 14(1), 4–58.

Answer key

Example B. (1) F3 = 46.5 mm. (2) Day, the block. (3) 49.00 moves to the day term, so the error mean square falls to 2.00 / 6 = 0.333. (4) F = 30.333 / 0.333 = 91. (5) F(2, 6) = 91.0, p < .001. (6) The design restricted the randomisation to within days, so the correct error term is the within-day variation; the unblocked analysis uses the wrong error and the wrong df, even when it happens to agree in direction.

Exercise C. (1) $I = ABCD$, resolution IV. (2) A = BCD, B = ACD, C = ABD, D = ABC; AB = CD, AC = BD, AD = BC. (3) A = 9, B = 5, C = -1, D = 0, AB + CD = 2 (AC + BD = 0, AD + BC = 0), verified in R with colSums(model.matrix(~ (A + B + C + D)^2, frac) * frac$height) / 4. (4) Carry A and B forward; hydration and sugar look inactive. Because AB and CD are aliased, reversing the sign of D in a second block of 8 runs (the complementary half, $I = -ABCD$) separates them; a full fold-over of every column would not, because $ABCD$ keeps its sign. Effect heredity suggests AB is the more likely source, since A and B are active and C and D are not.

Self-test. (1) One unit per method: the class received the method, so the 60 students are measurements within a unit (pseudo-replication); you need several classes per method. (2) Randomisation (random run order), replication (several batches per recipe), blocking (day as block), control (same mixer and scale). (3) High-A runs: 14 + 17 + 15 + 18 = 64; low-A runs: 10 + 11 + 10 + 12 = 43; effect = 21 / 4 = 5.25. (4) $I = ABD = ACE = BCDE$; the shortest word has three letters, so resolution III. (5) When two crossing nuisance sources (oven and shift) must be blocked and the number of treatments equals the levels of each. (6) Neither: the generators differ by design. Use one, record it, follow it.

Revision history