第 3 章　完全随机化实验与 Fisher 随机化检验

潜在结果框架，与随机化实验有着天然的亲近。先把各类随机化实验里的因果推断想清楚，再去看更复杂的非实验研究，往往会轻松许多。

本书第二部分专门谈随机化实验。这一章从最简单的一种开始：完全随机化实验（completely randomized experiment, CRE）。

3.1　完全随机化实验（CRE）

设想一项实验有 \(n\) 个单元，其中 \(n_1\) 个接受处理，\(n_0\) 个接受对照。我们可以从处理分配机制来定义 CRE。1

定义 3.1（CRE）　固定 \(n_1\) 与 \(n_0\)，且 \(n=n_1+n_0\)。若处理分配机制为

\[ \operatorname{pr}(Z=z)=\binom{n}{n_1}^{-1}, \]

其中 \(z=(z_1,\ldots,z_n)\) 满足 \(\sum_{i=1}^{n}z_i=n_1\)、\(\sum_{i=1}^{n}(1-z_i)=n_0\)，则称这是一次完全随机化实验。

在定义 3.1 里，我们把处理下的潜在结果向量 \(Y(1)=(Y_1(1),\ldots,Y_n(1))\) 与对照下的 \(Y(0)=(Y_1(0),\ldots,Y_n(0))\) 都看作固定的。即便把它们当作随机，也可以先条件于它们；于是分配机制变成

\[ \operatorname{pr}\{Z=z\mid Y(1),Y(0)\}=\binom{n}{n_1}^{-1}, \]

因为在 CRE 中 \(Z \perp\!\!\!\perp \{Y(1),Y(0)\}\)。直观上，\(Z\) 就是把 \(n_1\) 个 \(1\) 与 \(n_0\) 个 \(0\) 随机排一排。

Fisher (1935) 在《实验设计》里指出，随机化有两桩好处：

  1. 它让处理组与对照组平均而言彼此可比；
  2. 它给统计推断提供了一个「有道理的根基」（reasoned basis）。

第一点比较好懂：随机分配不会故意偏向处理或对照。第二点更细——Fisher 的意思是：随机化本身，就足以支撑一种统计检验。这种检验，如今叫做 Fisher 随机化检验（Fisher Randomization Test, FRT）。本章就在 CRE 下，把 FRT 的基本想法慢慢讲清楚。

3.2　Fisher 随机化检验（FRT）

Fisher (1935) 关心的是下面这个原假设：2

\[ H_{0\mathrm{f}}:\quad Y_i(1)=Y_i(0)\quad\text{对所有单元 }i=1,\ldots,n. \]

Rubin (1980) 叫它尖锐原假设（sharp null）：在它底下，观测数据就能把所有潜在结果定死——\(Y(1)=Y(0)=Y=(Y_1,\ldots,Y_n)\)。它也常被叫作强原假设（例如 Wu and Ding, 2021）。

概念上，在 \(H_{0\mathrm{f}}\) 下，FRT 对任意检验统计量

\[ T=T(Z,Y) \tag{3.1} \]

都成立——只要它是观测数据的函数。在 \(H_{0\mathrm{f}}\) 下，\(Y\) 是固定的，\(T\) 里唯一的随机成分就是处理向量 \(Z\)。实验者决定 \(Z\) 的分布，也就决定了 \(T\) 在 \(H_{0\mathrm{f}}\) 下的分布。\(p\) 值，便从这里长出来。

在 CRE 中，\(Z\) 在集合 \(\{z_1,\ldots,z_M\}\) 上均匀，其中 \(M=\binom{n}{n_1}\)，每个 \(z_m\) 都是恰有 \(n_1\) 个 \(1\)、\(n_0\) 个 \(0\) 的向量。例如 \(n=5\)、\(n_1=3\) 时，\(M=\binom{5}{3}=10\)，可以枚举如下：

