Fisher randomization tests

Every previous chapter estimated something. This one tests something, and it is the only place in the tutorial where the guarantee is exact: no model, no asymptotics, no condition on the design beyond its being a design, and no condition whatsoever on the test statistic. If the treatment does nothing to anybody, then rejecting when the randomization p-value is at most $\alpha$ has type-I error at most $\alpha$. In finite samples. Always.

The price of that guarantee is the null hypothesis it is a guarantee about, and most of this chapter is about reading that hypothesis carefully enough to see why the argument works.¹Fisher’s own example was the lady tasting tea: eight cups, four with the milk poured first, and a null under which the lady’s guesses are unrelated to the cups. The reference distribution is the $(\genfrac{}{}{0ex}{}{8}{4})=70$ ways the four could have been chosen — which is exactly what the design is here, at a smaller $n$.

Events under a design

A Design n is a probability mass function on the $2^n$ assignments, so a set of assignments has a probability: add up the mass of the assignments in it.

/-- The **probability of an event** `A` under the design: the total probability of those
assignments that satisfy `A`.

`[DecidablePred A]` is Lean bookkeeping, not mathematics: to *filter* a finite set by a
property, Lean wants to know the property can be tested.  For real inequalities such as
`T z₀ ≤ T z` Mathlib supplies the instance automatically (classically), so in practice this
argument is invisible.

This is definitionally `Design.probOf`; the two names coexist because `probEvent` reads
better in the statements of this file. -/
noncomputable def probEvent (D : Design n) (A : Assignment n → Prop) [DecidablePred A] : ℝ :=
  ∑ z ∈ univ.filter A, D.prob z

Two things to read here. First, univ.filter A is the finite set of assignments satisfying A, so the whole definition is a finite sum — the same move as Design.expect, and the reason nothing in this file needs measure theory. Second, [DecidablePred A] is a typeclass argument: a piece of data Lean supplies invisibly, saying that membership in the filtered set can be decided. For the events we care about — statements like T z₀ ≤ T z about real numbers — Mathlib provides that instance classically, and the argument is never written by hand. It is bookkeeping, not mathematics, and it is worth recognizing on sight so that it does not look like a hypothesis when you are auditing a statement.

The library proves the properties you would expect and nothing you would not: probEvent_nonneg, probEvent_le_one, monotonicity in the event (probEvent_mono), invariance under logically equivalent descriptions (probEvent_congr), the complement rule probEvent_compl, and the small fact that an event containing one assignment of positive probability has positive probability (probEvent_pos_of_mem). The last of these is what makes randomization p-values strictly positive, so it earns its keep below.

/-- Exercise (warm-up).  **A probability is nonnegative.**

`D.probEvent A` is by definition `∑ z ∈ univ.filter A, D.prob z`, and `D.nonneg z`
says each term is nonnegative.  One lemma finishes it: `Finset.sum_nonneg`, whose
argument is a function taking an index *and a proof that the index is in the set*
to nonnegativity of that term. -/
theorem ex_probEvent_nonneg (D : Design n) (A : Assignment n → Prop) [DecidablePred A] :
    0 ≤ D.probEvent A := by
  sorry

/-- Exercise (warm-up).  **A probability is nonnegative.**

`D.probEvent A` is by definition `∑ z ∈ univ.filter A, D.prob z`, and `D.nonneg z`
says each term is nonnegative.  One lemma finishes it: `Finset.sum_nonneg`, whose
argument is a function taking an index *and a proof that the index is in the set*
to nonnegativity of that term. -/
theorem ex_probEvent_nonneg (D : Design n) (A : Assignment n → Prop) [DecidablePred A] :
    0 ≤ D.probEvent A := by
  exact Finset.sum_nonneg fun z _ => D.nonneg z

/-- Exercise (core).  **Monotonicity of probability.**  If every assignment satisfying
`A` also satisfies `B`, then `P(A) ≤ P(B)`.

The assignments contributing to `P(A)` are a *subset* of those contributing to
`P(B)`, and the extra terms are nonnegative:
`Finset.sum_le_sum_of_subset_of_nonneg` says exactly that.  Its first argument is
the subset relation between the two filtered sets, which you prove by `intro z hz`
and unfolding membership with `Finset.mem_filter`. -/
theorem ex_probEvent_mono (D : Design n) (A B : Assignment n → Prop)
    [DecidablePred A] [DecidablePred B] (h : ∀ z, A z → B z) :
    D.probEvent A ≤ D.probEvent B := by
  sorry

/-- Exercise (core).  **Monotonicity of probability.**  If every assignment satisfying
`A` also satisfies `B`, then `P(A) ≤ P(B)`.

The assignments contributing to `P(A)` are a *subset* of those contributing to
`P(B)`, and the extra terms are nonnegative:
`Finset.sum_le_sum_of_subset_of_nonneg` says exactly that.  Its first argument is
the subset relation between the two filtered sets, which you prove by `intro z hz`
and unfolding membership with `Finset.mem_filter`. -/
theorem ex_probEvent_mono (D : Design n) (A B : Assignment n → Prop)
    [DecidablePred A] [DecidablePred B] (h : ∀ z, A z → B z) :
    D.probEvent A ≤ D.probEvent B := by
  refine Finset.sum_le_sum_of_subset_of_nonneg ?_ (fun z _ _ => D.nonneg z)
  intro z hz
  rw [Finset.mem_filter] at hz ⊢
  exact ⟨hz.1, h z hz.2⟩

Fisher’s sharp null

Neyman’s null is $\tau =0$: the treatment does nothing on average. Fisher’s null is that the treatment does nothing to anybody.

