2 Randomization

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

2.1 Randomization makes population inference possible

Let’s return to the question at the end of Lecture 1: why is randomization useful? We can observe at most one potential outcome for each person, so an individual treatment effect remains out of reach. But we can learn about an average treatment effect. Let’s see how, starting with a simple experiment.

2.1.1 A fixed population and a random assignment

Recall that the treatment assignment is a vector \(Z\in\mathcal Z\), and the potential outcome schedule \(Y(\cdot):\mathcal Z\to\mathbb R^n\) describes the outcomes under each possible assignment. In the Neyman–Rubin model, \(\mathcal Z=\mathcal A^n\) and the exposure map is \(A_i(z)=z_i\): unit \(i\)’s outcome depends on its own treatment, rather than the assignments to other units. Take \(\mathcal A=\{0,1\}\), so that each unit has potential outcomes \(Y_i(0)\) and \(Y_i(1)\), and consistency gives

\[Y_i=Z_iY_i(1)+(1-Z_i)Y_i(0).\]

For \(a=0,1\), define the potential-outcome means and variances

\[\bar Y(a)=\frac1n\sum_{i=1}^nY_i(a),\qquad S_a^2=\frac1{n-1}\sum_{i=1}^n\{Y_i(a)-\bar Y(a)\}^2.\]

The average treatment effect over these \(n\) units is \(\bar Y(1)-\bar Y(0)\). The variance of their individual treatment effects is

\[S_{01}^2=\frac1{n-1}\sum_{i=1}^n \bigl[\{Y_i(1)-Y_i(0)\}-\{\bar Y(1)-\bar Y(0)\}\bigr]^2.\]

These are ordinary sample means and variances, applied to the potential-outcome table. We cannot calculate them all from the observed data: each row of that table has a missing entry.

Under complete randomization, choose exactly \(n_1\) of the \(n\) units to receive treatment, with every such subset equally likely. Put \(n_0=n-n_1\), with both group sizes positive. The assignment law is

\[\operatorname{\mathsf P}\{Z=z\mid Y(\cdot)\}= \begin{cases} \binom n{n_1}^{-1},&z\in\{0,1\}^n,\ \sum_i z_i=n_1,\\ 0,&\text{otherwise}. \end{cases}\]

Is \(n_1\) random? No: it is part of this design. For example, with 30 people we might decide to treat exactly 10. Which 10 receive treatment is random; the number 10 is fixed.

2.1.2 The difference-in-means estimator

Take the average outcome among treated units, subtract the average among control units, and call the result \(\hat\beta\):

\[\hat\beta=\frac1{n_1}\sum_{i=1}^n Z_iY_i -\frac1{n_0}\sum_{i=1}^n(1-Z_i)Y_i.\]

Theorem. In the Neyman–Rubin model under complete randomization,

\[\operatorname{\mathsf E}\{\hat\beta\mid Y(\cdot)\}=\bar Y(1)-\bar Y(0),\]

\[\operatorname{Var}\{\hat\beta\mid Y(\cdot)\} =\frac{S_1^2}{n_1}+\frac{S_0^2}{n_0}-\frac{S_{01}^2}{n}.\]

Let’s prove the first identity. We are conditioning on the potential outcomes, but the estimator is written in terms of the observed outcomes. Consistency connects the two: if \(Z_i=1\), then \(Y_i=Y_i(1)\); if \(Z_i=0\), the term \(Z_iY_i\) contributes nothing. Thus \(Z_iY_i=Z_iY_i(1)\), and similarly \((1-Z_i)Y_i=(1-Z_i)Y_i(0)\).

Once we make that substitution, the potential outcomes are fixed and only the \(Z_i\)’s are random. Each unit has treatment probability \(n_1/n\), so \(\operatorname{\mathsf E}\{Z_i\mid Y(\cdot)\}=n_1/n\) and \(\operatorname{\mathsf E}\{1-Z_i\mid Y(\cdot)\}=n_0/n\). By linearity of conditional expectation,

\[\begin{aligned} \operatorname{\mathsf E}\{\hat\beta\mid Y(\cdot)\} &=\operatorname{\mathsf E}\left\{ \frac1{n_1}\sum_iZ_iY_i(1)-\frac1{n_0}\sum_i(1-Z_i)Y_i(0) \;\middle|\;Y(\cdot)\right\}\\ &=\frac1{n_1}\sum_i\frac{n_1}{n}Y_i(1) -\frac1{n_0}\sum_i\frac{n_0}{n}Y_i(0)\\ &=\bar Y(1)-\bar Y(0). \end{aligned}\]