> permutation10 = function(n, n1){
+   M = choose(n, n1)
+   treat.index = combn(n, n1)
+   Z = matrix(0, n, M)
+   for(m in 1:M){
+     treat = treat.index[, m]
+     Z[treat, m] = 1
+   }
+   Z
+ }
> permutation10(5, 3)

于是 \(T\) 在集合 \(\{T(z_1,Y),\ldots,T(z_M,Y)\}\)（可能有重复）上均匀。这个分布完全由 CRE 的设计决定，我们叫它 \(T\) 的随机化分布（randomization distribution）。

若 \(T\) 越大越「极端」，可用尾概率度量极端程度：3

\[ p_{\mathrm{frt}} = M^{-1}\sum_{m=1}^{M} I\{T(z_m,Y)\ge T(Z,Y)\}, \tag{3.2} \]

这就是 Fisher 意义上的 \(p\) 值。原书图 3.1 画了计算 \(p_{\mathrm{frt}}\) 的过程。

式 (3.2) 里的 \(p_{\mathrm{frt}}\)，对任意检验统计量、任意结果生成过程都成立；也自然推广到别的实验——后面几章会一再遇见它。更重要的是，它在有限样本下是精确的：4 在 \(H_{0\mathrm{f}}\) 下，

\[ \operatorname{pr}(p_{\mathrm{frt}}\le u)\le u \quad\text{对所有 }0\le u\le 1. \tag{3.3} \]

实践中 \(M\) 往往大得吓人（例如 \(n=100\)、\(n_1=50\) 时 \(M>10^{29}\)），枚举全部处理向量并不现实。我们常用 Monte Carlo 近似 \(p_{\mathrm{frt}}\)（附录 A.5 复习 Monte Carlo）。具体做法是：从所有可能的处理向量中独立随机抽取，或等价地随机置换 \(Z\)，用

\[ \hat p_{\mathrm{frt}} = R^{-1}\sum_{r=1}^{R} I\{T(z_r,Y)\ge T(Z,Y)\} \tag{3.4} \]

来近似，其中 \(z_r\) 是 \(Z\) 的 \(R\) 次随机置换。Monte Carlo 误差随 \(R\) 增大而迅速变小；见习题 3.2。因为计算涉及置换，在 CRE 语境下 FRT 有时也被叫作置换检验。但在更复杂的实验里，FRT 的想法比「置换」本身更宽广。

3.3　检验统计量的几种经典选择

从上面看，FRT 对任意检验统计量都能给出有限样本精确的 \(p\) 值——这是它宽厚的地方，却不该鼓励我们随便挑一个统计量。直觉上，统计量应当对 \(H_{0\mathrm{f}}\) 可能被违背的样子「敏感」。下面回顾几种经典选择。

例 3.1（均值差）　均值差统计量为

\[ \hat\tau=\hat{\bar Y}(1)-\hat{\bar Y}(0), \]

其中

\[ \hat{\bar Y}(1)=n_1^{-1}\sum_{Z_i=1}Y_i=n_1^{-1}\sum_{i=1}^{n}Z_i Y_i, \qquad \hat{\bar Y}(0)=n_0^{-1}\sum_{Z_i=0}Y_i=n_0^{-1}\sum_{i=1}^{n}(1-Z_i)Y_i \]

分别是处理组与对照组的样本均值。在 \(H_{0\mathrm{f}}\) 下，

\[ \mathrm{E}(\hat\tau)=0, \qquad \operatorname{var}(\hat\tau)=\frac{n}{n_1 n_0}s^2, \]

其中（由引理 C.2，关于简单随机抽样）

\[ \bar Y=n^{-1}\sum_{i=1}^{n}Y_i, \qquad s^2=(n-1)^{-1}\sum_{i=1}^{n}(Y_i-\bar Y)^2. \]

又由有限总体中心极限定理（引理 C.4），\(\hat\tau\) 的随机化分布近似正态：

\[ \frac{\hat\tau}{\sqrt{\dfrac{n}{n_1 n_0}s^2}}\to\mathcal{N}(0,1) \quad\text{依分布}. \tag{3.5} \]