/-- **Fisher's sharp null hypothesis of no effect whatsoever**: every unit's two potential
outcomes coincide.

Contrast with the weak (Neyman) null `P.tau = 0`, which only says the *average* of the
effects is zero and is compatible with large effects that cancel.  The sharp null is what
makes randomization inference possible, because it pins down every unobserved potential
outcome: if `y1 i = y0 i` then the outcome we did not see is equal to the one we did. -/
def SharpNull (P : Population n) : Prop := ∀ i, P.y1 i = P.y0 i

The two are not the same hypothesis, and the difference is the whole chapter. SharpNull P is a statement about $n$ pairs of numbers; $\tau =0$ is a statement about one average. Sharp implies weak (SharpNull.tau_eq_zero) and the converse fails: effects of $+1$ and $-1$ average to nothing.

What the sharp null buys is that the missing half of the table is not missing. If $y_1(i)=y_0(i)$ then the outcome you did not see equals the one you did, for every unit, so the observed vector does not depend on which assignment was drawn:

/-- Exercise (warm-up).  **Under the sharp null the observed outcome vector does not
depend on the assignment**: `Yobs P z = P.y0`, for every `z`.

This is the whole content of the sharp null, and everything else in the chapter
rests on it.  The conclusion is an equality of *functions*, so begin with
`funext i`; then `Yobs_eq_ite` turns the observed outcome into an `if`, the
hypothesis `h i : P.y1 i = P.y0 i` makes the two branches equal, and `ite_self`
collapses them. -/
theorem ex_SharpNull_Yobs_eq {P : Population n} (h : SharpNull P) (z : Assignment n) :
    Yobs P z = P.y0 := by
  sorry

/-- Exercise (warm-up).  **Under the sharp null the observed outcome vector does not
depend on the assignment**: `Yobs P z = P.y0`, for every `z`.

This is the whole content of the sharp null, and everything else in the chapter
rests on it.  The conclusion is an equality of *functions*, so begin with
`funext i`; then `Yobs_eq_ite` turns the observed outcome into an `if`, the
hypothesis `h i : P.y1 i = P.y0 i` makes the two branches equal, and `ite_self`
collapses them. -/
theorem ex_SharpNull_Yobs_eq {P : Population n} (h : SharpNull P) (z : Assignment n) :
    Yobs P z = P.y0 := by
  funext i
  rw [Yobs_eq_ite, h i, ite_self]

Imputation, and why it is not estimation

The experimenter never holds a Population. They hold one observed vector y. Under the sharp null that single vector determines the whole table:

/-- The population **imputed from an observed outcome vector under the sharp null**: both
potential outcomes of every unit are set equal to the outcome that was observed.

Under the sharp null this is not an estimate but the truth: if the two potential outcomes
agree then filling both of them in with the observed value reconstructs the population
exactly, which is the content of `SharpNull.imputeSharp_Yobs` below. -/
def imputeSharp (y : Fin n → ℝ) : Population n := ⟨y, y⟩

imputeSharp y is a Population whose two potential-outcome vectors are both y. It is what a randomization test computes with, and it has the property that makes replaying the experiment legitimate: whatever assignment you feed it, it gives back the data you started with.

/-- Exercise (warm-up).  **The imputed population reproduces the data under every
assignment**: `Yobs (imputeSharp y) z = y`, for all `z`.

Both potential outcomes of `imputeSharp y` are `y`, so the `if` that `Yobs`
performs has the same value in both branches.  The conclusion is an equality of
*functions*, so start with `funext i`; then `Yobs_eq_ite` exposes the `if`,
`imputeSharp_y1`/`imputeSharp_y0` replace the two branches, and `ite_self` closes
`(if c then a else a) = a`. -/
theorem ex_imputeSharp_Yobs (y : Fin n → ℝ) (z : Assignment n) :
    Yobs (imputeSharp y) z = y := by
  sorry

/-- Exercise (warm-up).  **The imputed population reproduces the data under every
assignment**: `Yobs (imputeSharp y) z = y`, for all `z`.

Both potential outcomes of `imputeSharp y` are `y`, so the `if` that `Yobs`
performs has the same value in both branches.  The conclusion is an equality of
*functions*, so start with `funext i`; then `Yobs_eq_ite` exposes the `if`,
`imputeSharp_y1`/`imputeSharp_y0` replace the two branches, and `ite_self` closes
`(if c then a else a) = a`. -/
theorem ex_imputeSharp_Yobs (y : Fin n → ℝ) (z : Assignment n) :
    Yobs (imputeSharp y) z = y := by
  funext i
  rw [Yobs_eq_ite, imputeSharp_y1, imputeSharp_y0, ite_self]

So absDiM y z (below) is not an approximation of anything: it is exactly the number the analyst would have reported had the design produced z instead of what it did.

That is the modest claim. The immodest one is that under the sharp null the imputed population is not a stand-in for the truth but is the truth:

/-- Exercise (core).  **Under the sharp null, imputation is exact**: the population
you get by filling both potential outcomes in with the observed vector *is* the
true population, whichever assignment was drawn.

`SharpNull.Yobs_eq` (the previous exercise) turns `Yobs P z₀` into `P.y0`, so the
goal becomes `imputeSharp P.y0 = P`, i.e. `⟨P.y0, P.y0⟩ = P`.  Take `P` apart with
`cases P with | mk y1 y0 => …` and finish by rewriting the first field with
`funext` and `h`. -/
theorem ex_imputeSharp_exact {P : Population n} (h : SharpNull P) (z₀ : Assignment n) :
    imputeSharp (Yobs P z₀) = P := by
  sorry

/-- Exercise (core).  **Under the sharp null, imputation is exact**: the population
you get by filling both potential outcomes in with the observed vector *is* the
true population, whichever assignment was drawn.