The calculation is short, but its meaning deserves some attention. We have conditioned on the whole schedule \(Y(\cdot)\): the causal law describing how treatment changes the outcomes. All the randomness comes from \(Z\). The experimenter generates it by rolling dice, for example, or using a pseudorandom number generator. These are devices we control, whose properties we can examine.

Accept the causal model and the specified randomization device, and we have an unbiased estimator of the average effect on these \(n\) units. We can also ask how precise it is. The variance calculation is left to the example sheet; if \(S_0^2\) and \(S_1^2\) stay bounded as both group sizes grow, the variance tends to zero. A larger experiment then gives a more precise estimate.

This is the sense in which randomization makes population inference possible. Here the population is the collection of \(n\) experimental units. We still cannot infer each person’s individual treatment effect, but we can estimate their average.

2.1.3 A practical example: Moderna’s vaccine trial

For much of 2020, there was a pandemic, we were staying at home, and we were rather bored. Towards the end of the year, several vaccine trials brought good news. Consider the question in Moderna’s trial: does vaccination prevent symptomatic COVID-19? How could we give a convincing answer, including to someone sceptical of vaccines?

In Moderna’s trial, participants were assigned to mRNA-1273 or placebo in a 1:1 ratio within predefined groups, or strata. They received two injections 28 days apart and did not know which treatment they had received. The primary outcome was symptomatic, laboratory-confirmed COVID-19 beginning at least 14 days after the second dose. The results for the per-protocol population were:

Group	Cases	Non-cases	Total	Risk
mRNA-1273	11	14,123	14,134	0.078%
Placebo	185	13,888	14,073	1.315%

Source: Baden et al. (2021; published online 30 December 2020), Figure 1 and primary efficacy results. Risks are calculated from the counts.

The risk ratio from this table is

\[\frac{11/14134}{185/14073}\approx0.0592,\]

so the relative risk reduction is approximately \(1-0.0592=94.1\%\). This is the scale of the widely reported “94% protection”.

To analyse a \(2\times2\) table like this, we can use Fisher’s exact test. In R:

tab <- rbind(vaccine = c(11, 14123),
             placebo = c(185, 13888))
fisher.test(tab)

The output gives \(p<2.2\times10^{-16}\), an estimated odds ratio of about \(0.0585\), and a 95% confidence interval approximately \((0.0287,0.1072)\) for the odds ratio.

We will return to how this test works. For now, notice what the trial’s design contributes to the argument: the vaccine and placebo were randomly assigned, and participants did not know their assignment. Randomized clinical trials have become a central way to gather medical evidence, with uses in the social sciences too.

Our first calculation gave a mathematical reason to study randomization; the vaccine trial gives a practical one. To examine its contribution more closely, we will follow Cox (2009), §2 and distinguish three roles:

Avoid hidden bias: identify causal effects using random assignment.
Justify analysis of variance: connect the design to measures of variation.
Construct exact tests: use the experiment’s known assignment distribution.

2.2 Role 1: avoid hidden bias

2.2.1 Causal identification

What do we mean mathematically by “avoiding hidden bias”? Consider units \((X_i,Y_i(0),Y_i(1),A_i)\), where \(X_i\) records pre-treatment measurements and \(A_i=Z_i\). For this calculation, we switch to a model in which the units are independent and identically distributed (i.i.d.). We can then drop the subscript \(i\) and work with one unit’s distribution. Our target is now a population average, rather than the average over a particular fixed collection of experimental units.

Allow randomization to depend on the pre-treatment measurements \(X\). Assume consistency, finite potential-outcome means, and

\[A\mathrel{\perp\!\!\!\perp}\{Y(0),Y(1)\}\mid X, \qquad 0<\operatorname{\mathsf P}(A=1\mid X)<1 \quad\text{almost surely}.\]

The first condition says that treatment is independent of the potential outcomes once we condition on \(X\). The second, called positivity, says that someone with those measurements has a chance of receiving either treatment. Assignment is never completely determined by \(X\).

Theorem (identification). Under these assumptions, the population average treatment effect satisfies

\[\begin{aligned} \operatorname{\mathsf E}\{Y(1)-Y(0)\} &=\operatorname{\mathsf E}\left\{ \operatorname{\mathsf E}(Y\mid A=1,X) -\operatorname{\mathsf E}(Y\mid A=0,X)\right\}\\ &=\operatorname{\mathsf E}\left\{ \frac{AY}{\operatorname{\mathsf P}(A=1\mid X)} -\frac{(1-A)Y}{\operatorname{\mathsf P}(A=0\mid X)}\right\}. \end{aligned}\]

