Fisher randomization tests

This chapter studies Fisher randomization tests. Under the sharp null of no effect for any unit, the missing potential outcomes can be imputed exactly. The design then supplies the reference distribution of any chosen test statistic. Rejecting when the randomization p-value is at most $\alpha$ has type-I error at most $\alpha$ in finite samples. The proof requires no asymptotic approximation, but it depends on the sharp null rather than only an average-effect null.¹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.

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 chapter 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.

theorem ex_probEvent_nonneg (D : Design n) (A : Assignment n → Prop) [DecidablePred A] :
    0 ≤ D.probEvent A := by
  sorry

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

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

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.

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. The sharp null constrains $n$ pairs of numbers; $\tau =0$ constrains one average.

Lemma. If $y_1(i)=y_0(i)$ for every unit $i$, then $\tau =0$. The converse fails.

This is SharpNull.tau_eq_zero in the library, and its Lean proof is those two sentences in order: a have hzero : ∀ i, P.effect i = 0 obtained from sub_eq_zero_of_eq (h i), then a simp only that unfolds Population.tau and mean, rewrites every summand to 0 with hzero, collapses the sum with Finset.sum_const_zero and divides with zero_div.

What the sharp null buys is that the missing half of the table is not missing.

Lemma. Under the sharp null the observed outcome vector does not depend on the assignment: $Y^{\mathrm{o}\mathrm{b}\mathrm{s}}(z)=y_0$ for every $z$.

So the experimenter can compute what would have been observed under every assignment the design could have produced, not only the one that happened. That is the whole of what the sharp null is for.

theorem ex_SharpNull_Yobs_eq {P : Population n} (h : SharpNull P) (z : Assignment n) :
    Yobs P z = P.y0 := by
  sorry

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:

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.

theorem ex_imputeSharp_Yobs (y : Fin n → ℝ) (z : Assignment n) :
    Yobs (imputeSharp y) z = y := by
  sorry

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.

Lemma. Under the sharp null, the population imputed from the observed data equals the true population, whichever assignment was drawn.

theorem ex_imputeSharp_exact {P : Population n} (h : SharpNull P) (z₀ : Assignment n) :
    imputeSharp (Yobs P z₀) = P := by
  sorry

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.

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:

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.

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

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.

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

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:

theorem ex_pValue_const (D : Design n) (c : ℝ) (z₀ : Assignment n) :
    pValue D (fun _ => c) z₀ = 1 := by
  sorry

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]

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

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

Write $\mathcal{D}$ for the design and $p(z)=\underset{\mathcal{D}}{Pr}(T(Z^{\prime })\ge T(z))$ for the p-value at the assignment $z$, the probability that a fresh draw $Z^{\prime }$ from the design is at least as extreme as $z$. The claim is that this quantity is a valid p-value in the exact, finite-sample sense.

Proposition (validity). Let $\mathcal{D}$ be any distribution on the assignments, $T$ any statistic, and $\alpha \ge 0$ any level. Then

Nothing in that argument mentions a distribution of $T$, a limit, or a model. It is a counting fact about an upward-closed set: if $z$ rejects and $T(z)\le T(z^{\prime })$ then $z^{\prime }$ rejects too, so the rejection region is contained in the tail above its own minimal element, and the probability of that tail is a p-value which rejects.

theorem pValue_valid {α : ℝ} (hα : 0 ≤ α) :
    D.probEvent (fun z => pValue D T z ≤ α) ≤ α := by
  sorry

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
  ·
    have h0 : D.probEvent (fun z => pValue D T z ≤ α) = 0 := by
      rw [Design.probEvent, hemp, Finset.sum_empty]
    rw [h0]
    exact hα
  ·
    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'

One thing to note about the Lean, since the hint above already names its lemmas: the upward closure never becomes a lemma of its own there. pValue_antitone, the exercise two pages back, is that property stated in the abstract, and pValue_valid does not quote it — minimality of zs hands the containment over directly, as hsub.

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.

Claim. The unconstrained statement — for every design, statistic and level $\alpha$, $Pr(p\le \alpha )\le \alpha$ — is false.

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

The Lean proof is that paragraph, in four lines.

1. intro h

   The goal is a negation, ¬ Q, which in Lean is Q → False. So intro h assumes the bad statement and leaves False to derive: the proof will be a refutation by counterexample.

   Goal

   h : ∀ (D : Design 4) (T : Assignment 4 → ℝ) (α : ℝ), (D.probEvent fun z => pValue D T z ≤ α) ≤ α
   ⊢ False

