3  Experimental design

These notes were automatically generated by AI and lightly reviewed by the lecturer. They may still contain mistakes and factual errors.

Let’s finish the argument about Fisher’s exact test, then turn to experimental design. Randomization gives us a way to analyse an experiment, but we also have choices about how to randomize. Those choices can make the analysis more precise, make a gradual change in practice possible, or let us stop a trial early.

3.1 Exact tests from the randomization

3.1.1 Why are the margins fixed?

Recall the \(2\times2\) table from Lecture 2 (Section 2.4.1). Let \(N_{ay}\) count units with exposure \(A=a\) and outcome \(Y=y\). The row totals \(n_0,n_1\) are fixed by complete randomization, and \(n=n_0+n_1\).

	\(Y=0\)	\(Y=1\)	Total
\(A=0\)	\(N_{00}\)	\(N_{01}\)	\(n_0\)
\(A=1\)	\(N_{10}\)	\(N_{11}\)	\(n_1\)
Total	\(N_{\cdot0}\)	\(N_{\cdot1}\)	\(n\)

In the two-binomials model, we test equality of the two success probabilities. Conditioning on the total number of successes removes the unknown common probability under the null:

\[\operatorname{\mathsf P}_{H_0} (N_{11}=n_{11}\mid N_{\cdot1}=n_{\cdot1}) =\frac{\binom{n_1}{n_{11}}\binom{n_0}{n_{\cdot1}-n_{11}}} {\binom n{n_{\cdot1}}}.\]

This is the hypergeometric distribution. Fisher’s exact test compares the observed count with it; unusually small or large values provide evidence against the null. The exponential-family argument in Lecture 2 (the optimality statement) also gives an optimality result for the appropriately calibrated randomized conditional test: it is uniformly most powerful among unbiased level-\(\alpha\) tests.

But why condition on the column total at all? In the two-binomials model, the row totals are fixed, but the total number of successes is random. George Barnard challenged Fisher on this choice: an unconditional test in the two-binomials problem can have more power in some cases. The disagreement became rather heated; the example sheet will explore the statistical argument, rather than the fight.

The randomization perspective gives a different reason for fixing the margins. In the Neyman–Rubin model, consider the sharp null

\[H_0:Y_i(0)=Y_i(1),\qquad i=1,\ldots,n.\]

This says that treatment makes no difference to any of these individuals. It is stronger than equality of two population success probabilities, and it does not refer to a superpopulation.

Under this null, consistency gives \(Y_i=Y_i(0)\) regardless of the assignment. Consequently,

\[N_{\cdot1}=\sum_{i=1}^nY_i=\sum_{i=1}^nY_i(0).\]

If we randomize the assignment again, the four cells can change, but the column totals remain fixed by the potential outcomes. This makes the fixed-margin comparison natural. The classical optimality theory gives a separate justification; it is not needed to obtain this distribution.

3.1.2 A general randomization test

We can use the same idea much more generally.

Definition 3.1.

Let \(X\) denote the covariates, \(Y(\cdot)\) the potential-outcome schedule, and \(Z\in\mathcal Z\) the assignment. Recall the two parts of the randomization assumption (Section 1.4.5):

\[Z\mathrel{\perp\!\!\!\perp}Y(\cdot)\mid X, \qquad \operatorname{\mathsf P}(Z=z\mid X)\text{ is known}.\]

Choose any real-valued test statistic \(T(X,Y(\cdot),Z)\). For now, take smaller values as evidence against the null. Let \(Z^*\) be an independent copy of \(Z\), conditional on \(X,Y(\cdot)\): a fresh assignment from the same experiment. Define the randomization \(p\)-value

\[Q=\operatorname{\mathsf P}\{T(X,Y(\cdot),Z^*)\le T(X,Y(\cdot),Z) \mid X,Y(\cdot),Z\}.\]

Inside this probability calculation, the observed assignment \(Z\) is fixed; only \(Z^*\) varies. For a discrete assignment space, we can write it out:

\[Q=\sum_{z\in\mathcal Z}\operatorname{\mathsf P}(Z=z\mid X)\, \mathbf1\{T(X,Y(\cdot),z)\le T(X,Y(\cdot),Z)\}.\]

Why can we use these known assignment probabilities? Initially the weight for an assignment \(z\) is

\[\operatorname{\mathsf P}(Z^*=z\mid X,Y(\cdot),Z) =\operatorname{\mathsf P}(Z=z\mid X,Y(\cdot)) =\operatorname{\mathsf P}(Z=z\mid X).\]