`SharpNull.Yobs_eq` (the previous exercise) turns `Yobs P z₀` into `P.y0`, so the
goal becomes `imputeSharp P.y0 = P`, i.e. `⟨P.y0, P.y0⟩ = P`.  Take `P` apart with
`cases P with | mk y1 y0 => …` and finish by rewriting the first field with
`funext` and `h`. -/
theorem ex_imputeSharp_exact {P : Population n} (h : SharpNull P) (z₀ : Assignment n) :
    imputeSharp (Yobs P z₀) = P := by
  rw [h.Yobs_eq z₀]
  cases P with
  | mk y1 y0 => exact congrArg (Population.mk · y0) (funext fun i => (h i).symm)

This is the sentence to remember when someone objects that randomization inference “assumes away” the missing data. It does not: it tests a hypothesis under which there is no missing data. The type of SharpNull.imputeSharp_Yobs says so — the conclusion is an equality of populations, not an approximation, a bound, or a limit.

Test statistics and the p-value

Fix any function of the assignment. Under the sharp null it is a known function: no unobserved potential outcome enters it, so its whole distribution can be enumerated.

/-- The **one-sided randomization p-value** of the observed assignment `z₀` for the test
statistic `T`: the probability, under the design, of drawing an assignment whose statistic is
at least as large as the one actually observed.

There is no null distribution to look up and no approximation: the reference distribution *is*
the design, the only randomness in the experiment.  In practice `T` is built from the observed
outcome vector via `imputeSharp`, e.g. `T = absDiM y` below; but the definition, and the
validity theorem, place no conditions on `T` at all. -/
noncomputable def pValue (D : Design n) (T : Assignment n → ℝ) (z₀ : Assignment n) : ℝ :=
  D.probEvent (fun z => T z₀ ≤ T z)

Unfolded, pValue D T z₀ is $\sum _{z:T(z_0)\le T(z)}Pr(z)$: a finite sum over assignments, with no reference distribution to look up. The statistic an analyst would actually use is built from the data through imputeSharp:

/-- The **two-sided difference-in-means statistic** attached to an observed outcome vector
`y`: `|DiM|` evaluated on the population imputed under the sharp null.

`absDiM y z` is the number the analyst would have computed if the assignment had been `z`.
Only `absDiM y z₀` at the realised assignment `z₀` is the actual estimate; the other values
are the counterfactual replications that make up the randomization distribution. -/
noncomputable def absDiM (y : Fin n → ℝ) (z : Assignment n) : ℝ := |DiM (imputeSharp y) z|

Note the shape of absDiM : (Fin n → ℝ) → Assignment n → ℝ. The first argument is the data — fixed once the experiment is over — and the second is the assignment, which ranges over everything the design could have done. Getting those two roles the right way round is the single most common error in an autoformalized randomization test, and it is visible in the type.

Two facts about pValue are worth proving before the theorem that matters.

/-- Exercise (core).  **A randomization p-value is never zero** when the observed
assignment could actually have been drawn.

The observed assignment is one of the assignments at least as extreme as itself,
because `T z₀ ≤ T z₀`.  `Design.probEvent_pos_of_mem` turns "this assignment
satisfies the event and has positive probability" into "the event has positive
probability"; the first argument is the witness that `z₀` satisfies the event,
which here is `le_refl (T z₀)`. -/
theorem ex_pValue_pos (D : Design n) (T : Assignment n → ℝ) (z₀ : Assignment n)
    (hz : 0 < D.prob z₀) : 0 < pValue D T z₀ := by
  sorry

/-- Exercise (core).  **A randomization p-value is never zero** when the observed
assignment could actually have been drawn.

The observed assignment is one of the assignments at least as extreme as itself,
because `T z₀ ≤ T z₀`.  `Design.probEvent_pos_of_mem` turns "this assignment
satisfies the event and has positive probability" into "the event has positive
probability"; the first argument is the witness that `z₀` satisfies the event,
which here is `le_refl (T z₀)`. -/
theorem ex_pValue_pos (D : Design n) (T : Assignment n → ℝ) (z₀ : Assignment n)
    (hz : 0 < D.prob z₀) : 0 < pValue D T z₀ := by
  exact D.probEvent_pos_of_mem (le_refl (T z₀)) hz

The consequence is worth stating in statistical language: the smallest p-value a randomization test can return is the probability of the assignment that happened. With four units and two treated there are six assignments, so no outcome of that experiment can produce , and nothing that happens in the data will make it significant at $0.05$. The design bounds the evidence before you look at anything.

/-- Exercise (core).  **The p-value is antitone in the observed statistic**: a more
extreme observation gets a smaller p-value.

Unfold nothing.  If `T z₀ ≤ T z₁` then every `z` with `T z₁ ≤ T z` also has
`T z₀ ≤ T z`, so the event defining `pValue D T z₁` is contained in the one
defining `pValue D T z₀`, and `Design.probEvent_mono` compares them.  The
containment is one `le_trans`. -/
theorem ex_pValue_antitone (D : Design n) (T : Assignment n → ℝ) {z₀ z₁ : Assignment n}
    (h : T z₀ ≤ T z₁) : pValue D T z₁ ≤ pValue D T z₀ := by
  sorry

/-- Exercise (core).  **The p-value is antitone in the observed statistic**: a more
extreme observation gets a smaller p-value.