2. have hle := h (completeRandomization 4 2 (by norm_num)) (absDiM exampleY) (-1)

   Instantiate the universal claim at one design, one statistic and one level. Nothing is being proved yet — h was assumed — but a ∀ statement is only as true as its worst instance, and this is the instance that will fail. The by norm_num fills the 2 ≤ 4 argument that completeRandomization demands. (The context still holds h; it is elided here and below.)

   Goal

   hle :
     ((completeRandomization 4 2 ⋯).probEvent fun z => pValue (completeRandomization 4 2 ⋯) (absDiM exampleY) z ≤ -1) ≤ -1
   ⊢ False

3. rw [probEvent_pValue_of_neg _ _ (by norm_num : (-1 : ℝ) < 0)] at hle

   probEvent_pValue_of_neg is the library lemma saying that at a negative level the rejection region is empty, so its probability is 0. Rewriting with it inside the hypothesis replaces the left-hand side of hle by 0, and hle now reads 0 ≤ -1.

   Goal

   hle : 0 ≤ -1
   ⊢ False

4. norm_num at hle

   norm_num evaluates the numeric claim 0 ≤ -1, finds it false, and closes the goal: a hypothesis that normalizes to False proves anything, False included.

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α

1. have hfixed : ∀ z : Assignment n, S (Yobs P z) = S P.y0 := fun z => by rw [hnull.Yobs_eq z]

   The one place the sharp null is spent. SharpNull.Yobs_eq says the observed vector is P.y0 for every assignment, so rewriting with it turns S (Yobs P z) into S P.y0: the analyst’s statistic is a single fixed function, not a family indexed by z. This is the first sentence of the proof above.

2. rw [D.probEvent_congr (B := fun z => pValue D (S P.y0) z ≤ α) fun z => by rw [hfixed z]]

   Design.probEvent_congr says logically equivalent events have the same probability. The two events here are “the p-value of S (Yobs P z) at z is at most α” and “the p-value of the fixed S P.y0 at z is at most α”; the inner by rw [hfixed z] proves them equivalent at each z. Rewriting with the resulting equation replaces the goal by the fixed-statistic version.

   Goal

   hfixed : ∀ (z : Assignment n), S (Yobs P z) = S P.y0
   ⊢ (D.probEvent fun z => pValue D (S P.y0) z ≤ α) ≤ α

3. exact pValue_valid D (S P.y0) hα

   Which is pValue_valid verbatim, at the statistic S P.y0. The theorem is a corollary of the proposition and one rewrite; all the mathematics was in the counting argument.

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.

Proposition. For these data, this design, and the two-sided statistic $T(z)=|\mathrm{D}\mathrm{i}\mathrm{M}(z)|$, the randomization p-value of the realized assignment is $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

The Lean is that proof with the enumeration done by the kernel rather than by the reader. The only real work is the middle step: a decidable description of the set “at least as extreme as observed”.

1. obtain ⟨v01, v02, v03, v12, v13, v23⟩ := exampleY_absDiM

   exampleY_absDiM is the six-way conjunction proved earlier — the table, as a theorem. obtain splits it into six named equations v01 : absDiM exampleY {0, 1} = 2 and so on, which are the numbers the rest of the proof reasons with.

2. intro z hz

   The have hiff opens here. Its statement is the bridge of the whole proof: for every assignment z in the design’s support, “as extreme as the observed value” is equivalent to the decidable disjunction z = {0, 1} ∨ z = {2, 3}. intro takes an arbitrary such z together with its membership proof.

   Goal

   v01 : absDiM exampleY {0, 1} = 2
   v02 : absDiM exampleY {0, 2} = 1
   v03 : absDiM exampleY {0, 3} = 0
   v12 : absDiM exampleY {1, 2} = 0
   v13 : absDiM exampleY {1, 3} = 1
   v23 : absDiM exampleY {2, 3} = 2
   z : Finset (Fin 4)
   hz : z ∈ powersetCard 2 univ
   ⊢ absDiM exampleY {0, 1} ≤ absDiM exampleY z ↔ z = {0, 1} ∨ z = {2, 3}

3. rw [powersetCard_two_fin_four] at hz

   powersetCard_two_fin_four is the by decide enumeration of the support: the six two-element subsets of Fin 4, written out as a literal Finset. Rewriting hz with it turns “z is a two-element subset” into “z belongs to this list of six”.