在 \(H_{0\mathrm{f}}\) 下 \(s^2\) 固定，因而在 FRT 里用

\[ \frac{\hat\tau}{\sqrt{\dfrac{n}{n_1 n_0}s^2}} \]

作检验统计量是等价的，也就可以据此算近似 \(p\) 值。

观测数据是 \(\{Y_i:Z_i=1\}\) 与 \(\{Y_i:Z_i=0\}\)，本质上是两样本问题。若假定结果独立同分布正态（见附录 A.4.1），经典等方差两样本 \(t\) 检验基于

\[ \frac{\hat\tau} {\sqrt{\dfrac{n}{n_1 n_0(n-2)} \Big[ \sum_{Z_i=1}\{Y_i-\hat{\bar Y}(1)\}^2 +\sum_{Z_i=0}\{Y_i-\hat{\bar Y}(0)\}^2 \Big] }} \sim t_{n-2}. \tag{3.6} \]

稍作代数（习题 3.8）可得展开：

\[ (n-1)s^2 = \sum_{Z_i=1}\{Y_i-\hat{\bar Y}(1)\}^2 +\sum_{Z_i=0}\{Y_i-\hat{\bar Y}(0)\}^2 +\frac{n_1 n_0}{n}\hat\tau^2. \tag{3.7} \]

当 \(n\) 很大时，\(\mathcal{N}(0,1)\) 与 \(t_{n-2}\)、以及 \(n-1\) 与 \(n-2\) 的差别都可忽略；在 \(H_{0\mathrm{f}}\) 下 \(\hat\tau\) 依概率趋于 0，\(\frac{n_1 n_0}{n}\hat\tau^2\) 也可渐近忽略。因此例 3.1 的近似 \(p\) 值，接近 t.test(..., var.equal=TRUE) 的 \(p\) 值。可在备择 \(\tau\ne 0\) 时，(3.7) 里多出来的那一项，可能让 FRT 比通常 \(t\) 检验更不强大；Ding (2016) 指出过这一点。

上面的讨论也说明：用 \(\hat\tau\) 做 FRT，实质上用的是合并方差，忽略了两组异方差。经典统计里，异方差正态两样本问题叫 Behrens–Fisher 问题（附录 A.4.1）。这时标准选择是下面的学生化统计量。

例 3.2（学生化统计量）　学生化统计量5 为

\[ t = \frac{\hat{\bar Y}(1)-\hat{\bar Y}(0)} {\sqrt{\dfrac{\hat S^2(1)}{n_1}+\dfrac{\hat S^2(0)}{n_0}}}, \]

其中

\[ \hat S^2(1)=(n_1-1)^{-1}\sum_{Z_i=1}\{Y_i-\hat{\bar Y}(1)\}^2, \qquad \hat S^2(0)=(n_0-1)^{-1}\sum_{Z_i=0}\{Y_i-\hat{\bar Y}(0)\}^2 \]

是两组样本方差。在 \(H_{0\mathrm{f}}\) 下，有限总体 CLT 再次给出 \(t\to\mathcal{N}(0,1)\)（依分布），近似 \(p\) 值接近 t.test(..., var.equal=FALSE)。

有一件事特别要紧：即便潜在分布并不正态，FRT 仍然为传统 \(t\) 检验（无论 var.equal 真或假）提供了理由。教科书常靠正态假设引入 \(t\) 检验，可现实里正态往往太强。只要有限总体 CLT 成立，这些程序依然可用；即便不信 CLT，也可以把 \(\hat\tau\) 与 \(t\) 放进 FRT，得到有限样本精确的 \(p\) 值。

第 8 章会从另一个角度说明：在 FRT 里用 \(t\)，对两组异方差更稳健。

下面这个统计量，则对重尾结果带来的离群点更温柔。

例 3.3（Wilcoxon 秩和）　\(\hat\tau\) 用原始结果的样本均值，抽样分布依赖方差；\(t\) 还用到样本方差，抽样分布依赖更高阶矩——因而两者都容易被离群点拽歪。