Unfold nothing.  If `T z₀ ≤ T z₁` then every `z` with `T z₁ ≤ T z` also has
`T z₀ ≤ T z`, so the event defining `pValue D T z₁` is contained in the one
defining `pValue D T z₀`, and `Design.probEvent_mono` compares them.  The
containment is one `le_trans`. -/
theorem ex_pValue_antitone (D : Design n) (T : Assignment n → ℝ) {z₀ z₁ : Assignment n}
    (h : T z₀ ≤ T z₁) : pValue D T z₁ ≤ pValue D T z₀ := by
  refine D.probEvent_mono _ _ fun z hz => ?_
  exact le_trans h hz

Validity, as we are about to see, is a property of the randomization and not of the statistic. It is easy to under-read that as “the choice of statistic does not matter”. It matters entirely — for power. A statistic that ignores the data is perfectly valid and perfectly useless:

/-- Exercise (warm-up).  **A statistic that ignores the data never rejects.**  If `T`
is constant then the p-value is `1`, whatever the design and whatever was observed.

The event `{z : T z₀ ≤ T z}` is `{z : c ≤ c}`, which is all of them.
`Design.probEvent_congr` replaces an event by a logically equivalent one — here by
`fun _ => True` — and `Design.probEvent_true` says the sure event has probability
one.  (Validity is not the same as usefulness: this test is valid at every `α < 1`
because it never rejects at all.) -/
theorem ex_pValue_const (D : Design n) (c : ℝ) (z₀ : Assignment n) :
    pValue D (fun _ => c) z₀ = 1 := by
  sorry

/-- Exercise (warm-up).  **A statistic that ignores the data never rejects.**  If `T`
is constant then the p-value is `1`, whatever the design and whatever was observed.

The event `{z : T z₀ ≤ T z}` is `{z : c ≤ c}`, which is all of them.
`Design.probEvent_congr` replaces an event by a logically equivalent one — here by
`fun _ => True` — and `Design.probEvent_true` says the sure event has probability
one.  (Validity is not the same as usefulness: this test is valid at every `α < 1`
because it never rejects at all.) -/
theorem ex_pValue_const (D : Design n) (c : ℝ) (z₀ : Assignment n) :
    pValue D (fun _ => c) z₀ = 1 := by
  have hcongr : ∀ z : Assignment n, ((fun _ => c) z₀ ≤ (fun _ => c) z) ↔ True := by
    intro z
    simp
  rw [pValue, D.probEvent_congr hcongr, D.probEvent_true]

/-- Exercise (core).  **A one-sided test never gives a larger p-value than the
two-sided test built from the same statistic**, provided the observed value points
in the direction the one-sided test is looking (`0 ≤ T z₀`).

If `T z₀ ≤ T z` then `|T z₀| = T z₀ ≤ T z ≤ |T z|`, so the one-sided event is
contained in the two-sided one and `Design.probEvent_mono` applies.  The two facts
about absolute values you need are `abs_of_nonneg` and `le_abs_self`; a `calc`
chain reads best here. -/
theorem ex_pValue_one_sided_le (D : Design n) (T : Assignment n → ℝ) (z₀ : Assignment n)
    (h₀ : 0 ≤ T z₀) : pValue D T z₀ ≤ pValue D (fun z => |T z|) z₀ := by
  sorry

/-- Exercise (core).  **A one-sided test never gives a larger p-value than the
two-sided test built from the same statistic**, provided the observed value points
in the direction the one-sided test is looking (`0 ≤ T z₀`).

If `T z₀ ≤ T z` then `|T z₀| = T z₀ ≤ T z ≤ |T z|`, so the one-sided event is
contained in the two-sided one and `Design.probEvent_mono` applies.  The two facts
about absolute values you need are `abs_of_nonneg` and `le_abs_self`; a `calc`
chain reads best here. -/
theorem ex_pValue_one_sided_le (D : Design n) (T : Assignment n → ℝ) (z₀ : Assignment n)
    (h₀ : 0 ≤ T z₀) : pValue D T z₀ ≤ pValue D (fun z => |T z|) z₀ := by
  refine D.probEvent_mono _ _ fun z hz => ?_
  calc |T z₀| = T z₀ := abs_of_nonneg h₀
    _ ≤ T z := hz
    _ ≤ |T z| := le_abs_self (T z)

Validity

Informally: for every design, every statistic and every level $\alpha \ge 0$,

/-- Exercise (stretch).  **Validity of the Fisher randomization test.**

For every design `D`, every test statistic `T`, and every level `α ≥ 0`,

`P( pValue D T z ≤ α ) ≤ α`.

Read statistically: the rule "reject the sharp null when the randomization p-value is at most
`α`" has type-I error at most `α`.  There is no model, no asymptotic approximation, no
condition on `T`, and no condition on the design beyond its being a design.  The guarantee is
exact and finite-sample.  (The inequality can be strict — with only six possible assignments a
test at `α = 0.2` really rejects with probability `1/6`, because p-values live on a coarse
grid.)

The proof is a neat finite argument.  Let `S = {z : pValue D T z ≤ α}` be the rejection
region.  If `S` is empty its probability is `0 ≤ α`.  Otherwise pick `z*` in `S` minimising
`T` over `S` (`Finset.exists_min_image`, which needs `S` finite — it is).  Every `z ∈ S` has
`T z* ≤ T z`, so `S` is contained in the event `{T z* ≤ T ·}` whose probability is by
definition `pValue D T z*`, and that is `≤ α` because `z* ∈ S`. -/
theorem pValue_valid {α : ℝ} (hα : 0 ≤ α) :
    D.probEvent (fun z => pValue D T z ≤ α) ≤ α := by
  sorry

/-- Exercise (stretch).  **Validity of the Fisher randomization test.**

For every design `D`, every test statistic `T`, and every level `α ≥ 0`,

`P( pValue D T z ≤ α ) ≤ α`.