4. simp only [Finset.mem_insert, Finset.mem_singleton] at hz

   Membership in a literal Finset is a chain of inserts, and these two lemmas unfold it: Finset.mem_insert turns z ∈ insert a s into z = a ∨ z ∈ s, and Finset.mem_singleton turns z ∈ {a} into z = a. What was a membership is now a six-way disjunction of equations.

   Goal

   hz : z = {0, 1} ∨ z = {0, 2} ∨ z = {0, 3} ∨ z = {1, 2} ∨ z = {1, 3} ∨ z = {2, 3}
   ⊢ absDiM exampleY {0, 1} ≤ absDiM exampleY z ↔ z = {0, 1} ∨ z = {2, 3}

5. rcases hz with rfl | rfl | rfl | rfl | rfl | rfl

   Case-split on that disjunction. Writing rfl in each branch rather than a name tells Lean to substitute: in the first branch z is replaced by {0, 1} everywhere, in the second by {0, 2}, and so on. Six goals, each about a concrete assignment.

6. all_goals simp only [v01, v02, v03, v12, v13, v23]

   Rewrite with the table in all six branches at once. Both sides of the inequality become numerals — this is where real analysis stops and arithmetic begins.

   Goal

   case inl
   ⊢ 2 ≤ 2 ↔ True ∨ {0, 1} = {2, 3}
   
   case inr.inl
   ⊢ 2 ≤ 1 ↔ {0, 2} = {0, 1} ∨ {0, 2} = {2, 3}
   
   case inr.inr.inl
   ⊢ 2 ≤ 0 ↔ {0, 3} = {0, 1} ∨ {0, 3} = {2, 3}
   
   case inr.inr.inr.inl
   ⊢ 2 ≤ 0 ↔ {1, 2} = {0, 1} ∨ {1, 2} = {2, 3}
   
   case inr.inr.inr.inr.inl
   ⊢ 2 ≤ 1 ↔ {1, 3} = {0, 1} ∨ {1, 3} = {2, 3}
   
   case inr.inr.inr.inr.inr
   ⊢ 2 ≤ 2 ↔ {2, 3} = {0, 1} ∨ True

7. all_goals norm_num

   norm_num decides each numeric inequality and simplifies the iff around it. Two branches close outright (2 ≤ 2 against a disjunction containing True); the four where the observed value is not attained leave the residue “this assignment is neither {0, 1} nor {2, 3}” — a statement about Finset equality, which norm_num has no opinion about.

   Goal

   case inr.inl
   ⊢ ¬{0, 2} = {0, 1} ∧ ¬{0, 2} = {2, 3}
   
   case inr.inr.inl
   ⊢ ¬{0, 3} = {0, 1} ∧ ¬{0, 3} = {2, 3}
   
   case inr.inr.inr.inl
   ⊢ ¬{1, 2} = {0, 1} ∧ ¬{1, 2} = {2, 3}
   
   case inr.inr.inr.inr.inl
   ⊢ ¬{1, 3} = {0, 1} ∧ ¬{1, 3} = {2, 3}

8. all_goals decide

   Equality of literal Finsets over Fin 4 is decidable, so decide evaluates the four residual goals to True and closes them. hiff is proved.

9. have hcard : #({x ∈ ((univ : Finset (Fin 4)).powersetCard 2) | x = ({0, 1} : Finset (Fin 4)) ∨ x = ({2, 3} : Finset (Fin 4))}) = 2 := by decide

   The count. Now that the filtering predicate is decidable, the size of the rejection region is a computation, and decide runs it: two assignments out of the six.

10. rw [pValue, completeRandomization_probEvent, Finset.filter_congr hiff, hcard, show Nat.choose 4 2 = 6 from by decide]

    Five rewrites, left to right. pValue unfolds the definition to a probEvent; completeRandomization_probEvent turns that probability into “count divided by $(\genfrac{}{}{0ex}{}{4}{2})$”; Finset.filter_congr hiff swaps the undecidable filtering predicate for the decidable one, which is legitimate because hiff holds on the set being filtered; hcard replaces the count by 2; and the show term computes $(\genfrac{}{}{0ex}{}{4}{2})=6$.

    Goal

    ⊢ ↑2 * (↑6)⁻¹ = 1 / 3

11. norm_num

    Arithmetic in ℝ: $2\cdot 6^{-1}=1/3$.

Three of those moves are worth having in your vocabulary on their own.

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:

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

#eval ((univ : Finset (Fin 4)).powersetCard 2)

example : #({z ∈ (univ : Finset (Fin 4)).powersetCard 2 |
    z = ({0, 1} : Finset (Fin 4)) ∨ z = ({2, 3} : Finset (Fin 4))}) = 2 := by
  decide

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 $|\{\text{rejecting assignments}\}|/(\genfrac{}{}{0ex}{}{n}{n_1})$, which is the calculation a reader would do by hand.