另一种流行选择，是基于合并样本秩。令 \(R_i\) 为 \(Y_i\) 在合并样本中的秩：

\[ R_i=\#\{j:Y_j\le Y_i\}. \]

Wilcoxon 秩和统计量是处理组秩之和：

\[ W=\sum_{i=1}^{n}Z_i R_i. \]

为代数简单，先假定没有结。6 合并样本秩之和固定为 \(n(n+1)/2\)，因而 \(W\) 等价于两组秩均值之差。在 \(H_{0\mathrm{f}}\) 下 \(R_i\) 固定，于是

\[ \mathrm{E}(W)=\frac{n_1(n+1)}{2}, \qquad \operatorname{var}(W)=\frac{n_1 n_0(n+1)}{12}, \]

（方差推导用到引理 C.2）。再由有限总体 CLT，

\[ \frac{\sum_{i=1}^{n}Z_i R_i-\dfrac{n_1(n+1)}{2}} {\sqrt{\dfrac{n_1 n_0(n+1)}{12}}} \to\mathcal{N}(0,1) \quad\text{依分布}. \tag{3.8} \]

据此可做渐近检验。R 中 wilcox.test 可基于 \(W-n_1(n_1+1)/2\) 计算精确与渐近 \(p\) 值。Lehmann (1975) 的渐近分析表明：用 \(W\) 做 FRT，在很宽的一类数据生成过程上，功效都还不错。

例 3.4（Kolmogorov–Smirnov 统计量）　处理改变结果的方式可以很多。很自然地，用经验分布来概括两组：

\[ \hat F_1(y)=n_1^{-1}\sum_{i=1}^{n}Z_i I(Y_i\le y), \qquad \hat F_0(y)=n_0^{-1}\sum_{i=1}^{n}(1-Z_i)I(Y_i\le y). \]

比较二者，得到著名的 Kolmogorov–Smirnov 统计量

\[ D=\max_y \bigl|\hat F_1(y)-\hat F_0(y)\bigr|. \]

\(D\) 的精确分布是道不轻松的数学题。当 \(n_1\to\infty\)、\(n_0\to\infty\) 时，

\[ \operatorname{pr}\!\left(\sqrt{\frac{n_1 n_0}{n}}D\le x\right) \to \frac{\sqrt{2\pi}}{x}\sum_{j=1}^{\infty} \exp\!\bigl(-(2j-1)^2\pi^2/(8x^2)\bigr), \]

据此可算渐近 \(p\) 值（Van der Vaart, 2000）。R 中 ks.test 可算精确与渐近 \(p\) 值。

3.4　LaLonde 实验数据的案例分析

我用 LaLonde (1986) 的实验数据来演示 FRT。数据在 Matching 包中（Sekhon, 2011）：

> library(Matching)
> data(lalonde)
> z = lalonde$treat
> y = lalonde$re78

这里 z 是是否随机进入就业培训项目的二值处理，y 是 1978 年真实收入。原书图 3.2 画出了处理组与对照组结果的直方图。

用现成函数计算观测到的检验统计量：

> tauhat = t.test(y[z == 1], y[z == 0], var.equal = TRUE)$statistic
> student = t.test(y[z == 1], y[z == 0], var.equal = FALSE)$statistic
> W = wilcox.test(y[z == 1], y[z == 0])$statistic
> D = ks.test(y[z == 1], y[z == 0])$statistic

再随机置换处理向量，用 Monte Carlo 近似四个统计量的随机化分布（存入 Tauhat、Student、Wilcox、Ks）：

> MC = 10^4
> Tauhat = Student = Wilcox = Ks = rep(0, MC)
> for(mc in 1:MC) {
+   zperm = sample(z)
+   Tauhat[mc]  = t.test(y[zperm == 1], y[zperm == 0], var.equal = TRUE)$statistic
+   Student[mc] = t.test(y[zperm == 1], y[zperm == 0], var.equal = FALSE)$statistic
+   Wilcox[mc]  = wilcox.test(y[zperm == 1], y[zperm == 0])$statistic
+   Ks[mc]      = ks.test(y[zperm == 1], y[zperm == 0])$statistic
+ }

基于 FRT 的单侧 \(p\) 值都小于 \(0.05\)：

> exact.pv = c(mean(Tauhat >= tauhat), mean(Student >= student),
+              mean(Wilcox >= W), mean(Ks >= D))
> round(exact.pv, 3)
[1] 0.002 0.002 0.006 0.040

不用 Monte Carlo，也可算渐近 \(p\) 值，同样都小于 \(0.05\)：

> asym.pv = c(
+   t.test(y[z == 1], y[z == 0], var.equal = TRUE)$p.value,
+   t.test(y[z == 1], y[z == 0], var.equal = FALSE)$p.value,
+   wilcox.test(y[z == 1], y[z == 0])$p.value,
+   ks.test(y[z == 1], y[z == 0])$p.value)
> round(asym.pv, 3)
[1] 0.005 0.008 0.011 0.046

两边的差别，来自渐近近似，也来自 t.test 与 wilcox.test 默认是双侧检验——公平比较时，前三个 \(p_{\mathrm{frt}}\) 宜乘以 2。

原书图 3.3 展示了四个统计量的随机化分布直方图，以及各自的观测值。尽管图 3.2 里结果显示远非正态，前三个统计量的正态近似仍然相当温柔地成立。一般而言，像图 3.3 这样的图，往往比单独一个 \(p\) 值更能把尖锐原假设讲清楚。近来 Bind and Rubin (2020) 在论文标题里提议：「若有可能，请报告 Fisher 精确 \(p\) 值，并展示其背后的零随机化分布。」我同意。

3.5　随机化实验与 FRT 的一点历史

3.5.1　James Lind 的实验

James Lind（1716—1794）是苏格兰医生，皇家海军卫生领域的先驱。那时坏血病是水手的主要死因之一。他做了最早有清晰记录的随机化实验之一，并在维生素 C 被发现之前，就得出柑橘能治坏血病的结论。

在 Lind (1753) 里，他描述了把 12 名坏血病患者分到六组的实验（略作简化）：

  1. 两人每天喝一夸脱苹果酒；
  2. 两人每天三次、每次二十五滴硫酸；
  3. 两人每天三次、每次两匙醋；
  4. 两人每天半品脱海水；
  5. 两人每天两个橙子、一个柠檬；
  6. 两人每天一份辛辣糊，外加大麦水。

六天后，第五组康复了，其余组没有。若把处理简化为

\[ Z_i=1(\text{单元 }i\text{ 接受了柑橘}), \]

结果简化为

\[ Y_i=1(\text{单元 }i\text{ 六天后康复}), \]

便得到一张 \(2\times 2\) 表：

\(Y=1\) \(Y=0\)
\(Z=1\) 2 0
\(Z=0\) 0 10

这是该设计下能看见的最极端表格，数据强烈暗示柑橘对坏血病有正向作用。可统计上，我们如何度量证据的强度？

顺着 FRT 的逻辑：若处理完全无效（在 \(H_{0\mathrm{f}}\) 下），这张极端表出现的概率是

\[ \binom{12}{2}^{-1}=\frac{1}{66}=0.015, \]

这就是 \(p_{\mathrm{frt}}\)。在 \(H_{0\mathrm{f}}\) 下，这像一场小小的惊讶——在 \(0.05\) 水平上，我们轻松拒绝尖锐原假设。

3.5.2　品茶女士

Fisher (1935) 描述过著名的「品茶女士」实验。7 一位女士声称，她能分辨奶茶的两种做法：先加奶，或先加茶。听起来对多数人有些古怪。作为统计学家，Fisher 设计实验，看看她究竟能不能分辨。

他做了 8 杯茶：4 杯先加奶，4 杯先加茶，以随机顺序端给她，请她挑出先加奶的 4 杯。结果可汇总为：

女士判「先奶」 女士判「先茶」 列合计
Fisher：先奶 \(X\) \(4-X\) 4
Fisher：先茶 \(4-X\) \(X\) 4
行合计 4 4 8

\(X\) 可取 \(0,1,2,3,4\)。真实实验里 \(X=4\)——最极端的数据，强烈暗示她真能分辨。证据强度又该如何度量？

在「她完全分辨不出」的原假设下，\(\binom{8}{4}=70\) 种可能顺序里，只有一种给出 \(X=4\)。于是

\[ p_{\mathrm{frt}}=\frac{1}{70}=0.014. \]

在 \(0.05\) 水平上，我们拒绝原假设。

3.5.3　Fisher 实验设计的两条原则

上面两个例子里，\(p_{\mathrm{frt}}\) 都由实验的随机化所支撑。这点明了 Fisher 实验设计的第一条原则：随机化。

同时，这两场实验，在某种意义上已是能给出统计意义的「最小」实验。若 Lind 每组只放一人，最小 \(p\) 值是

\[ \binom{6}{1}^{-1}=\frac{1}{6}=0.167; \]

若 Fisher 只做 6 杯茶（3 杯先奶、3 杯先茶），最小 \(p\) 值是

\[ \binom{6}{3}^{-1}=\frac{1}{20}=0.05. \]

于是我们几乎永远无法在 \(0.05\) 水平上拒绝原假设。这点明了第二条原则：重复（replication）——要让 FRT 有功效，实验里必须有足够多的单元。

第 5 章会谈到第三条原则：区组（blocking）。

3.6　讨论

3.6.1　其他尖锐原假设与置信区间

上面聚焦于 \(H_{0\mathrm{f}}\)。其实 FRT 的逻辑也适用于别的尖锐原假设。例如，对已知向量 \(\tau=(\tau_1,\ldots,\tau_n)\)，可检验

\[ H_0(\tau):\quad Y_i(1)-Y_i(0)=\tau_i\quad\text{对所有 }i=1,\ldots,n. \]

在 \(H_0(\tau)\) 下个体效应已知，可用观测数据补全所有缺失潜在结果；于是任意 \(T=T(Z,Y(1),Y(0))\) 的分布完全由分配机制决定，从而可算作为 \(\tau\) 函数的 \(p_{\mathrm{frt}}(\tau)\)。若能列出所有可能的 \(\tau\)，由假设检验与置信集的对偶（附录 A.2.5），可得到平均因果效应的 \((1-\alpha)\) 水平置信集：

\[ \Bigl\{\tau=n^{-1}\sum_{i=1}^{n}\tau_i : p_{\mathrm{frt}}(\tau)\ge\alpha\Bigr\}. \]

概念上直接，实践上却因 \(\tau\) 的可能数目太多而复杂。对二值结果，Rigdon and Hudgens (2015) 与 Li and Ding (2016) 提出过基于 FRT、计算上可行的置信区间构造。对一般无界结果，这策略常常算不动。

一种经典简化（Rosenbaum, 2002b, 2010）是只考虑常数个体效应的子类：

\[ H_0(c):\quad Y_i(1)-Y_i(0)=c\quad\text{对所有 }i, \]

对已知常数 \(c\)。给定 \(c\) 可算 \(p_{\mathrm{frt}}(c)\)，再由对偶得到

\[ \{c:p_{\mathrm{frt}}(c)\ge\alpha\}. \]

这只涉及一维搜索，计算可行。可「个体效应处处相等」太强——尤其对二值结果，除非所有人效应都是 \(0\)、\(-1\) 或 \(1\)，否则它根本不成立。一般地，该假设有可检验含义，可被观测数据拒绝；Ding et al. (2016) 提出过正式的统计检验。

3.6.2　其他检验统计量

FRT 是一种很宽的策略：适用于任何随机化实验、任何检验统计量。3.3 节已给了几个例子。其实统计量还可以更一般。例如，若有预处理协变量矩阵 \(X\)（第 \(i\) 行为单元 \(i\) 的 \(X_i\)），8 可允许 \(T(Z,Y,X)\) 同时依赖处理、结果与协变量。习题 3.6 有一个例子。

3.6.3　结语

对一般实验，\(Z\) 的分布未必在「\(n_1\) 个 \(1\) 与 \(n_0\) 个 \(0\) 的所有置换」上均匀；可它的分布对实验者完全已知。因此我们总能模拟它，从而得到尖锐原假设下任意检验统计量的分布。有限样本精确 \(p\) 值为

\[ p_{\mathrm{frt}}=\operatorname{pr}'\{T(Z',Y)\ge T(Z,Y)\}, \]