Look at what is on each side. On the left is a quantity involving potential outcomes, which we cannot observe together. On the right are quantities involving \(X\), \(A\) and \(Y\), which we do observe. This is causal identification: the assumptions let us express the causal quantity in terms of the distribution of observed variables. We can then estimate it using empirical data. This idea will recur throughout the course.

2.2.2 Why the identification formula holds

Fix \(a\in\{0,1\}\). Begin with the law of total expectation:

\[\operatorname{\mathsf E}\{Y(a)\} =\operatorname{\mathsf E}\bigl[\operatorname{\mathsf E}\{Y(a)\mid X\}\bigr].\]

Conditional independence lets us add \(A=a\) inside; positivity makes this a possible treatment within the covariate group. Consistency then replaces \(Y(a)\) by \(Y\):

\[\begin{aligned} \operatorname{\mathsf E}\{Y(a)\} &=\operatorname{\mathsf E}\bigl[\operatorname{\mathsf E}\{Y(a)\mid X,A=a\}\bigr]\\ &=\operatorname{\mathsf E}\bigl\{\operatorname{\mathsf E}(Y\mid X,A=a)\bigr\}. \end{aligned}\]

Subtract the expression for \(a=0\) from that for \(a=1\) to obtain the first identification formula.

For the second equality, the trick is to multiply and divide by the treatment probability:

\[\begin{aligned} \operatorname{\mathsf E}\{Y(a)\} &=\operatorname{\mathsf E}\left\{ \frac{\operatorname{\mathsf P}(A=a\mid X)\operatorname{\mathsf E}(Y\mid X,A=a)} {\operatorname{\mathsf P}(A=a\mid X)}\right\}\\ &=\operatorname{\mathsf E}\left\{ \frac{\operatorname{\mathsf E}[\mathbf1\{A=a\}Y\mid X]} {\operatorname{\mathsf P}(A=a\mid X)}\right\}\\ &=\operatorname{\mathsf E}\left\{ \frac{\mathbf1\{A=a\}Y}{\operatorname{\mathsf P}(A=a\mid X)}\right\}. \end{aligned}\]

Why does the second line hold? Condition on \(A\) inside the numerator. The indicator \(\mathbf1\{A=a\}\) is zero unless \(A=a\), leaving the probability of that event times the conditional mean of \(Y\). For the last line, the denominator is a function of \(X\), so we can move it inside the inner expectation and apply total expectation again.

This is inverse-probability weighting (IPW). We use outcomes from units receiving treatment \(a\), weighting each by the inverse of its probability of receiving that treatment. The formula makes “avoiding hidden bias” precise: there may be outcome predictors we never measure, but identification still holds when assignment is conditionally independent of the potential outcomes.

2.3 Role 2: justify analysis of variance

2.3.1 One-way analysis of variance

Let’s recall one-way analysis of variance (ANOVA), beginning with the usual textbook model. Take binary exposures \(A_i\) under complete randomization, with fixed group sizes \(n_1,n_0\), and suppose

\[Y_i=\alpha+\beta A_i+\varepsilon_i,\qquad \varepsilon_i\overset{\mathrm{i.i.d.}}{\sim}N(0,\sigma^2),\quad \sigma^2>0,\]

with the errors independent of the assignment. Write the observed group means as \(\bar Y_1,\bar Y_0\) and the overall mean as \(\bar Y=n^{-1}\sum_iY_i\). The subscripts distinguish these observed group means from the potential-outcome means \(\bar Y(1),\bar Y(0)\).

When we put the two groups together, we see variation in their outcomes. How much is between the groups, and how much is within them? The between-group variation may reflect a treatment effect; within-group variation might arise from measurement error, for example. ANOVA separates these two components.

Start with the total sum of squares, \(S_T=\sum_i(Y_i-\bar Y)^2\). Split the sum into the two groups, then add and subtract each group’s mean:

\[\begin{aligned} S_T &=\sum_i A_i(Y_i-\bar Y_1+\bar Y_1-\bar Y)^2\\ &\quad+\sum_i(1-A_i)(Y_i-\bar Y_0+\bar Y_0-\bar Y)^2. \end{aligned}\]

Expand the squares. The cross terms vanish because, within each group, \(\sum_{i:A_i=a}(Y_i-\bar Y_a)=0\). Hence

\[\underbrace{\sum_i(Y_i-\bar Y)^2}_{S_T} =\underbrace{\sum_i(Y_i-\bar Y_{A_i})^2}_{S_E} +\underbrace{n_1(\bar Y_1-\bar Y)^2+n_0(\bar Y_0-\bar Y)^2}_{S_A}.\]