The first equality uses the conditional independence and identical distributions of \(Z\) and \(Z^*\). The second uses the randomization assumption. Knowing the assignment law lets us evaluate these weights, so the calculation uses both parts of that assumption.

Theorem 3.2 (validity).

Under the randomization assumptions, for every \(0\le\alpha\le1\),

\[\operatorname{\mathsf P}\{Q\le\alpha\mid Y(\cdot),X\}\le\alpha.\]

Proof.

The proof uses a general fact about cumulative distribution functions. If a real-valued random variable \(T\) has CDF \(F\), then

\[\operatorname{\mathsf P}\{F(T)\le\alpha\}\le\alpha, \qquad 0\le\alpha\le1.\]

For a continuous distribution, \(F(T)\) is uniform on \([0,1]\). With atoms, it is super-uniform: its probability of falling below \(\alpha\) is at most \(\alpha\). The example sheet asks you to establish this fact.

Conditional on \(X,Y(\cdot)\), the observed and reference statistics have the same distribution. Our \(Q\) is its CDF evaluated at the observed statistic, so applying this fact conditionally proves the theorem.

3.1.3 When can we calculate it?

There is still an apparent problem: \(Q\) depends on \(Y(\cdot)\), and the fundamental problem of causal inference (Section 1.4.4) says that we do not observe the whole schedule. Under a sharp null, however, we can impute it. For example, under \(Y_i(0)=Y_i(1)\), consistency gives

\[Y_i(0)=Y_i(1)=Y_i.\]

Each missing potential outcome equals the observed outcome. We can therefore evaluate the statistic under every possible assignment and calculate \(Q\). Rejecting when \(Q\le\alpha\) gives a level-\(\alpha\) test, conditional on the covariates and potential outcomes of these particular units. No large-sample approximation is involved; discreteness can make the test conservative.

3.2 More precise inference: adjustment and blocking

So far, our inference has depended on the randomization assumption. Let’s now turn to the design itself. In The Design of Experiments (1935), Fisher emphasized three principles:

1. Randomization: assign treatments randomly to experimental units.
2. Replication: apply an experimental treatment to multiple independent experimental units, rather than trying to infer its effect from one observation.
3. Blocking: group similar units before assigning treatments, then randomize within those groups.

Blocking aims to reduce variation by comparing more homogeneous units. Before we see how it works, consider what we can do with a completely randomized experiment and a better analysis.

3.2.1 Three least-squares estimators

Suppose exactly \(n_1\) out of \(n\) units receive treatment, with \(A_i=Z_i\). We observe a pre-treatment covariate \(X_i\). Take \(X_i\) to be one-dimensional to keep the notation simple; the same idea extends to several covariates.

The ANOVA estimator fits an intercept and a treatment term:

\[\begin{aligned} (\hat\alpha_1,\hat\beta_1) &=\underset{\alpha,\beta}{\operatorname{argmin}} \sum_{i=1}^n(Y_i-\alpha-\beta A_i)^2. \end{aligned}\]

Here \(\hat\beta_1\) is just the difference between the treated and control means. In R formula notation, the fit is Y ~ 1 + A.

If \(X\) predicts \(Y\), some outcome variation can be explained by \(X\). Complete randomization assigns treatment independently of these covariates, so we may be able to remove some of that variation from the comparison. Analysis of covariance (ANCOVA) adds a covariate term:

\[\begin{aligned} (\hat\alpha_2,\hat\beta_2,\hat\gamma_2) &=\underset{\alpha,\beta,\gamma}{\operatorname{argmin}} \sum_{i=1}^n(Y_i-\alpha-\beta A_i-\gamma X_i)^2. \end{aligned}\]

Its formula is Y ~ 1 + A + X.

In ANCOVA, are we trying to test \(\gamma\) rather than \(\beta\)? No: our target is still the treatment effect. For now we are constructing estimates, and \(\hat\beta_2\), like \(\hat\beta_1\), estimates that effect. The covariate term helps adjust the comparison.

A third approach adds an interaction between treatment and the covariate. Write \(\bar X=n^{-1}\sum_iX_i\). Analysis of heterogeneous covariance (ANHECOVA) solves

\[\begin{aligned} &(\hat\alpha_3,\hat\beta_3,\hat\gamma_3,\hat\delta_3)\\ &\quad=\underset{\alpha,\beta,\gamma,\delta}{\operatorname{argmin}} \sum_{i=1}^n\bigl\{Y_i-\alpha-\beta A_i-\gamma(X_i-\bar X) -\delta A_i(X_i-\bar X)\bigr\}^2. \end{aligned}\]