其中 \(\operatorname{pr}'\) 是条件于数据、对 \(Z'\) 的分布取平均。后续章节会谈别的实验；我想强调：FRT 的适用范围，比本书具体讨论的那些实验更广。

FRT 接纳任何检验统计量，却不回答实践里「该选哪一个」。若目标是相对尖锐原假设去找惊讶，就希望统计量在备择下有高功效。可一般没有「处处压倒别人」的统计量——功效取决于备择。3.3 节四个统计量，各自被不同备择轻轻推着：\(\hat\tau\) 与 \(t\) 对应非零平均处理效应；\(W\) 对应重尾结果下的常数因果效应。指定一个工作备择（working alternative），常有助于构造统计量；它不必精确，也不妨碍 FRT 的有效性。习题 3.6、3.7 会演示如何用工作备择或统计模型来构造检验统计量。

3.7　习题

3.1　\(p_{\mathrm{frt}}\) 的精确性
证明 (3.3)。

3.2　\(\hat p_{\mathrm{frt}}\) 的 Monte Carlo 误差
给定数据时，\(p_{\mathrm{frt}}\) 是固定的，而 (3.4) 中的 Monte Carlo 估计 \(\hat p_{\mathrm{frt}}\) 是随机的。证明

\[ \mathrm{E}_{\mathrm{mc}}(\hat p_{\mathrm{frt}})=p_{\mathrm{frt}}, \qquad \operatorname{var}_{\mathrm{mc}}(\hat p_{\mathrm{frt}})\le \frac{1}{4R}, \]