Here \(S_T\) is the total sum of squares, \(S_E\) the within-group (residual) sum of squares, and \(S_A\) the between-group (treatment) sum of squares. Because \(\bar Y=(n_1\bar Y_1+n_0\bar Y_0)/n\), the last term simplifies to

\[S_A=\frac{n_1n_0}{n}(\bar Y_1-\bar Y_0)^2 =\frac{n_1n_0}{n}\hat\beta^2.\]

2.3.2 Gaussian justification of the \(F\)-test

Under the Gaussian model, conditional on the assignment,

\[\frac{\hat\beta-\beta}{\sigma\sqrt{1/n_1+1/n_0}}\sim N(0,1), \qquad \frac{S_E}{\sigma^2}\sim\chi^2_{n-2},\]

and these two quantities are independent. Under \(H_0:\beta=0\), the square of the first is \(S_A/\sigma^2\sim\chi^2_1\). Dividing the independent chi-squared quantities by their degrees of freedom gives

\[F=\frac{S_A/1}{S_E/(n-2)}\sim F_{1,n-2},\qquad n>2.\]

Source	Sum of squares	Degrees of freedom	Mean square
Treatment	\(S_A\)	\(1\)	\(S_A\)
Residual	\(S_E\)	\(n-2\)	\(S_E/(n-2)\)
Total	\(S_T\)	\(n-1\)

Before computers, this table organized the calculation: compute the sums of squares, divide by their degrees of freedom, then divide the treatment mean square by the residual mean square to get \(F\). The same distributional result also gives the usual two-sample \(t\)-test.

2.3.3 Randomization justification

Portrait of the young Ronald Aylmer Fisher in 1913. — R. A. Fisher, 1913.¹

This is the textbook justification, but the historical development was different. Fisher first studied the decomposition of genetic variation. When he moved to agricultural experiments, he saw an analogy: use random allocation of fertilizers or crops across a field, much as nature randomizes genes, and apply the same analysis of variance. For Fisher, randomization was also a way to justify the analysis itself.

Return to complete randomization with the potential-outcome schedule fixed. The sharp null says that treatment has no effect on any unit:

\[Y_i(1)=Y_i(0)\quad\text{for every }i.\]

Theorem. Under complete randomization and the sharp null,

\[\operatorname{\mathsf E}\{S_A\mid Y(\cdot)\}=S^2, \qquad \operatorname{\mathsf E}\{S_E\mid Y(\cdot)\}=(n-2)S^2,\]

where \(S^2=S_0^2=S_1^2\) is the sample variance of the common potential outcomes \(Y_i(0)=Y_i(1)\).

We already have what we need to prove this. Under the sharp null, \(\bar Y(1)-\bar Y(0)=0\), so \(\hat\beta\) has conditional mean zero. Its conditional second moment is therefore its variance. Using the first theorem,

\[\begin{aligned} \operatorname{\mathsf E}\{S_A\mid Y(\cdot)\} &=\frac{n_1n_0}{n}\operatorname{\mathsf E}\{\hat\beta^2\mid Y(\cdot)\}\\ &=\frac{n_1n_0}{n}\left(\frac{S_1^2}{n_1}+\frac{S_0^2}{n_0} -\frac{S_{01}^2}{n}\right)\\ &=S^2. \end{aligned}\]

In the last line, \(S_{01}^2=0\) because every individual effect is zero, and \(S_0^2=S_1^2=S^2\). The total sum of squares is \(S_T=(n-1)S^2\), since the observed outcomes do not change with assignment. Subtracting \(S_A\) gives

\[\operatorname{\mathsf E}\{S_E\mid Y(\cdot)\}=(n-2)S^2.\]

Divide by the degrees of freedom, and both mean squares have expectation \(S^2\). This matches the Gaussian calculation without assuming Gaussian outcomes. Higher moments can also be compared to understand how closely the randomization distribution resembles the Gaussian reference. The expectations above agree exactly; an exact \(F\) distribution does not follow from them.

Fisher also had a geometric way of thinking about this. His daughter, Joan Fisher Box, described it in 1980:

“He could picture the distribution of n results as a pattern in n-dimensional space, and he could see that randomization would produce a symmetry in that pattern rather like that produced by a kaleidoscope, and which approximated the required spherical symmetry available, in particular, from standard normal theory assumptions.”

— Joan Fisher Box, R. A. Fisher and the Design of Experiments, 1922–1926.

The kaleidoscope image suggests how randomization could supply approximate symmetry of the kind used in normal-theory calculations.