Read statistically: the rule "reject the sharp null when the randomization p-value is at most
`α`" has type-I error at most `α`.  There is no model, no asymptotic approximation, no
condition on `T`, and no condition on the design beyond its being a design.  The guarantee is
exact and finite-sample.  (The inequality can be strict — with only six possible assignments a
test at `α = 0.2` really rejects with probability `1/6`, because p-values live on a coarse
grid.)

The proof is a neat finite argument.  Let `S = {z : pValue D T z ≤ α}` be the rejection
region.  If `S` is empty its probability is `0 ≤ α`.  Otherwise pick `z*` in `S` minimising
`T` over `S` (`Finset.exists_min_image`, which needs `S` finite — it is).  Every `z ∈ S` has
`T z* ≤ T z`, so `S` is contained in the event `{T z* ≤ T ·}` whose probability is by
definition `pValue D T z*`, and that is `≤ α` because `z* ∈ S`. -/
theorem pValue_valid {α : ℝ} (hα : 0 ≤ α) :
    D.probEvent (fun z => pValue D T z ≤ α) ≤ α := by
  rcases Finset.eq_empty_or_nonempty (univ.filter (fun z => pValue D T z ≤ α)) with hemp | hne
  · -- nobody rejects: the probability is `0`
    have h0 : D.probEvent (fun z => pValue D T z ≤ α) = 0 := by
      rw [Design.probEvent, hemp, Finset.sum_empty]
    rw [h0]
    exact hα
  · -- somebody rejects: take the least extreme rejecting assignment `z*`
    obtain ⟨zs, hzs, hmin⟩ := Finset.exists_min_image _ T hne
    have hzs' : pValue D T zs ≤ α := (Finset.mem_filter.mp hzs).2
    have hsub : ∀ z, pValue D T z ≤ α → T zs ≤ T z := fun z hz =>
      hmin z (Finset.mem_filter.mpr ⟨Finset.mem_univ z, hz⟩)
    calc D.probEvent (fun z => pValue D T z ≤ α)
        ≤ D.probEvent (fun z => T zs ≤ T z) := D.probEvent_mono _ _ hsub
      _ = pValue D T zs := rfl
      _ ≤ α := hzs'

The proof deserves reading even if you do not write it. Everything hinges on the rejection region being upward closed in T: if z rejects and T z ≤ T z' then z' rejects, by pValue_antitone. An upward-closed set of assignments is contained in the tail above its own minimal element, and the probability of that tail is a p-value — the p-value at the minimal element, which rejects, so it is at most $\alpha$. The counting is the whole argument; nothing about distributions, and nothing about $T$, is used.

Notice which line does the work in Lean: _ = pValue D T zs := rfl. The probability of the event “at least as extreme as zs” is the p-value at zs by definition — rfl means the two sides are literally the same term after unfolding. When you are reading an autoformalization, a rfl step like this one is a place to slow down: it is where a definition is silently doing the job the statement claims a theorem is doing. Here it is honest, because pValue was defined as that probability three screens earlier.

Why 0 ≤ α is needed

The informal statement above is what I first wrote in the brief for this chapter, minus the hypothesis: for any $\alpha$, $Pr(p\le \alpha )\le \alpha$. It is false, and Lean will not let you prove it.

/-- **The unconditional form of the validity theorem is false.**  Drop `0 ≤ α` and
the statement claims, at `α = -1`, that a probability is at most `-1`.

Probabilities are nonnegative, and at a negative level the rejection region is
empty, so the left-hand side is exactly `0`: the "theorem" would say `0 ≤ -1`.
`probEvent_pValue_of_neg` is the library lemma recording that degenerate case, and
`norm_num` finishes off the arithmetic absurdity. -/
theorem not_pValue_valid_without_hypothesis :
    ¬ ∀ (D : Design 4) (T : Assignment 4 → ℝ) (α : ℝ),
        D.probEvent (fun z => pValue D T z ≤ α) ≤ α := by
  intro h
  have hle := h (completeRandomization 4 2 (by norm_num)) (absDiM exampleY) (-1)
  rw [probEvent_pValue_of_neg _ _ (by norm_num : (-1 : ℝ) < 0)] at hle
  norm_num at hle

At $\alpha =-1$ the rejection region is empty — p-values are nonnegative — so the left-hand side is exactly $0$, and the claim reduces to $0\le -1$. The mathematics is unharmed; nobody tests at a negative level. But the statement was wrong, and this is the most common way an autoformalization goes wrong: not by saying something absurd, but by being almost right, with a degenerate corner that the informal version never contemplated because a human reader silently supplies ”$\alpha$ is a probability”.

The honest headline

pValue_valid is stated for a fixed statistic T. That is not quite what an analyst does: they see the data, and then compute a statistic from it. The statistic is a function of the assignment through the observed outcomes, which is a different object, and the theorem has to say so.

theorem randomization_test_valid {P : Population n} (hnull : SharpNull P) (D : Design n)
    (S : (Fin n → ℝ) → Assignment n → ℝ) {α : ℝ} (hα : 0 ≤ α) :
    D.probEvent (fun z => pValue D (S (Yobs P z)) z ≤ α) ≤ α := by
  have hfixed : ∀ z : Assignment n, S (Yobs P z) = S P.y0 := fun z => by rw [hnull.Yobs_eq z]
  rw [D.probEvent_congr (B := fun z => pValue D (S P.y0) z ≤ α) fun z => by rw [hfixed z]]
  exact pValue_valid D (S P.y0) hα