其中下标 “mc” 表示 Monte Carlo 带来的随机性（\(z_r\) 是从所有可能 \(Z\) 中独立抽取的 \(R\) 次）。

注： \(p_{\mathrm{frt}}\) 因 \(Z\) 随机而随机；本题条件于数据，故 \(p_{\mathrm{frt}}\) 固定。\(\hat p_{\mathrm{frt}}\) 因置换而随机。本题说明它对 \(p_{\mathrm{frt}}\) 无偏，并给出方差上界。Luo et al. (2021, Theorem 2) 有更细致的误差界。

3.3　有限样本有效的 Monte Carlo 近似
尽管 \(\hat p_{\mathrm{frt}}\) 无偏，有限 \(R\) 下它未必满足 \(\operatorname{pr}(\hat p_{\mathrm{frt}}\le u)\le u\)。下面的修正始终是有限样本有效的 \(p\) 值（Phipson and Smyth, 2010 在置换检验里指出过这个技巧）：

\[ \tilde p_{\mathrm{frt}} = \frac{1+\sum_{r=1}^{R} I\{T(z_r,Y)\ge T(Z,Y)\}}{1+R}. \]

证明对任意 \(R\)，\(\tilde p_{\mathrm{frt}}\) 始终满足 \(\operatorname{pr}(\tilde p_{\mathrm{frt}}\le u)\le u\)（对所有 \(u\in(0,1)\)）。