2.4 Role 3: construct exact tests

2.4.1 The two-binomials problem

Return to a binary outcome such as a COVID-19 case. Let \(N_{ay}\) count units with exposure \(A=a\) and outcome \(Y=y\). A dot denotes summation over that index:

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

Under complete randomization, the row totals \(n_0,n_1\) are fixed. Once we know the number of cases in a row, we also know its number of non-cases, so only two counts are free to vary. There are several ways to model this table. In the two-binomials problem, we take

\[N_{01}\sim\operatorname{Bin}(n_0,\pi_0),\qquad N_{11}\sim\operatorname{Bin}(n_1,\pi_1),\qquad N_{01}\mathrel{\perp\!\!\!\perp}N_{11}.\]

Each participant has a Bernoulli outcome, with case probability \(\pi_0\) or \(\pi_1\) according to their treatment group. We test whether those probabilities are the same:

\[H_0:\pi_0=\pi_1\qquad\text{against}\qquad H_1:\pi_0\ne\pi_1.\]

In a joint distribution of exposure and binary outcome, equality of these conditional probabilities says that \(A\) and \(Y\) are independent.

2.4.2 Why condition on the total number of cases?

The likelihood is

\[L(\pi_0,\pi_1)= \binom{n_0}{N_{01}}\pi_0^{N_{01}}(1-\pi_0)^{N_{00}} \binom{n_1}{N_{11}}\pi_1^{N_{11}}(1-\pi_1)^{N_{10}}.\]

For \(0<\pi_0,\pi_1<1\), we can write this in exponential form:

\[\begin{aligned} L(\pi_0,\pi_1) &=\binom{n_0}{N_{01}}\binom{n_1}{N_{11}} (1-\pi_0)^{n_0}(1-\pi_1)^{n_1}\\ &\quad\times\exp\left\{ N_{01}\log\frac{\pi_0}{1-\pi_0} +N_{11}\log\frac{\pi_1}{1-\pi_1}\right\}. \end{aligned}\]

Now collect the terms in the exponent:

\[\begin{aligned} &N_{01}\log\frac{\pi_0}{1-\pi_0}+N_{11}\log\frac{\pi_1}{1-\pi_1}\\ &\quad=(N_{01}+N_{11})\log\frac{\pi_0}{1-\pi_0} +N_{11}\left(\log\frac{\pi_1}{1-\pi_1}-\log\frac{\pi_0}{1-\pi_0}\right). \end{aligned}\]

This is an exponential family with sufficient statistics \(N_{\cdot1}=N_{01}+N_{11}\) and \(N_{11}\). The natural parameter multiplying \(N_{11}\) is the log odds ratio. Testing \(\pi_0=\pi_1\) amounts to testing whether it is zero. The other parameter, the control log odds, is a nuisance parameter: we do not know it, but it is not what we are testing.

The classical approach is to condition on the sufficient statistic for that nuisance parameter. Here this means conditioning on the total number of cases, \(N_{\cdot1}\). Under \(H_0\), the conditional distribution is

\[\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}}},\]

for counts compatible with the fixed row totals. This is a hypergeometric distribution, with no unknown case probability in its formula.

Fisher’s exact test compares the observed treated-case count with this conditional null distribution. Unusually small or large counts provide evidence against equal case probabilities. “Exact” refers to using this finite-sample distribution rather than a large-sample approximation.

2.4.3 An optimality property

Classical exponential-family theory gives a two-sided conditional test that, with suitably chosen cutoffs and randomization at boundary counts, is uniformly most powerful among unbiased level-\(\alpha\) tests (UMPU).²

What does that mean? An unbiased level-\(\alpha\) test satisfies

\[\operatorname{\mathsf P}_{\pi_0,\pi_1}(\text{reject})\le\alpha \quad\text{whenever }\pi_0=\pi_1,\]

and

\[\operatorname{\mathsf P}_{\pi_0,\pi_1}(\text{reject})\ge\alpha \quad\text{whenever }\pi_0\ne\pi_1.\]

So its probability of rejection is at most \(\alpha\) when the null is true, and at least \(\alpha\) under every alternative. “Uniformly most powerful” means that, for every competing unbiased level-\(\alpha\) test and every alternative,

\[\begin{aligned} &\operatorname{\mathsf P}_{\pi_0,\pi_1}(\text{UMPU test rejects})\\ &\qquad\ge \operatorname{\mathsf P}_{\pi_0,\pi_1}(\text{competing test rejects}), \qquad \pi_0\ne\pi_1. \end{aligned}\]