Read the quantifiers. S : (Fin n → ℝ) → Assignment n → ℝ is the analyst’s rule, chosen in advance but arbitrary; Yobs P z is the data the assignment z would produce; pValue D (S (Yobs P z)) z is the p-value that analyst would report on that data. The event is over z, and z appears three times in it: in the data, in the statistic built from the data, and as the realized assignment whose extremeness is being judged. That triple occurrence is the thing a superficial formalization drops.

And it is exactly where the sharp null is spent. Under SharpNull P, Yobs P z does not depend on z, so S (Yobs P z) collapses to the single fixed function S P.y0 and pValue_valid applies verbatim. Under the weak null the collapse fails, the statistic genuinely varies with the assignment, and this argument gives nothing at all — which is why randomization tests of $\tau =0$ are a separate and much harder (and asymptotic) subject.

The instance a data analyst actually runs is randomization_test_absDiM_valid, which is randomization_test_valid with S := absDiM and no further proof.

The worked example

Four units, two treated, outcomes $y=(1,2,3,4)$, and the coin treated units 0 and 1. There are $(\genfrac{}{}{0ex}{}{4}{2})=6$ equally likely assignments; since every one of them treats two units and controls two, the statistic is the gap between two pair averages.

The observed value is $2$; two of the six assignments reach it; the p-value is $2/6=1/3$. The library proves exactly that from the definitions; every number in the table above is derived rather than asserted.

/-- **The worked example.**  With outcomes `y = (1,2,3,4)` and the realised assignment
`{0,1}`, the two-sided randomization p-value is exactly `1/3`. -/
theorem example_pValue :
    pValue (completeRandomization 4 2 (by norm_num)) (absDiM exampleY) {0, 1} = 1 / 3 := by
  obtain ⟨v01, v02, v03, v12, v13, v23⟩ := exampleY_absDiM
  have hiff : ∀ z ∈ ((univ : Finset (Fin 4)).powersetCard 2),
      (absDiM exampleY {0, 1} ≤ absDiM exampleY z ↔ (z = {0, 1} ∨ z = {2, 3})) := by
    intro z hz
    rw [powersetCard_two_fin_four] at hz
    simp only [Finset.mem_insert, Finset.mem_singleton] at hz
    rcases hz with rfl | rfl | rfl | rfl | rfl | rfl
    all_goals simp only [v01, v02, v03, v12, v13, v23]
    all_goals norm_num
    all_goals decide
  have hcard : #({x ∈ ((univ : Finset (Fin 4)).powersetCard 2) |
      x = ({0, 1} : Finset (Fin 4)) ∨ x = ({2, 3} : Finset (Fin 4))}) = 2 := by decide
  rw [pValue, completeRandomization_probEvent, Finset.filter_congr hiff, hcard,
    show Nat.choose 4 2 = 6 from by decide]
  norm_num

Three tactic moves make that proof go through, and all three are worth having in your vocabulary.

decide evaluates a decidable proposition and checks the answer is True. Membership, subsets and cardinalities of Finset (Fin 4) are all computable, so the enumeration of the design’s support and the count of the rejection region are both one word:

/-- The support of the design: the six two-element subsets of `Fin 4`, found by
brute force.  `decide` evaluates the `Finset` on both sides and compares them; at
this size that is instantaneous. -/
example : ((univ : Finset (Fin 4)).powersetCard 2)
    = ({{0, 1}, {0, 2}, {0, 3}, {1, 2}, {1, 3}, {2, 3}} : Finset (Finset (Fin 4))) := by
  decide

-- The same six assignments, printed rather than proved.
#eval ((univ : Finset (Fin 4)).powersetCard 2)

/-- How many of them are at least as extreme as the observed one, once the
real-valued comparison has been replaced by a decidable description of the same
set.  This is the count that becomes the numerator of the p-value. -/
example : #({z ∈ (univ : Finset (Fin 4)).powersetCard 2 |
    z = ({0, 1} : Finset (Fin 4)) ∨ z = ({2, 3} : Finset (Fin 4))}) = 2 := by
  decide

/-- And the denominator. -/
example : Nat.choose 4 2 = 6 := by decide

Finset.filter_congr replaces the filtering predicate by an equivalent one, given that the two agree on the set being filtered. This is the step that lets decide run at all: the real inequality absDiM exampleY {0,1} ≤ absDiM exampleY z is not decidable — Lean’s instance for it is classical, and computes nothing — so the proof first shows that on the six assignments of the design it is equivalent to z = {0,1} ∨ z = {2,3}, which is decidable. Real analysis on the left, finite combinatorics on the right, and the bridge is hiff.