提示： 可用下面两条概率事实。

引理 3.1　若 \(X_1\sim\mathrm{Binomial}(R,p_1)\)、\(X_2\sim\mathrm{Binomial}(R,p_2)\) 且 \(p_1\ge p_2\)，则对一切 \(x\) 有 \(\operatorname{pr}(X_1\le x)\le\operatorname{pr}(X_2\le x)\)。

引理 3.2　若 \(p\sim\mathrm{Uniform}(0,1)\) 且 \(X\mid p\sim\mathrm{Binomial}(R,p)\)，则边际上 \(X\) 在 \(\{0,1,\ldots,R\}\) 上均匀。

3.4　Fisher 精确检验
考虑二值结果的 CRE，数据汇总为 \(2\times 2\) 表（行列合计分别为 \(n_1,n_0\) 等）。在 \(H_{0\mathrm{f}}\) 下，证明任意 \(T(n_{11},n_{10},n_{01},n_{00})\) 都是 \(n_{11}\) 与其他非随机常数的函数，且 \(n_{11}\) 的精确分布是超几何分布；请写出参数。

注： Barnard (1947)、Ding and Dasgupta (2016) 指出：二值结果 CRE 下，Fisher 精确检验（附录 A.3.1）与 FRT 等价。

3.5　品茶女士的更多细节
回顾 3.5.2 节，计算 \(\operatorname{pr}(X=k)\)，\(k=0,1,2,3,4\)。