The fit allows a different covariate slope in each treatment group. Why centre \(X\)? The fitted treated-minus-control contrast at \(X_i\) is \(\hat\beta_3+\hat\delta_3(X_i-\bar X)\). Averaging gives

\[\frac1n\sum_i\{\hat\beta_3+\hat\delta_3(X_i-\bar X)\} =\hat\beta_3,\]

because the centred covariates sum to zero. Thus the treatment coefficient itself is the average fitted contrast.

3.2.2 When adjustment helps—and when it hurts

Consider a simulation with 100 treated units and 400 controls. Generate independent standard normal \(X,\epsilon_0,\epsilon_1\) for each unit and set

\[Y(0)=X+\epsilon_0,\qquad Y(1)=2+X+\epsilon_1.\]

The population average treatment effect is 2. Each repetition draws a new sample and randomizes treatment. Look first at the upper panels, which show the same single experiment fitted in three ways. ANOVA uses two horizontal lines, one at each treatment group’s mean. ANCOVA fits parallel lines whose vertical separation estimates the treatment effect. ANHECOVA allows the slopes to differ and compares the fitted lines at \(\bar X\). Here the true slopes are equal, so the last two fits are similar.

The lower panels ask a different question: how much does each estimate vary when we repeat the experiment? They show the distributions over 50,000 experiments.

Parallel slopes: ANOVA fits horizontal lines, ANCOVA fits parallel lines, and ANHECOVA allows separate slopes. The lower panels show repeated-experiment distributions, with the population ATE of 2 marked.

On a small screen, scroll the figure horizontally to read the labels.

All three distributions are centred near 2, but adjustment reduces their spread. The standard deviations are 0.157 for ANOVA and 0.112 for each adjusted estimator.

Now change just the treated slope:

\[Y(0)=X+\epsilon_0,\qquad Y(1)=2-X+\epsilon_1.\]

The population ATE is still 2, but the effect varies with \(X\). ANHECOVA can pick up that heterogeneity; ANCOVA forces the slopes to agree.

Opposing slopes: forcing parallel fitted lines can make adjustment less precise than the unadjusted comparison. Allowing separate slopes retains a precision gain in this example.

On a small screen, scroll the figure horizontally to read the labels.

The standard deviations are now 0.158, 0.197 and 0.143. Adjusting for a covariate can make the estimate less precise. This is a deliberately stark example; ANCOVA often works well, but its efficiency gain is not guaranteed.

Large-sample comparison. Under the usual regularity conditions for regression adjustment in randomized experiments, the asymptotic variance of sample-centred ANHECOVA is no larger than that of either ANOVA or ANCOVA. There is no universal ordering between ANOVA and ANCOVA. This is an asymptotic, non-strict comparison, not a promise that ANHECOVA wins in every finite sample; see Reluga, Ye and Zhao (2022), §3.2.

Fisher developed ANOVA and ANCOVA, and both became standard tools. The lack of a general efficiency guarantee for ordinary ANCOVA received renewed attention in the 2000s, motivating the use of treatment–covariate interactions.

3.2.3 Blocking: improve the design first

The appealing feature of blocking is that we can obtain a more precise comparison with a simple analysis. The name comes from agricultural experiments: divide a field into groups of similar plots and compare treatments within those groups, reducing the contribution of between-group variation.

A Latin square blocks on two factors at once. Imagine four rows and four columns of field plots, with four treatments labelled A–D. Arrange them so that each treatment appears exactly once in each row and once in each column, then randomize within these constraints.

A four-by-four Latin square. Every treatment occurs once in each row and each column.

Rows and columns may each affect crop yield. A full factorial arrangement of row, column and treatment would contain \(4\times4\times4=64\) combinations. The Latin square uses 16, balancing treatment across both blocking factors. In this sense it is a special fractional factorial design.

Fisher’s Latin-square window at Caius.¹

Fisher’s memorial window at Gonville and Caius College depicted a seven-by-seven Latin square. The college removed the window in 2020 following controversy over his involvement in eugenics. This is a complicated history: eugenics did not mean entirely the same thing then as it does today. The college’s Fisher in the 21st Century discussion and talks explore that history.

Return to the opposing-slopes simulation. This time, change the assignment:

1. Observe the 500 covariate values and sort the units by \(X\).
2. Form 100 blocks of five consecutive units.
3. Within each block, randomly assign one unit to treatment and four to control.
4. Estimate the effect by the ordinary difference in means.

The same opposing-slopes setting, now comparing two designs with the same unadjusted estimator. Blocking reduces the standard deviation from 0.158 to 0.143.

On a small screen, scroll the figure horizontally to read the labels.

In this example, blocking followed by a difference in means reaches essentially the same precision as ANHECOVA under complete randomization. We have changed the design rather than complicated the analysis.

When explaining the result to a medical regulator, a careful blocked design and a simple estimate can be more compelling than a less careful design followed by regression adjustments and their asymptotic justification. Better design and simpler analysis can make the practical argument easier too.

3.3 Stepped wedge designs

3.3.1 A change introduced over time

Sometimes cost or logistics prevent an intervention from being introduced everywhere at once. A stepped wedge design randomizes the order in which groups change condition. It can also be used to withdraw a service.

Consider the Australian study by Haines et al. (2017). The study first withdrew an existing weekend allied health service, then introduced a new service model. We focus on the first trial at Dandenong Hospital: six wards, 7,986 admissions, and the outcome \(Y=\log(\text{patient length of stay})\). The practical concern was whether withdrawing weekend services would prolong patients’ stays.

The study’s original design figure, with Dandenong trial 1 highlighted and enlarged. White cells retain the existing weekend service; light grey cells have it withdrawn. Cell entries count weekday/weekend service occasions, not admissions.²

On a small screen, scroll the figure horizontally to read the labels.

3.3.2 Assignment is not the same as exposure

Let \(Z_w\) be the crossover month for ward \(w\). One ward switches in each of months 2–7, and the crossover order is randomized:

\[Z\text{ is uniform over the }6!=720 \text{ permutations of }(2,3,4,5,6,7).\]

For a patient \(i\) entering ward \(w_i\) in month \(t_i\), use the simple exposure map

\[A_i(Z)=\mathbf1\{t_i\ge Z_{w_i}\},\]

where \(A_i=1\) means that the ward’s routine weekend service has been withdrawn at admission.

A patient admitted just before the crossover might stay through it and experience both conditions. There are also some weekend service occasions after withdrawal. Perhaps wards needed time to adjust; exact admission timing and time since withdrawal might matter. Our simple exposure map ignores these possibilities. The assignment is specified by randomization, but translating it into exposure requires a modelling decision.

3.3.3 Randomize the crossover order

Under the sharp null, every patient’s outcome stays fixed under every possible schedule. An exact randomization test (Theorem 3.2) therefore proceeds by changing the crossover order, recomputing each patient’s exposure, and recalculating the statistic. All 720 schedules have equal weight.

The reanalysis by Zhang and Zhao (2023), §2 uses the coefficient of \(A\) in two least-squares fits:

\[Y\sim1+A,\qquad Y\sim1+A+\text{ward},\]

where ward is a categorical variable. The second fit adjusts for differences between wards. Here larger coefficients are evidence in the direction of longer stays after withdrawal, so we use upper-tail randomization \(p\)-values, counting statistics at least as large as the observed one.

Exact randomization distributions over all 720 schedules under the sharp null. Rust lines mark the observed coefficients.

On a small screen, scroll the figure horizontally to read the labels.

Statistic	Observed coefficient	Upper-tail \(p\)-value
Unadjusted	0.290	\(60/720\approx0.0833\)
Ward-adjusted	0.164	\(3/720\approx0.0042\)

There is a curious feature: all 720 ward-adjusted coefficients are positive. Length of stay increased over time, which was correlated with service withdrawal. The study runs into the Australian winter, so perhaps seasonal pressures contributed to longer stays; this is only a conjecture. The positive coefficients do not mean that withdrawal must have a positive treatment effect: a randomization distribution need not be centred at zero.

Could we instead permute the patients’ admission times, or which wards they enter? We could calculate such quasi-randomization tests, but those permutations are not justified by the trial’s randomization. The wards are also quite different. The example sheet will let you explore these alternatives.

3.4 Group sequential designs

3.4.1 Why not stop when the evidence is clear?

Trials can be long and expensive. If the evidence already shows that a treatment is effective, why keep spending money? There may also be an ethical reason to make an effective treatment available sooner.

The difficulty is that looking repeatedly gives us repeated opportunities to reject. A test calibrated for one look does not automatically control error when we keep looking and stop as soon as a result is significant. Let’s first set up the general multiple-testing problem.

Definition 3.3.