completeRandomization_probEvent turns the abstract probability into $\mathrm{\#}\{\text{rejecting assignments}\}/(\genfrac{}{}{0ex}{}{n}{n_1})$, which is the calculation a reader would do by hand.

/-- Exercise (stretch).  **The same data, a different realised assignment.**  Had the
coin treated units `0` and `2` instead, the observed statistic would have been `1`
rather than `2`, and four of the six assignments reach `1`.  So the p-value is
`4/6 = 2/3`.

The proof is the proof of `example_pValue` with a longer list.  `exampleY_absDiM`
supplies the six values; `hiff` replaces the real inequality
`absDiM exampleY {0,2} ≤ absDiM exampleY z` by an explicit list of assignments,
which is decidable; `Finset.filter_congr` swaps one predicate for the other inside
the count; and `decide` does the counting. -/
theorem ex_example_pValue_02 :
    pValue (completeRandomization 4 2 (by norm_num)) (absDiM exampleY) {0, 2} = 2 / 3 := by
  sorry

/-- Exercise (stretch).  **The same data, a different realised assignment.**  Had the
coin treated units `0` and `2` instead, the observed statistic would have been `1`
rather than `2`, and four of the six assignments reach `1`.  So the p-value is
`4/6 = 2/3`.

The proof is the proof of `example_pValue` with a longer list.  `exampleY_absDiM`
supplies the six values; `hiff` replaces the real inequality
`absDiM exampleY {0,2} ≤ absDiM exampleY z` by an explicit list of assignments,
which is decidable; `Finset.filter_congr` swaps one predicate for the other inside
the count; and `decide` does the counting. -/
theorem ex_example_pValue_02 :
    pValue (completeRandomization 4 2 (by norm_num)) (absDiM exampleY) {0, 2} = 2 / 3 := by
  obtain ⟨v01, v02, v03, v12, v13, v23⟩ := exampleY_absDiM
  have hiff : ∀ z ∈ ((univ : Finset (Fin 4)).powersetCard 2),
      (absDiM exampleY {0, 2} ≤ absDiM exampleY z ↔
        (z = {0, 1} ∨ z = {0, 2} ∨ z = {1, 3} ∨ z = {2, 3})) := by
    intro z hz
    rw [powersetCard_two_fin_four] at hz
    simp only [Finset.mem_insert, Finset.mem_singleton] at hz
    rcases hz with rfl | rfl | rfl | rfl | rfl | rfl
    all_goals simp only [v01, v02, v03, v12, v13, v23]
    all_goals norm_num
    all_goals decide
  have hcard : #({z ∈ ((univ : Finset (Fin 4)).powersetCard 2) |
      z = ({0, 1} : Finset (Fin 4)) ∨ z = ({0, 2} : Finset (Fin 4)) ∨
        z = ({1, 3} : Finset (Fin 4)) ∨ z = ({2, 3} : Finset (Fin 4))}) = 4 := by decide
  rw [pValue, completeRandomization_probEvent, Finset.filter_congr hiff, hcard,
    show Nat.choose 4 2 = 6 from by decide]
  norm_num

Change the data

The p-value is a function of the data as much as of the design. Same design, same realized assignment, a different outcome vector:

/-- **Change these four numbers and re-run.**  The same design, the same realised
assignment `{0,1}`, different data.

With `y = (5,1,2,4)` the treated pair `{0,1}` has mean `3` and the control pair
`{2,3}` has mean `3`, so the observed statistic is `0` — the *least* extreme value
the statistic can take.  Every assignment is then at least as extreme, and the
p-value is exactly `1`. -/
noncomputable def myY : Fin 4 → ℝ := ![5, 1, 2, 4]

/-- Exercise (stretch).  With `myY` the observed difference in means is zero, so the
randomization p-value is `1`.

No enumeration is needed for this one: `absDiM` is an absolute value, hence
nonnegative, so once you know the observed value is `0` the event "at least as
extreme as observed" is the sure event.  `absDiM_pair` computes the observed value;
`Design.probEvent_congr` and `Design.probEvent_true` finish. -/
theorem ex_myY_pValue :
    pValue (completeRandomization 4 2 (by norm_num)) (absDiM myY) {0, 1} = 1 := by
  sorry

/-- **Change these four numbers and re-run.**  The same design, the same realised
assignment `{0,1}`, different data.

With `y = (5,1,2,4)` the treated pair `{0,1}` has mean `3` and the control pair
`{2,3}` has mean `3`, so the observed statistic is `0` — the *least* extreme value
the statistic can take.  Every assignment is then at least as extreme, and the
p-value is exactly `1`. -/
noncomputable def myY : Fin 4 → ℝ := ![5, 1, 2, 4]

/-- Exercise (stretch).  With `myY` the observed difference in means is zero, so the
randomization p-value is `1`.

No enumeration is needed for this one: `absDiM` is an absolute value, hence
nonnegative, so once you know the observed value is `0` the event "at least as
extreme as observed" is the sure event.  `absDiM_pair` computes the observed value;
`Design.probEvent_congr` and `Design.probEvent_true` finish. -/
theorem ex_myY_pValue :
    pValue (completeRandomization 4 2 (by norm_num)) (absDiM myY) {0, 1} = 1 := by
  have h0 : absDiM myY {0, 1} = 0 := by
    rw [absDiM_pair myY 0 1 2 3 (by decide) (by decide) (by decide)]
    norm_num [myY, Matrix.cons_val_two, Matrix.cons_val_three, Matrix.tail_cons,
      Matrix.head_cons]
  have hcongr : ∀ z : Assignment 4, (absDiM myY {0, 1} ≤ absDiM myY z) ↔ True := by
    intro z
    rw [h0]
    simp only [absDiM, iff_true]
    exact abs_nonneg _
  rw [pValue, Design.probEvent_congr _ hcongr, Design.probEvent_true]

If you change myY to something else, that proof will break, and the break is informative: the h0 line is the only place the data enters, and Lean will tell you the observed statistic is no longer 0. To explore rather than to prove, run the calculator instead.

A calculator you can trust

pValue is real-valued, hence noncomputable: Lean will not run it. The same number can be computed over the integers, though, because every assignment in this design splits four units two and two, so comparing $|\mathrm{D}\mathrm{i}\mathrm{M}|$ across assignments is comparing the integer gap “treated total minus control total”.

/-- Treated total minus control total, over the integers.  Twice the difference in
means, for a design in which every assignment splits four units two and two. -/
def gap (y : Fin 4 → ℤ) (z : Finset (Fin 4)) : ℤ := (∑ i ∈ z, y i) - (∑ i ∈ zᶜ, y i)

/-- How many of the six equally likely assignments are at least as extreme as `z₀`.
Integers are decidable, so this is a number Lean can compute. -/
def rejectCount (y : Fin 4 → ℤ) (z₀ : Finset (Fin 4)) : ℕ :=
  #{z ∈ (univ : Finset (Fin 4)).powersetCard 2 | |gap y z₀| ≤ |gap y z|}