theorem ex_example_pValue_02 :
    pValue (completeRandomization 4 2 (by norm_num)) (absDiM exampleY) {0, 2} = 2 / 3 := by
  sorry

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:

noncomputable def myY : Fin 4 → ℝ := ![5, 1, 2, 4]

theorem ex_myY_pValue :
    pValue (completeRandomization 4 2 (by norm_num)) (absDiM myY) {0, 1} = 1 := by
  sorry

noncomputable def myY : Fin 4 → ℝ := ![5, 1, 2, 4]

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”.

def gap (y : Fin 4 → ℤ) (z : Finset (Fin 4)) : ℤ := (∑ i ∈ z, y i) - (∑ i ∈ zᶜ, y i)

def rejectCount (y : Fin 4 → ℤ) (z₀ : Finset (Fin 4)) : ℕ :=
  #{z ∈ (univ : Finset (Fin 4)).powersetCard 2 | |gap y z₀| ≤ |gap y z|}

def pValueCalc (y : Fin 4 → ℤ) (z₀ : Finset (Fin 4)) : ℚ := (rejectCount y z₀ : ℚ) / 6

#eval (((univ : Finset (Fin 4)).powersetCard 2).val.map fun z => |gap ![1, 2, 3, 4] z|)

#eval pValueCalc ![1, 2, 3, 4] {0, 1}
#eval pValueCalc ![1, 2, 3, 4] {0, 2}
#eval pValueCalc ![1, 2, 3, 4] {0, 3}
#eval pValueCalc ![5, 1, 2, 4] {0, 1}

#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. Two claims are needed, and the first says the integer calculator computes the right statistic up to a known factor.

Lemma. Let $y\in {\mathbb{Z}}^4$ and let $z$ treat exactly two of the four units. Write $g(z)=\sum _{i\in z}{y}_i-\sum _{i\notin z}{y}_i$ for the integer gap between the treated and control totals. Then $|\mathrm{D}\mathrm{i}\mathrm{M}(z)|=|g(z)|/2$.

Proposition. With four units, two treated, the randomization p-value of the two-sided difference in means at any two-element $z_0$ equals the rational number the calculator returns.

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

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

The first proof is the one place in the chapter where casts between ℤ and ℝ have to be managed by hand.

1. have hc : #zᶜ = 2 := by rw [Finset.card_compl, hz, Fintype.card_fin]

   The control group also has two units, which the statement does not say and the arithmetic needs. Finset.card_compl gives #zᶜ = Fintype.card (Fin 4) - #z; hz replaces #z by 2 and Fintype.card_fin replaces the card of Fin 4 by 4, leaving 4 - 2 = 2, which rw closes by reflexivity.

2. simp only [absDiM, DiM, imputeSharp_y1, imputeSharp_y0, gap, hz, hc]

   Five unfoldings and two substitutions, in one pass. absDiM becomes |DiM …|; DiM becomes the difference of two group means; imputeSharp_y1 and imputeSharp_y0 replace both potential-outcome vectors of the imputed population by the data, which is how absDiM depends on y alone; gap unfolds the right-hand side into its two integer sums; and hz, hc put the literal 2 in both denominators.

   Goal

   hz : #z = 2
   hc : #zᶜ = 2
   ⊢ |(∑ i ∈ z, ↑(y i)) / ↑2 - (∑ i ∈ zᶜ, ↑(y i)) / ↑2| = |↑(∑ i ∈ z, y i - ∑ i ∈ zᶜ, y i)| / 2

3. push_cast

   The two sides live in ℝ but the right one is a cast of an integer subtraction. push_cast moves casts toward the leaves: ↑(∑ y i - ∑ y i) becomes ∑ ↑(y i) - ∑ ↑(y i), and the cast ↑2 in the denominators becomes the real numeral 2. Now both sides are real expressions in the same atoms.

   Goal

   ⊢ |(∑ i ∈ z, ↑(y i)) / 2 - (∑ i ∈ zᶜ, ↑(y i)) / 2| = |∑ x ∈ z, ↑(y x) - ∑ x ∈ zᶜ, ↑(y x)| / 2