Suppose the data have distribution \(\mathsf P_\theta\), where \(\theta\in\Theta\) is unknown. Represent the hypotheses by functions

\[H_1,\ldots,H_m:\Theta\to\{0,1\},\]

with \(H_j(\theta)=0\) when the \(j\)th null hypothesis is true and \(H_j(\theta)=1\) otherwise. The \(p\)-value \(Q_j\) is valid if

\[\operatorname{\mathsf P}_\theta(Q_j\le u)\le u, \qquad 0\le u\le1,\quad \theta\text{ such that }H_j(\theta)=0.\]

Definition 3.4.

Reject hypothesis \(j\) when \(Q_j\le\alpha_j\). The family-wise error rate (FWER) is the probability of at least one false rejection:

\[\operatorname{FWER}_\theta =\operatorname{\mathsf P}_\theta\left\{ \sum_{j=1}^m[1-H_j(\theta)]\mathbf1\{Q_j\le\alpha_j\}\ge1\right\}.\]

The factor \(1-H_j(\theta)\) selects the true nulls.

3.4.2 Bonferroni’s correction

Proposition 3.5.

If each \(Q_j\) is valid for its own null, taking \(\alpha_j=\alpha/m\) controls the FWER at \(\alpha\), for every \(\theta\in\Theta\).

Proof.

Fix \(\theta\). A false rejection occurs when we reject at least one of the nulls for which \(H_j(\theta)=0\). Apply the union bound to precisely those nulls:

\[\begin{aligned} \operatorname{FWER}_\theta &=\operatorname{\mathsf P}_\theta \left(\bigcup_{j:H_j(\theta)=0}\{Q_j\le\alpha/m\}\right)\\ &\le\sum_{j:H_j(\theta)=0} \operatorname{\mathsf P}_\theta(Q_j\le\alpha/m)\\ &\le\sum_{j:H_j(\theta)=0}\frac\alpha m\le\alpha. \end{aligned}\]

There are at most \(m\) true nulls. We need neither independence of the tests nor an assumption that all the nulls are true.

3.4.3 Five looks: new data or all the data?

Consider independent \(Y_i\sim N(\mu,1)\), with five planned looks, and test \(H_0:\mu\le0\). At look \(k\), use a standard-normal statistic under \(\mu=0\) and an upper-tail \(p\)-value

\[Q_k=1-\Phi(T_k),\]

where \(\Phi\) is the standard normal CDF. Reject and stop at the first \(Q_k\le0.05/5=0.01\).

We could use only the new observation at each look:

\[T_k=Y_k.\]

These five tests are independent. By the fifth test, their overall rejection probability is 4.9% under \(\mu=0\), and their power is 38.4% under \(\mu=1\).

But why throw away data we already have? Alternatively, accumulate all the observations so far:

\[T_k=\frac{Y_1+\cdots+Y_k}{\sqrt{k}},\qquad k=1,\ldots,5.\]

Use the same 1% cutoff at every look. Now the overall rejection probability is 3.1% under \(\mu=0\), and power is 54.2% under \(\mu=1\).

Bonferroni with five independent tests versus five looks at accumulating data. Each curve gives the probability of at least one rejection by that look; rust is FWER at mean zero and teal is power at mean one.

On a small screen, scroll the figure horizontally to read the labels.

Using all the data improves power, but the tests are positively correlated because they reuse observations. Their rejection events overlap, making the union bound in the proof of Proposition 3.5 conservative. We use only 3.1 percentage points of the permitted 5% error rate. Can we use the remaining allowance to gain power?

3.4.4 Spend error probability on the first rejection

In a group sequential design, each look tests the same hypothesis about the same treatment and population. Let \(\Theta_0\) be the null set. In the notation above,

\[H_1(\theta)=\cdots=H_m(\theta)=\mathbf1\{\theta\notin\Theta_0\}.\]

Observe statistics \(T_1,\ldots,T_m\). At look \(k\), reject and stop if \(T_k>b_k\). If no boundary is crossed by look \(m\), finish without rejecting.

To reject first at look \(k\), we must have continued through all earlier looks:

\[T_1\le b_1,\ldots,T_{k-1}\le b_{k-1},\quad T_k>b_k.\]

These first-rejection events are disjoint, so their probabilities add without the overcounting in the union bound.

Definition 3.6.

Choose a non-decreasing \(\alpha\)-spending function \(\alpha(t)\) with

\[\alpha(0)=0,\qquad \alpha(1)=\alpha.\]