/-- The two-sided randomization p-value of the worked design, as a rational number
you can look at. -/
def pValueCalc (y : Fin 4 → ℤ) (z₀ : Finset (Fin 4)) : ℚ := (rejectCount y z₀ : ℚ) / 6

-- The randomization distribution of the worked example, as integer gaps:
-- `4, 2, 0, 0, 2, 4`, i.e. twice `2, 1, 0, 0, 1, 2`.  (`.val` keeps the
-- multiplicities; `Finset.image` would silently collapse the two zeros.)
#eval (((univ : Finset (Fin 4)).powersetCard 2).val.map fun z => |gap ![1, 2, 3, 4] z|)

-- Change the outcome vector, change the realised assignment, re-run.
#eval pValueCalc ![1, 2, 3, 4] {0, 1}   -- 1/3, the worked example
#eval pValueCalc ![1, 2, 3, 4] {0, 2}   -- 2/3
#eval pValueCalc ![1, 2, 3, 4] {0, 3}   -- 1
#eval pValueCalc ![5, 1, 2, 4] {0, 1}   -- 1: the observed split is the least extreme

#eval is not a proof. A computation that prints 1/3 says nothing about pValue unless something connects the two, and the connection is a theorem like any other:

/-- With four units and two treated, the statistic is the integer gap over two. -/
theorem absDiM_intCast_eq (y : Fin 4 → ℤ) {z : Finset (Fin 4)} (hz : #z = 2) :
    absDiM (fun i => (y i : ℝ)) z = |(gap y z : ℝ)| / 2 := by
  have hc : #zᶜ = 2 := by
    rw [Finset.card_compl, hz, Fintype.card_fin]
  simp only [absDiM, DiM, imputeSharp_y1, imputeSharp_y0, gap, hz, hc]
  push_cast
  rw [← sub_div, abs_div]
  norm_num

/-- **The calculator is correct.**  For any integer outcome vector and any realised
assignment of the right size, the rational number `pValueCalc` computes is the real
number `pValue` defines.

So `#eval pValueCalc ![1,2,3,4] {0,1}` printing `1/3` is not a coincidence beside
`example_pValue`; it is the same theorem, evaluated. -/
theorem pValue_eq_pValueCalc (y : Fin 4 → ℤ) {z₀ : Finset (Fin 4)} (hz₀ : #z₀ = 2) :
    pValue (completeRandomization 4 2 (by norm_num)) (absDiM fun i => (y i : ℝ)) z₀
      = (pValueCalc y z₀ : ℝ) := by
  have hiff : ∀ z ∈ ((univ : Finset (Fin 4)).powersetCard 2),
      (absDiM (fun i => (y i : ℝ)) z₀ ≤ absDiM (fun i => (y i : ℝ)) z
        ↔ |gap y z₀| ≤ |gap y z|) := by
    intro z hz
    have hz2 : #z = 2 := (Finset.mem_powersetCard.mp hz).2
    rw [absDiM_intCast_eq y hz₀, absDiM_intCast_eq y hz2,
      div_le_div_iff_of_pos_right (by norm_num : (0:ℝ) < 2), ← Int.cast_abs, ← Int.cast_abs,
      Int.cast_le]
  rw [pValue, completeRandomization_probEvent, Finset.filter_congr hiff]
  simp only [pValueCalc, rejectCount, Rat.cast_div, Rat.cast_natCast, Rat.cast_ofNat]
  rw [show Nat.choose 4 2 = 6 from by decide]
  push_cast
  ring

Now the #eval output is a statement about pValue, evaluated. Without pValue_eq_pValueCalc it would be a statement about pValueCalc and nothing else — which is a general lesson about auto-formalized code that mixes definitions with computations: ask what has been proved to agree with what.

Reading autoformalized Lean

Two plausible-looking formalizations of this chapter’s theorem, both wrong.

The level is unconstrained

-- WRONG: no hypothesis on α
theorem pValue_valid' (D : Design n) (T : Assignment n → ℝ) (α : ℝ) :
    D.probEvent (fun z => pValue D T z ≤ α) ≤ α := by
  sorry

This is the version I wrote in English first, and it is false at every negative $\alpha$ — refuted above by not_pValue_valid_without_hypothesis. The tell is structural rather than mathematical: a real parameter with no constraint on it, in a statement whose informal reading assumes it is a probability. When you audit an autoformalized theorem, list the free variables and ask, for each, what the English took for granted about its range.

The observed assignment does not count itself

-- WRONG: strict inequality
noncomputable def pValue' (D : Design n) (T : Assignment n → ℝ) (z₀ : Assignment n) : ℝ :=
  D.probEvent (fun z => T z₀ < T z)

theorem pValue'_valid (D : Design n) (T : Assignment n → ℝ) {α : ℝ} (hα : 0 ≤ α) :
    D.probEvent (fun z => pValue' D T z ≤ α) ≤ α := by
  sorry

One character. pValue' counts the assignments strictly more extreme than the observed one, so it can be zero — take the worked example, observe {0,1}, and pValue' = 0 because nothing beats a statistic of 2. Then the rejection region at $\alpha =0$ contains both {0,1} and {2,3}, and the type-I error is . Validity fails.

The flaw is invisible in the informal phrase “the probability of a result as extreme as the one observed”, which is exactly why it is worth checking against a hand-computable example: at $n=4$ the whole distribution fits in the table above, and a formalization that gets the boundary wrong is caught in seconds. pValue_pos is the library’s insurance, and its proof — the observed assignment satisfies T z₀ ≤ T z₀ — is the same character.