第 4 章 完全随机化实验中的 Neyman 重复抽样推断
Neyman (1923) 那篇奠基性论文，不只提出了潜在结果的记号，还在 CRE 下为平均因果效应给出了严谨的推断结果。Fisher 关心的是尖锐原假设下的 \(p\) 值；Neyman 则另辟蹊径：给出无偏点估计，再基于点估计的抽样分布，构造保守的置信区间。这一章会介绍 Neyman (1923) 的基本结果——第二部分后面的章节，几乎都要从这里经过。
4.1 有限总体量
考虑一项 CRE：\(n\) 个单元，其中 \(n_1\) 个接受处理，\(n_0\) 个接受对照。对单元 \(i=1,\ldots,n\)，有潜在结果 \(Y_i(1)\)、\(Y_i(0)\)，以及个体效应 \(\tau_i=Y_i(1)-Y_i(0)\)。潜在结果有有限总体均值
\[ \bar Y(1)=n^{-1}\sum_{i=1}^{n}Y_i(1), \qquad \bar Y(0)=n^{-1}\sum_{i=1}^{n}Y_i(0), \]
方差1
\[ S^2(1)=(n-1)^{-1}\sum_{i=1}^{n}\{Y_i(1)-\bar Y(1)\}^2, \qquad S^2(0)=(n-1)^{-1}\sum_{i=1}^{n}\{Y_i(0)-\bar Y(0)\}^2, \]
以及协方差
\[ S(1,0)=(n-1)^{-1}\sum_{i=1}^{n}\{Y_i(1)-\bar Y(1)\}\{Y_i(0)-\bar Y(0)\}. \]
个体效应的均值为
\[ \tau=n^{-1}\sum_{i=1}^{n}\tau_i=\bar Y(1)-\bar Y(0), \]
方差为
\[ S^2(\tau)=(n-1)^{-1}\sum_{i=1}^{n}(\tau_i-\tau)^2. \]
方差与协方差之间有一个温柔而有用的关系。
引理 4.1 \(2S(1,0)=S^2(1)+S^2(0)-S^2(\tau)\)。
这是基本结果，后面会用到。证明只需初等代数，留作习题 4.1。
这些固定量，都是科学表 \(\{Y_i(1),Y_i(0)\}_{i=1}^n\) 的函数。我们的兴趣是：根据 CRE 得到的数据 \((Z_i,Y_i)_{i=1}^n\)，去估计平均因果效应 \(\tau\)。
4.2 Neyman (1923) 定理
根据观测结果，可以算样本均值
\[ \hat{\bar Y}(1)=n_1^{-1}\sum_{i=1}^{n}Z_i Y_i, \qquad \hat{\bar Y}(0)=n_0^{-1}\sum_{i=1}^{n}(1-Z_i)Y_i, \]
以及样本方差
\[ \hat S^2(1)=(n_1-1)^{-1}\sum_{i=1}^{n}Z_i\{Y_i-\hat{\bar Y}(1)\}^2, \qquad \hat S^2(0)=(n_0-1)^{-1}\sum_{i=1}^{n}(1-Z_i)\{Y_i-\hat{\bar Y}(0)\}^2. \]
可对 \(S(1,0)\) 与 \(S^2(\tau)\)，却没有样本版本——因为对每个单元 \(i\)，\(Y_i(1)\) 与 \(Y_i(0)\) 从未被同时看见。Neyman (1923) 证明了下面的定理。
定理 4.1 在 CRE 下：
均值差估计量 \(\hat\tau=\hat{\bar Y}(1)-\hat{\bar Y}(0)\) 对 \(\tau\) 无偏：\(\mathrm{E}(\hat\tau)=\tau\)；
\(\hat\tau\) 的方差为
\[ \operatorname{var}(\hat\tau) = \frac{S^2(1)}{n_1}+\frac{S^2(0)}{n_0}-\frac{S^2(\tau)}{n} \tag{4.1} \]
\[ = \frac{n_0}{n_1 n}S^2(1)+\frac{n_1}{n_0 n}S^2(0)+\frac{2}{n}S(1,0); \tag{4.2} \]
- 方差估计量
\[ \hat V=\frac{\hat S^2(1)}{n_1}+\frac{\hat S^2(0)}{n_0} \]
对估计 \(\operatorname{var}(\hat\tau)\) 是保守的，因为
\[ \mathrm{E}(\hat V)-\operatorname{var}(\hat\tau)=\frac{S^2(\tau)}{n}\ge 0, \]
等号成立当且仅当对所有单元都有 \(\tau_i=\tau\)。
证明放在 4.3 节。先说清定理 4.1 里 \(\mathrm{E}(\cdot)\) 与 \(\operatorname{var}(\cdot)\) 的含义：潜在结果都是固定数字，只有处理指示 \(Z_i\) 是随机的——期望与方差，全是对 \(Z_i\) 的随机性取的，而 \(Z\) 是 \(n_1\) 个 \(1\) 与 \(n_0\) 个 \(0\) 的随机置换。原书图 4.1 画了 \(\hat\tau\) 的随机性：它在由 \(M=\binom{n}{n_1}\) 种可能分配诱导的 \(\{\hat\tau_1,\ldots,\hat\tau_M\}\) 上离散均匀。把图 4.1 与图 3.1 对照，能看见 FRT 与 Neyman (1923) 定理的关键差别：
- FRT 适用于任意检验统计量；Neyman (1923) 定理主要关于均值差。虽可为别的估计量做类似推导，对一般估计量这往往相当困难。
- 图 3.1 里观测结果向量 \(Y\) 是固定的；图 4.1 里观测结果向量 \(Y(z_m)\) 随 \(z_m\) 改变。
- 图 3.1 里的 \(T(z_m,Y)\) 都能由观测数据算出来；图 4.1 里的 \(\hat\tau_m\) 却是假想值——因为并非所有潜在结果都已知。
点估计 \(\hat\tau\) 本身很标准，可在潜在结果框架与 CRE 下，它的方差并不平凡。公式 (4.1) 不同于经典均值差方差公式2：它不仅依赖潜在结果的有限总体方差，还依赖个体效应的有限总体方差——等价地，依赖潜在结果的有限总体协方差。
可惜的是，\(S^2(\tau)\) 与 \(S(1,0)\) 无法从数据中识别——\(Y_i(1)\) 与 \(Y_i(0)\) 从未被同时观测。
公式 (4.1) 有一点让人困惑：个体效应越异质，\(\hat\tau\) 的变异反而越小。4.5.1 节会用数值例子核实 (4.1)。直觉在哪里？我用等价形式 (4.2) 来解释。比较潜在结果正相关与负相关两种情形。处理组虽是从 \(n\) 个单元中抽出的简单随机样本，某次实现里仍可能碰巧看到相对较大的处理潜在结果。假定这发生了：
- 若 \(S(1,0)>0\)，这些被处理单元的对照潜在结果也相对较大，于是对照组观测到的结果相对较小，\(\hat\tau\) 就会偏大；
- 若 \(S(1,0)<0\)，这些被处理单元的对照潜在结果相对较小，于是对照组观测结果相对较大，\(\hat\tau\) 就会偏小。
若某次实现里处理潜在结果相对较小，情形会反过来。整体上：\(\hat\tau\) 的无偏性不依赖潜在结果的相关；但在 \(S(1,0)>0\) 时，更容易看见更极端的 \(\hat\tau\)。因此，潜在结果正相关时，\(\hat\tau\) 的方差更大。
Li and Ding (2017, Theorem 5 and Proposition 3) 进一步基于有限总体 CLT，证明了 \(\hat\tau\) 的渐近正态性。
定理 4.2 令 \(n\to\infty\)、\(n_1\to\infty\)。若 \(n_1/n\) 的极限落在 \((0,1)\) 内，\(\{S^2(1),S^2(0),S(1,0)\}\) 有极限，且
\[ \max_{1\le i\le n}\{Y_i(1)-\bar Y(1)\}^2/n\to 0, \qquad \max_{1\le i\le n}\{Y_i(0)-\bar Y(0)\}^2/n\to 0, \]
则
\[ \frac{\hat\tau-\tau}{\sqrt{\operatorname{var}(\hat\tau)}}\to\mathcal{N}(0,1) \quad\text{依分布}, \]
并且 \(\hat S^2(1)\to S^2(1)\)、\(\hat S^2(0)\to S^2(0)\)（依概率）。
定理 4.2 的证明偏技术，超出本书范围。它保证：在大样本与一些正则条件下，\(\hat\tau\) 的抽样分布可用正态近似；也保证样本方差对总体方差相合，从而 Neyman (1923) 方差估计量的概率极限不小于 \(\hat\tau\) 的真实方差。这也让一类保守的大样本置信区间站得住脚：
\[ \hat\tau \pm z_{1-\alpha/2}\sqrt{\hat V}, \]
其中 \(z_{1-\alpha/2}\) 是标准正态的 \(1-\alpha/2\) 上分位数。它与标准两样本问题的置信区间渐近相同（见附录 A.4.1）。样本量足够大时，该区间覆盖 \(\tau\) 的概率至少约为 \(1-\alpha\)。由对偶，它对应检验
\[ H_{0\mathrm{n}}:\ \tau=0, \]
称为弱原假设（weak null hypothesis）。
由于「总有一个潜在结果缺失」这一根本困难，我们最多只能得到保守的方差估计。在统计里，置信区间的定义允许过度覆盖，因而也允许方差估计偏保守（见附录 A.2）。若实践中「少报处理效应」问题不大，保守性往往可以接受；可有时也会伤人——比如结果度量的是副作用：在医学实验里，少报新药副作用，会对患者健康造成严重后果。
4.3 证明
下面证明定理 4.1。
首先，\(\hat\tau\) 的无偏性来自表示
\[ \hat\tau = n_1^{-1}\sum_{i=1}^{n}Z_i Y_i(1) - n_0^{-1}\sum_{i=1}^{n}(1-Z_i)Y_i(0) \]
以及期望的线性性：
\[ \begin{aligned} \mathrm{E}(\hat\tau) &= n_1^{-1}\sum_{i=1}^{n}\mathrm{E}(Z_i)Y_i(1) - n_0^{-1}\sum_{i=1}^{n}\mathrm{E}(1-Z_i)Y_i(0)\\ &= n_1^{-1}\sum_{i=1}^{n}\frac{n_1}{n}Y_i(1) - n_0^{-1}\sum_{i=1}^{n}\frac{n_0}{n}Y_i(0) = \tau. \end{aligned} \]
其次，可把 \(\hat\tau\) 写成
\[ \hat\tau = \sum_{i=1}^{n}Z_i\Bigl\{\frac{Y_i(1)}{n_1}+\frac{Y_i(0)}{n_0}\Bigr\} - n_0^{-1}\sum_{i=1}^{n}Y_i(0). \]
由简单随机抽样的引理 C.2，
\[ \begin{aligned} \operatorname{var}(\hat\tau) &= \frac{n_1 n_0}{n(n-1)} \sum_{i=1}^{n} \Bigl\{ \frac{Y_i(1)}{n_1}+\frac{Y_i(0)}{n_0} -\frac{\bar Y(1)}{n_1}-\frac{\bar Y(0)}{n_0} \Bigr\}^2\\ &= \frac{n_0}{n_1 n}S^2(1) +\frac{n_1}{n_0 n}S^2(0) +\frac{2}{n}S(1,0). \end{aligned} \]
再由引理 4.1，也可写成
\[ \operatorname{var}(\hat\tau) = \frac{S^2(1)}{n_1}+\frac{S^2(0)}{n_0}-\frac{S^2(\tau)}{n}. \]
第三，处理组是大小为 \(n_1\) 的简单随机样本，引理 C.3 保证 \(Y_i(1)\) 的样本方差对其总体方差无偏：\(\mathrm{E}\{\hat S^2(1)\}=S^2(1)\)。同理 \(\mathrm{E}\{\hat S^2(0)\}=S^2(0)\)。因此 \(\hat V\) 对 (4.1) 的前两项无偏。
4.4 CRE 的回归分析
实践中，人们常用回归来推断平均因果效应 \(\tau\)。标准做法是：把结果对处理指示（含截距）做普通最小二乘（OLS）
\[ (\hat\alpha,\hat\beta)=\arg\min_{(a,b)}\sum_{i=1}^{n}(Y_i-a-bZ_i)^2, \]
并用处理系数 \(\hat\beta\) 作为平均因果效应的估计。可以证明
\[ \hat\beta=\hat\tau. \tag{4.3} \]
可 OLS 通常的方差估计（见附录 B 的 (B.4)，例如 R 里 lm 的输出）等于
\[ \hat V_{\mathrm{ols}} = \frac{n(n_1-1)}{(n-2)n_1 n_0}\hat S^2(1) + \frac{n(n_0-1)}{(n-2)n_1 n_0}\hat S^2(0) \approx \frac{\hat S^2(1)}{n_0}+\frac{\hat S^2(0)}{n_1}, \tag{4.4} \]
近似在 \(n_1,n_0\) 很大时成立。即便大样本下，它也不同于 \(\hat V\)。
幸运的是，Eicker–Huber–White（EHW）稳健方差估计（见附录 B 的 (B.3)）接近 \(\hat V\)：
\[ \hat V_{\mathrm{ehw}} = \frac{\hat S^2(1)}{n_1}\cdot\frac{n_1-1}{n_1} + \frac{\hat S^2(0)}{n_0}\cdot\frac{n_0-1}{n_0} \approx \frac{\hat S^2(1)}{n_1}+\frac{\hat S^2(0)}{n_0}, \tag{4.5} \]
大样本下几乎就是 \(\hat V\)。所谓 HC2 变体，则与 \(\hat V\) 完全相同。car 包的 hccm 函数可返回 EHW 稳健方差及其 HC2 变体。
习题 4.3 为 (4.3)–(4.5) 提供更多技术细节。
4.5 例子
4.5.1 模拟
先取 \(n=100\)，其中 60 个处理、40 个对照，并生成个体效应为常数的潜在结果：
n = 100
n1 = 60
n0 = 40
y0 = rexp(n)
y0 = sort(y0, decreasing = TRUE)
y1 = y0 + 1科学表固定后，反复生成 CRE，并应用定理 4.1 得到点估计、保守方差估计，以及基于正态近似的置信区间。原书图 4.2 的 (1,1) 面板是潜在结果散点图，(1,2) 面板是 \(\hat\tau-\tau\) 的直方图。
再把对照潜在结果按相反顺序排序：
y0 = sort(y0, decreasing = FALSE)重复上述模拟——对应图 4.2 第二行。
最后随机打乱对照潜在结果：
y0 = sample(y0)再重复——对应图 4.2 第三行。
重要的是：三组模拟里，潜在结果的相关不同，但边际分布相同。下表比较真实方差、平均估计方差，以及 95% 置信区间覆盖率：
| 常数效应 | 负相关 | 独立 | |
|---|---|---|---|
| 真实方差 | 0.036 | 0.007 | 0.020 |
| 估计方差 | 0.036 | 0.036 | 0.036 |
| 覆盖率 | 0.947 | 1.000 | 0.989 |
真实方差依赖潜在结果的相关：正相关对应更大的抽样方差——这核实了 (4.2)。估计方差几乎相同，因为 \(\hat V\) 的公式只依赖潜在结果的边际分布。真实与估计之间的落差，让三组覆盖率彼此不同。只有常数因果效应时，估计方差才与真实方差相同——核实定理 4.1 的第 3 点。
图 4.2 也画出了基于 CLT 的正态密度曲线，与模拟直方图非常接近，核实了定理 4.2。
4.5.2 重尾结果与正态近似的失败
定理 4.2 中 \(\hat\tau\) 的 CLT 依赖一些正则条件；潜在结果很重尾时，这些条件会被打破。可把上面的模拟稍作修改来说明：假定个体效应仍为常数，但对照潜在结果以概率 \(0.1\)、\(0.3\) 或 \(0.5\) 被 Cauchy 成分污染。例如污染概率 \(0.1\)：
eps = rbinom(n, 1, 0.1)
y0 = (1 - eps)*rexp(n) + eps*rcauchy(n)
y1 = y0 + 1原书图 4.3 与 4.4 展示了 \(\hat\tau-\tau\) 直方图及其正态近似的两次实现。重尾时，正态近似相当差；而且与图 4.2 不同，直方图对随机种子非常敏感。
4.5.3 应用
第 3.4 节用 lalonde 数据演示过 FRT。现在用同一数据演示本章理论：
> library(Matching)
> data(lalonde)
> z = lalonde$treat
> y = lalonde$re78由定理 4.1 的公式，容易算点估计与标准误：
> n1 = sum(z)
> n0 = length(z) - n1
> tauhat = mean(y[z==1]) - mean(y[z==0])
> vhat = var(y[z==1])/n1 + var(y[z==0])/n0
> sehat = sqrt(vhat)
> tauhat
[1] 1794.343
> sehat
[1] 670.9967实践中也常用 OLS 估计平均因果效应，并顺带得到标准误：
> olsfit = lm(y ~ z)
> summary(olsfit)$coef[2, 1:2]
Estimate Std. Error
1794.3431 632.8536可这个标准误，看起来比定理 4.1 给出的偏小。改用 EHW 稳健标准误即可缓解：
> library(car)
> sqrt(hccm(olsfit)[2, 2])
[1] 672.6823
> sqrt(hccm(olsfit, type = "hc0")[2, 2])
[1] 669.3155
> sqrt(hccm(olsfit, type = "hc2")[2, 2])
[1] 670.9967注意：HC2 恰好回到定理 4.1 的 \(\hat V\)——像两位老朋友在门口重新会合。
4.6 习题
4.1 引理 4.1 的证明
证明引理 4.1。
4.2 定理 4.1 的另一证明
在 CRE 下，计算 \(\operatorname{var}\{\hat{\bar Y}(1)\}\)、\(\operatorname{var}\{\hat{\bar Y}(0)\}\)、\(\operatorname{cov}\{\hat{\bar Y}(1),\hat{\bar Y}(0)\}\)，并用它们计算 \(\operatorname{var}(\hat\tau)\)。
注： 可用附录 C 的结果。
4.3 Neyman 推断与 OLS
证明 (4.3)–(4.5)。并证明 EHW 稳健方差的 HC2 变体恰好等于 \(\hat V\)。
注： 附录 B 复习了 OLS 的若干重要技术结果。
4.4 处理效应异质性
证明 \(S^2(\tau)=0\) 蕴含 \(S^2(1)=S^2(0)\)；并给出反例：\(S^2(1)=S^2(0)\) 但 \(S^2(\tau)\ne 0\)。
再证明 \(S^2(1)<S^2(0)\) 蕴含
\[ S(Y(0),\tau)=(n-1)^{-1}\sum_{i=1}^{n}\{Y_i(0)-\bar Y(0)\}(\tau_i-\tau)<0; \]
并给出反例：\(S^2(1)>S^2(0)\) 但 \(S(Y(0),\tau)<0\)。
注： 第一点说：没有处理效应异质性，则处理与对照潜在结果方差相等；逆命题不真。第二点说：若处理潜在结果方差小于对照，则个体处理效应与对照潜在结果负相关；逆命题也不真。Gerber and Green (2012, p.293) 与 Ding et al. (2019, Appendix B.3) 有相关讨论。
4.5 方差公式的更好上界
Neyman (1923) 的保守方差估计，本质上用了
\[ \operatorname{var}(\hat\tau) = \frac{S^2(1)}{n_1}+\frac{S^2(0)}{n_0}-\frac{S^2(\tau)}{n} \le \frac{S^2(1)}{n_1}+\frac{S^2(0)}{n_0}, \]
即利用 \(S^2(\tau)\ge 0\)。请证明更紧的上界
\[ \operatorname{var}(\hat\tau) \le \frac{1}{n} \Bigl\{ \sqrt{\frac{n_0}{n_1}}S(1) + \sqrt{\frac{n_1}{n_0}}S(0) \Bigr\}^2. \tag{4.6} \]
等号何时成立？
上界 (4.6) 动机另一个保守方差估计
\[ \hat V' = \frac{1}{n} \Bigl\{ \sqrt{\frac{n_0}{n_1}}\hat S(1) + \sqrt{\frac{n_1}{n_0}}\hat S(0) \Bigr\}^2. \]
4.5.1 节的模拟用了 \(\hat V\)。请再跑一遍模拟，额外比较 \(\hat V'\) 及其置信区间。
注： 证明可参考附录 A.1.4。(4.6) 还可再改进；Aronow et al. (2014) 用 Fréchet–Hoeffding 不等式给出了 \(\operatorname{var}(\hat\tau)\) 的尖锐上界。这些改进实践中很少用，主要有两个原因：比 \(\hat V\) 更复杂，而 \(\hat V\) 可用 OLS 方便实现；基于 \(\hat V\) 的置信区间在其他设定（例如结果对处理的真线性模型）下也成立，那些改进却未必。理论上有趣，实践影响有限。
4.6 Neyman (1923) 的向量版本
经典结果针对标量结果。实践中常有多个结果。现把潜在结果扩展为向量，考虑对 \(V\in\mathbb{R}^K\) 的平均因果效应
\[ \tau_V = \frac{1}{n}\sum_{i=1}^{n}\{V_i(1)-V_i(0)\}, \]
Neyman 型估计量为处理组与对照组样本均值向量之差 \(\hat\tau_V=\bar V_1-\bar V_0\)。
在 CRE 下，证明 \(\hat\tau_V\) 对 \(\tau_V\) 无偏；求出 \(\hat\tau_V\) 的协方差矩阵；并给出（可能保守的）方差估计量。
4.7 BRE 中的推断
本书对 BRE 采用如下定义。
定义 4.1（BRE） 处理指示 \(Z_i\) 独立同分布于 \(\mathrm{Bernoulli}(\pi)\)，\(n_1=\sum_{i=1}^{n}Z_i\) 人接受处理，\(n_0=\sum_{i=1}^{n}(1-Z_i)\) 人接受对照。
第一，可用 FRT 分析 BRE。在 BRE 中如何检验 \(H_{0\mathrm{f}}\)？若真实实验是 BRE，还能不能用与 CRE 相同的 FRT 程序？能则说明理由，不能则解释为什么。
第二，也可像 Neyman (1923) 对 CRE 所做的那样，为 \(\tau\) 找点估计与方差估计：
- \(\hat\tau\) 对 \(\tau\) 无偏吗？相合吗？
- 找一个对 \(\tau\) 无偏的估计量。
- 比较上述无偏估计量的方差，与 \(\hat\tau\) 的渐近方差。
注： 在 BRE 下，\(\hat\tau\) 没有有限方差，但其渐近分布的方差有限。
4.8 推荐阅读
Ding (2016) 比较了分析 CRE 的 Fisher 路径与 Neyman 路径。