3.6　协变量调整的 FRT
为本节 3.6.2 补充细节。3.4 节用四个统计量重分析了 LaLonde 实验数据。有额外协变量时，至少还有两种策略（假定潜在结果与协变量皆固定）：

  1. 先把结果对协变量回归，取残差作「伪结果」，再对残差定义那四个统计量并做 FRT，报告 \(p\) 值；
  2. 把「结果对处理与协变量回归」的处理系数当作检验统计量，做 FRT，报告 \(p\) 值。

为什么上面五种 \(p\) 值都是有限样本精确的？请说明理由。

3.7　广义线性模型下的 FRT
仍用习题 3.6 的数据，但把结果改成 re78 是否为正的二值指示。对处理与协变量做逻辑回归：处理系数显著吗？\(p\) 值是多少？再以该系数为检验统计量做 FRT，计算 \(p\) 值。

3.8　一个代数细节
验证 (3.7)。

3.9　推荐阅读
Bind and Rubin (2020) 是一篇较新的文章，倡导在分析复杂实验时报告 \(p\) 值，并展示相应的随机化分布。

  1. 有人会以为 CRE 是每个 \(Z_i\) 独立同分布于 \(\mathrm{Bernoulli}(\pi)\)，于是 \(n_1\sim\mathrm{Binomial}(n,\pi)\)。那其实叫 Bernoulli 随机化实验（BRE）；条件于 \((n_1,n_0)\) 之后，它就回到 CRE。第 4 章习题 4.7 会再谈 BRE。↩︎

  2. Fisher (1935) 本人并没用这种 \(H_{0\mathrm{f}}\) 写法——他那时还没用潜在结果记号。这种形式属于 Rubin (1980)。↩︎

  3. 因此 (3.2) 是单侧的。若需要双侧，一种做法是：当 \(T\) 大致以 0 为中心时，改用 \(|T|\)。↩︎

  4. 这是数理统计里对 \(p\) 值的标准定义。不等式常来自统计量的离散性；若等号成立，则原假设下 \(p\) 值服从 \(\mathrm{Uniform}(0,1)\)。离散性带来一些技术细节，证明里会出现不等式而非等式；细节留给习题 3.1。↩︎

  5. 用 \(t\) 这个记号是有意的，因为它与 \(t\) 分布有关；但不要与自由度 \(\nu\) 的 \(t_\nu\) 记号混淆。↩︎

  6. 有结时 FRT 照样可用：可给原始结果加一点小噪声，或对打结的观测用平均秩。详见 Lehmann (1975, 第 1 章第 4 节)。↩︎

  7. 这件事后来成了 Salsburg (2001) 那本现代统计史书的书名。↩︎

  8. 在因果推断里，若 \(X_i\) 不受处理影响——即若它有潜在结果 \(X_i(1)\)、\(X_i(0)\)，则必有 \(X_i(1)=X_i(0)\)——我们才叫它协变量。传统统计书常不区分处理与协变量，因为它们都出现在结果回归的右边，统称 covariates。本书区分它们，因为它们在因果故事里扮演不同角色。↩︎