4. rw [← sub_div, abs_div]

   Two rewrites that reshape the left-hand side to match the right. sub_div states (a - b) / c = a / c - b / c; used right-to-left it collects the two halves into one quotient. abs_div then states |a / b| = |a| / |b|, moving the modulus inside.

   Goal

   ⊢ |∑ i ∈ z, ↑(y i) - ∑ i ∈ zᶜ, ↑(y i)| / |2| = |∑ x ∈ z, ↑(y x) - ∑ x ∈ zᶜ, ↑(y x)| / 2

5. norm_num

   What is left is |…| / |2| = |…| / 2 with identical numerators, so all that remains is $|2|=2$; norm_num evaluates it.

The second proof reuses the machinery of the worked example, with hiff now proved once and for all rather than by enumeration.

1. have hz2 : #z = 2 := (Finset.mem_powersetCard.mp hz).2

   Inside hiff. Finset.mem_powersetCard says z ∈ s.powersetCard k is “z ⊆ s and #z = k”; the second component is the cardinality fact, which the lemma above needs in order to apply to z as well as to z₀.

2. rw [absDiM_intCast_eq y hz₀, absDiM_intCast_eq y hz2]

   Apply the lemma on both sides of the inequality: each |DiM| becomes the absolute integer gap over 2. Both sides now differ only by the cast and the division.

   Goal

   ⊢ |↑(gap y z₀)| / 2 ≤ |↑(gap y z)| / 2 ↔ |gap y z₀| ≤ |gap y z|

3. rw [div_le_div_iff_of_pos_right (by norm_num : (0:ℝ) < 2)]

   div_le_div_iff_of_pos_right states that for 0 < c, a / c ≤ b / c ↔ a ≤ b. Dividing both sides by the same positive constant is what the informal proof called “changes nothing”, and this is the lemma that says so.

   Goal

   ⊢ |↑(gap y z₀)| ≤ |↑(gap y z)| ↔ |gap y z₀| ≤ |gap y z|

4. rw [← Int.cast_abs, ← Int.cast_abs]

   Int.cast_abs states ↑|m| = |↑m|. Used right-to-left, twice, it moves each modulus inside the cast, so both sides of the real inequality are now casts of integer absolute values.

   Goal

   ⊢ ↑|gap y z₀| ≤ ↑|gap y z| ↔ |gap y z₀| ≤ |gap y z|

5. rw [Int.cast_le]

   Int.cast_le states (↑m : ℝ) ≤ ↑k ↔ m ≤ k: the cast into ℝ is an order embedding. Rewriting with it turns the goal into |gap y z₀| ≤ |gap y z| ↔ |gap y z₀| ≤ |gap y z|, which rw closes by reflexivity. hiff is proved, with no case analysis and no reference to the particular data.

6. rw [pValue, completeRandomization_probEvent, Finset.filter_congr hiff]

   The same three rewrites as in the worked example: unfold the p-value, turn the probability into a count over $(\genfrac{}{}{0ex}{}{4}{2})$, and swap the real-analytic filtering predicate for the integer one that hiff licenses. The left-hand side is now the calculator’s count, divided by a binomial coefficient.

   Goal

   ⊢ ↑(#({x ∈ powersetCard 2 univ | |gap y z₀| ≤ |gap y x|})) * (↑(Nat.choose 4 2))⁻¹ = ↑(pValueCalc y z₀)

7. simp only [pValueCalc, rejectCount, Rat.cast_div, Rat.cast_natCast, Rat.cast_ofNat]

   Unfold the right-hand side. pValueCalc is rejectCount / 6 in ℚ and rejectCount is the count itself; the three Rat.cast_* lemmas push the cast ℚ → ℝ through the division, through the natural-number count, and through the numeral 6, so that the right-hand side is a real expression in the same count.

   Goal

   ⊢ ↑(#({x ∈ powersetCard 2 univ | |gap y z₀| ≤ |gap y x|})) * (↑(Nat.choose 4 2))⁻¹ =
       ↑(#({x ∈ powersetCard 2 univ | |gap y z₀| ≤ |gap y x|})) / 6

8. rw [show Nat.choose 4 2 = 6 from by decide]

   Evaluate the binomial coefficient on the left, by kernel computation.

   Goal

   ⊢ ↑(#({x ∈ powersetCard 2 univ | |gap y z₀| ≤ |gap y x|})) * (↑6)⁻¹ =
       ↑(#({x ∈ powersetCard 2 univ | |gap y z₀| ≤ |gap y x|})) / 6

9. push_cast

   Normalize the remaining casts of the count so that the two sides are literally the same real atom multiplied and divided by numerals.

10. ring

    x * 6⁻¹ = x / 6 is a commutative-ring identity, and ring verifies it.

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

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

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.