Here \(\alpha\) without an argument is the total error allowance, such as 0.05. Let \(t_k\) denote the information fraction at look \(k\), with \(0=t_0<t_1<\cdots<t_m=1\). For equally informative observations, this fraction can be proportional to the sample size; in our five-look example, \(t_k=k/5\).

Choose the boundaries using the joint distribution of the statistics so that

\[\begin{aligned} &\sup_{\theta\in\Theta_0}\operatorname{\mathsf P}_\theta (T_1\le b_1,\ldots,T_{k-1}\le b_{k-1},T_k>b_k)\\ &\hspace{4em}\le\alpha(t_k)-\alpha(t_{k-1}). \end{aligned}\tag{3.1}\]

When possible, calibrate with equality. The increment is the allowance for the first rejection at look \(k\), not the marginal probability \(\operatorname{\mathsf P}_\theta(T_k>b_k)\) with earlier looks ignored.

Proposition 3.7 (error control).

Boundaries satisfying Equation (3.1) control the probability of rejecting by the final look at \(\alpha\), for every \(\theta\in\Theta_0\).

Proof.

Indeed, the possible first-rejection times partition the rejection event. For any \(\theta\in\Theta_0\),

\[\begin{aligned} \operatorname{\mathsf P}_\theta(\text{reject}) &=\sum_{k=1}^m\operatorname{\mathsf P}_\theta(\text{first reject at look }k)\\ &\le\sum_{k=1}^m\{\alpha(t_k)-\alpha(t_{k-1})\}\\ &=\alpha(1)-\alpha(0)=\alpha. \end{aligned}\]

The sum telescopes. Rather than separately bounding overlapping events, we have used their joint distribution to allocate error probability to disjoint first crossings.

3.4.5 Two choices of spending function

There is no universally best spending function: the choice depends on what we want the design to do. Two popular choices are Pocock-type and O’Brien–Fleming-type spending.

For the accumulating normal statistics, the boundaries below are calibrated so that the cumulative probability of crossing by look \(k\) is

\[\operatorname{\mathsf P}_0\left\{ \max_{j\le k}(T_j-b_j)>0\right\}=\alpha(k/5), \qquad \alpha(1)=0.05.\]

Here 0 denotes \(\mu=0\), the least favourable null value for an upper-boundary crossing. The left panel shows cumulative error spending \(\alpha(t)\); the right shows the resulting thresholds \(b_k\) on the statistic’s scale.

Spending functions and their corresponding boundaries for five looks. O’Brien–Fleming-type spending reserves almost all the allowance early on, requiring much stronger evidence for an early rejection.

On a small screen, scroll the figure horizontally to read the labels.

Pocock-type spending permits more error probability early and gives relatively even thresholds. O’Brien–Fleming-type spending reserves more of the allowance for later: an early rejection requires much stronger evidence, but the threshold at the final look is lower.

Cumulative error and power for the two spending rules in the same normal experiment. Both finish at 5% FWER; their power differs across looks. Values come from deterministic normal-probability calculations.

On a small screen, scroll the figure horizontally to read the labels.

At \(\mu=1\), Pocock-type spending gives more power at the first four looks. O’Brien–Fleming-type spending overtakes at the fifth, finishing at 70.9% power compared with 63.8%. Both use the full 5% error allowance. The comparison makes the design choice visible: an opportunity to stop earlier versus more power at the final look.

Further reading

• Fisher, The Design of Experiments (1935), especially Chapter V on Latin squares.
• Reluga, Ye and Zhao (2022), A unified analysis of regression adjustment in randomized experiments, §§2–3.
• Zhang and Zhao (2023), What is a randomization test?, §1.3 on Fisher’s exact test and §§2–3 on the hospital example and the general framework.
• O’Brien and Fleming (1979), A Multiple Testing Procedure for Clinical Trials, and Lan and DeMets (1983), Discrete Sequential Boundaries for Clinical Trials.

© 2026 Qingyuan Zhao. Except where otherwise noted, these lecture notes are licensed under Creative Commons Attribution 4.0 International (CC BY 4.0). Third-party material is excluded from this licence; see the individual credits and rights notices. Public-domain material remains in the public domain.

---

1. Glass by Maria McClafferty; photograph by Schutz (2006), straightened by McGeddon. Image source, CC BY-SA 2.5. The college’s 29 July 2020 notice records the removal.

2. Haines et al. (2017), Figure 1, CC BY 4.0. Highlight and enlarged crop added for these teaching materials.