第 5 章 随机化实验中的分层与事后分层
能区组的就区组，不能的就随机化。
——Box et al. (1978, p.103)
这是 George Box 第二句最有名的话。1 这一章，会慢慢讲清它的意思。
5.1 分层
CRE 有时会「碰巧」生出不太理想的处理分配。先看一个带离散协变量 \(X_i\in\{1,\ldots,K\}\) 的 CRE。记 \(n_{[k]}=\#\{i:X_i=k\}\)、\(\pi_{[k]}=n_{[k]}/n\) 为第 \(k\) 层的人数与比例。CRE 把 \(n_1\) 人分到处理、\(n_0\) 人分到对照，于是第 \(k\) 层内有
\[ n_{[k]1}=\#\{i:X_i=k,Z_i=1\}, \qquad n_{[k]0}=\#\{i:X_i=k,Z_i=0\} \]
人。以正概率，某个 \(n_{[k]1}\) 或 \(n_{[k]0}\) 会变成零——某一层可能只剩处理组或只剩对照组。即便都不为零，也很可能
\[ \frac{n_{[k]1}}{n_1}-\frac{n_{[k]0}}{n_0}\ne 0, \tag{5.1} \]
而且差距可以相当大。于是第 \(k\) 层在两组中的比例不同——尽管平均而言差为零（见习题 5.1）：
\[ \mathrm{E}\!\left(\frac{n_{[k]1}}{n_1}-\frac{n_{[k]0}}{n_0}\right)=0. \tag{5.2} \]
当某些层上 \(n_{[k]1}/n_1-n_{[k]0}/n_0\) 很大时，处理组与对照组会出现不太舒服的协变量失衡。这种失衡会损伤实验质量：结果差异究竟该算在处理头上，还是算在协变量失衡头上，就变得难解释了。
怎样主动避免失衡？可以做分层随机化实验（stratified randomized experiment, SRE）。
定义 5.1（SRE） 固定各 \(n_{[k]1}\)（或 \(n_{[k]0}\)）。在离散协变量 \(X\) 的 \(K\) 个层内，各自独立做一次 CRE。
农业实验里，SRE 也叫随机区组设计，层叫区组；分层随机化也叫区组随机化。2 SRE 中可行随机化的总数为
\[ \prod_{k=1}^{K}\binom{n_{[k]}}{n_{[k]1}}, \]
每种可行随机化等概率。第 \(k\) 层内接受处理的比例为
\[ e_{[k]}=\frac{n_{[k]1}}{n_{[k]}}, \]
也叫倾向得分——第三部分会把它放到舞台中央（见定义 11.1）。SRE 与 CRE 不同：其一，SRE 的可行随机化是 CRE 可行随机化的子集，故
\[ \prod_{k=1}^{K}\binom{n_{[k]}}{n_{[k]1}} < \binom{n}{n_1}; \]
其二，\(e_{[k]}\) 在 SRE 中固定，在 CRE 中随机。
对每个单元 \(i\)，有潜在结果 \(Y_i(1)\)、\(Y_i(0)\) 与个体效应 \(\tau_i=Y_i(1)-Y_i(0)\)。第 \(k\) 层的层内平均因果效应为
\[ \tau_{[k]}=n_{[k]}^{-1}\sum_{X_i=k}\tau_i. \]
总体平均因果效应
\[ \tau=n^{-1}\sum_{i=1}^{n}\tau_i=\sum_{k=1}^{K}\pi_{[k]}\tau_{[k]} \]
是层内平均效应的加权平均。
若关心 \(\tau_{[k]}\)，可把第 \(k\) 层当作一次 CRE，直接用第 3、4 章的方法。下面讨论对 \(\tau\) 的统计推断。
5.2 FRT
5.2.1 理论
与 CRE 平行，我们先谈 SRE 中的 FRT。尖锐原假设仍是
\[ H_{0\mathrm{f}}:\ Y_i(1)=Y_i(0)\quad\text{对所有 }i=1,\ldots,n. \]
FRT 的根本想法适用于任何随机化实验：可用任意检验统计量 \(T=T(Z,Y,X)\)。在 SRE 与 \(H_{0\mathrm{f}}\) 下，\(Z\) 分布已知，因而 \(T\) 分布已知。但有两件细事要小心：
- 模拟处理向量时，必须按定义 5.1，在 \(X\) 的各层内置换处理指示——这样的 FRT 有时叫条件随机化检验或条件置换检验；
- 检验统计量最好能反映 SRE 的结构。
下面给几种经典选择。
例 5.1（分层估计量） 出于估计 \(\tau\) 的动机（见 5.3 节），可在 FRT 中使用
\[ \hat\tau_S=\sum_{k=1}^{K}\pi_{[k]}\hat\tau_{[k]}, \]
其中
\[ \hat\tau_{[k]} = n_{[k]1}^{-1}\sum_{i=1}^{n}I(X_i=k,Z_i=1)Y_i - n_{[k]0}^{-1}\sum_{i=1}^{n}I(X_i=k,Z_i=0)Y_i \]
是第 \(k\) 层内的均值差。
例 5.2（学生化分层估计量） 类比两样本学生化统计量，可对分层估计量用
\[ t_S=\frac{\hat\tau_S}{\sqrt{\hat V_S}}, \qquad \hat V_S = \sum_{k=1}^{K}\pi_{[k]}^2 \Bigl( \frac{\hat S_{[k]}^2(1)}{n_{[k]1}} + \frac{\hat S_{[k]}^2(0)}{n_{[k]0}} \Bigr), \]
其中 \(\hat S_{[k]}^2(1)\)、\(\hat S_{[k]}^2(0)\) 是第 \(k\) 层内处理组与对照组的样本方差。这个形式来自 5.3 节的 Neyman 视角。
例 5.3（组合 Wilcoxon 秩和） 先在第 \(k\) 层内算 Wilcoxon 秩和 \(W_{[k]}\)（回想例 3.3），再组合为
\[ W_S=\sum_{k=1}^{K}c_{[k]}W_{[k]}. \]
Van Elteren (1960) 基于不同渐近框架与最优性准则，提出两种权重：
\[ c_{[k]}=\frac{1}{n_{[k]1}n_{[k]0}} \quad\text{或}\quad c_{[k]}=\frac{1}{n_{[k]}+1}. \]
动机偏技术；别的权重也可能合理。
例 5.4（Hodges–Lehmann 对齐秩统计量） Van Elteren (1960) 的统计量在「少量大层」时表现不错；但「许多小层」时比较不够充分，可能浪费信息。Hodges and Lehmann (1962) 建议先把结果按层中心化
\[ \tilde Y_i=Y_i-\bar Y_{[k]} \quad\text{（若 }X_i=k\text{）}, \]
再对合并后的 \((\tilde Y_1,\ldots,\tilde Y_n)\) 取秩 \((\tilde R_1,\ldots,\tilde R_n)\)，最后构造
\[ \tilde W=\sum_{i=1}^{n}Z_i\tilde R_i. \]
以上统计量在 SRE 下都可模拟精确分布；也可算均值与方差，用正态近似得 \(p\) 值。
我找了一阵，没找到关于 SRE 下 Kolmogorov–Smirnov 统计量的详细讨论。下面是我的一个提议。
例 5.5（Kolmogorov–Smirnov 统计量） 先在第 \(k\) 层内算处理与对照经验分布的最大差 \(D_{[k]}\)，再取
\[ D_S=\sum_{k=1}^{K}c_{[k]}D_{[k]} \quad\text{或}\quad D_{\max}=\max_{1\le k\le K}c_{[k]}D_{[k]}, \]
其中 \(c_{[k]}=\sqrt{n_{[k]1}n_{[k]0}/n_{[k]}}\)，由 \(n_{[k]1},n_{[k]0}\to\infty\) 时 \(D_{[k]}\) 的极限分布动机（见例 3.4）。\(D_S\) 与 \(D_{\max}\) 更适合各层都较大的情形。另一个合理选择是
\[ D=\max_y\Bigl|\sum_{k=1}^{K}\pi_{[k]}\{\hat F_{[k]1}(y)-\hat F_{[k]0}(y)\}\Bigr|, \]
它对「大层」与「许多小层」都更友善一些。
5.2.2 一个应用
Penn Bonus 实验可用来演示 SRE 中的 FRT。Koenker and Xiao (2002) 使用的数据来自按季度分层的就业培训项目，结果是就业前等待时间：
> penndata = read.table("Penn46_ascii.txt")
> z = penndata$treatment
> y = log(penndata$duration)
> block = penndata$quarter
> table(penndata$treatment, penndata$quarter)我聚焦 \(\hat\tau_S\) 与 \(W_S\)，其余统计量的 FRT 留作习题 5.6。下面函数计算 \(\hat\tau_S\) 与 \(W_S\)：
stat_SRE = function(z, y, x) {
xlevels = unique(x); K = length(xlevels)
PiK = TauK = WK = rep(0, K)
for(k in 1:K) {
xk = xlevels[k]; zk = z[x == xk]; yk = y[x == xk]
PiK[k] = length(zk)/length(z)
TauK[k] = mean(yk[zk == 1]) - mean(yk[zk == 0])
WK[k] = wilcox.test(yk[zk == 1], yk[zk == 0])$statistic
}
return(c(sum(PiK*TauK), sum(WK/PiK)))
}再按观测数据，在各层内生成一次随机处理分配：
zRandomSRE = function(z, x) {
xlevels = unique(x); K = length(xlevels)
zrandom = z
for(k in 1:K) {
xk = xlevels[k]
zrandom[x == xk] = sample(z[x == xk])
}
return(zrandom)
}有了这些，就可以模拟随机化分布并算 \(p\) 值：
> stat.obs = stat_SRE(z, y, block)
> MC = 10^3
> ... # 层内置换并重复计算
> mean(statSREMC[, 1] <= stat.obs[1])
[1] 0.002
> mean(statSREMC[, 2] <= stat.obs[2])
[1] 0.001这里用左尾概率，因为处理对结果有负向效应。原书图 5.1 展示了观测统计量与随机化分布。
5.3 Neyman 推断
5.3.1 点估计与区间估计
SRE 的推断，建立在一个温柔的事实之上：它本质上是 \(K\) 个独立的 CRE。由此可把 Neyman (1923) 的结果推广到 SRE。第 \(k\) 层内，均值差 \(\hat\tau_{[k]}\) 对 \(\tau_{[k]}\) 无偏，方差为
\[ \operatorname{var}(\hat\tau_{[k]}) = \frac{S_{[k]}^2(1)}{n_{[k]1}} + \frac{S_{[k]}^2(0)}{n_{[k]0}} - \frac{S_{[k]}^2(\tau)}{n_{[k]}}, \]
其中 \(S_{[k]}^2(1)\)、\(S_{[k]}^2(0)\)、\(S_{[k]}^2(\tau)\) 是层内潜在结果与个体效应的方差。因此分层估计量 \(\hat\tau_S=\sum_{k=1}^{K}\pi_{[k]}\hat\tau_{[k]}\) 对 \(\tau=\sum_{k=1}^{K}\pi_{[k]}\tau_{[k]}\) 无偏，且
\[ \operatorname{var}(\hat\tau_S)=\sum_{k=1}^{K}\pi_{[k]}^2\operatorname{var}(\hat\tau_{[k]}). \]
若 \(n_{[k]1}\ge 2\) 且 \(n_{[k]0}\ge 2\)，可构造保守方差估计
\[ \hat V_S = \sum_{k=1}^{K}\pi_{[k]}^2 \Bigl( \frac{\hat S_{[k]}^2(1)}{n_{[k]1}} + \frac{\hat S_{[k]}^2(0)}{n_{[k]0}} \Bigr). \]
基于 \(\hat\tau_S\) 的正态近似，可得 Wald 型 \(1-\alpha\) 置信区间
\[ \hat\tau_S\pm z_{1-\alpha/2}\sqrt{\hat V_S}. \]
从检验角度看，在 \(H_{0\mathrm{n}}:\tau=0\) 下，可将 \(t_S=\hat\tau_S/\sqrt{\hat V_S}\) 与标准正态分位数比较，得渐近 \(p\) 值。例 5.2 里的 \(t_S\) 正是它。第 8 章会说明：在 FRT 中用 \(t_S\)，在 \(H_{0\mathrm{f}}\) 下有限样本精确，在 \(H_{0\mathrm{n}}\) 下渐近有效。
\(\hat\tau_S\) 的 CLT 技术细节此处从略；Liu and Yang (2020) 给出了证明，涵盖「少量大层」与「许多小层」两种情形。下面用数值例子轻轻摸一摸。
5.3.2 数值例子
函数 Neyman_SRE 计算 SRE 下的 Neyman 点估计与方差估计（对各层算均值差与 \(\hat S^2/n\)，再按 \(\pi_{[k]}\) 加权）。
第一组模拟：\(K=5\)，每层 80 人（50 处理、30 对照）。\(10^4\) 次模拟后，
> var(TauHat)
[1] 0.002248925
> mean(VarHat)
[1] 0.002266396图 5.2 上方面板显示点估计直方图对称、钟形、围绕真参数。因个体效应为常数，平均方差估计几乎等于估计量的方差。
第二组模拟：\(K=50\)，每层 8 人（5 处理、3 对照）。同样，直方图钟形，平均方差估计与真实方差几乎重合。
最后用 Penn Bonus 实验演示：
> est = Neyman_SRE(z, y, block)
> est[1]
[1] -0.08990646
> sqrt(est[2])
[1] 0.03079775就业培训项目显著缩短了（对数）就业前等待时间。
5.3.3 比较 SRE 与 CRE
SRE 相对 CRE 有什么好处？前面从协变量平衡动机了它；下面会看到：更好的平衡，往往也带来对平均因果效应更精确的估计。为公平比较，假定对所有 \(k\) 有 \(e_{[k]}=e\)，这保证均值差等于分层估计量：
\[ \hat\tau=\hat\tau_S. \tag{5.3} \]
证明留作习题 5.2。
再比抽样方差。经典方差分析把总方差拆成层内与层间之和，于是
\[ S^2(1) = \sum_{k=1}^{K} \Bigl[ \frac{n_{[k]}-1}{n-1}S_{[k]}^2(1) + \frac{n_{[k]}}{n-1}\{\bar Y_{[k]}(1)-\bar Y(1)\}^2 \Bigr], \]
对 \(S^2(0)\)、\(S^2(\tau)\) 有类似分解。CRE 下 \(\hat\tau\) 的方差于是拆成「层内项」与「层间项」；在各 \(n_{[k]}\) 较大时，近似为层内贡献加上层间贡献。
常倾向得分假定下 \(\pi_{[k]}/n_{[k]1}=1/(ne)\) 等成立，SRE 下 \(\operatorname{var}_{\mathrm{SRE}}(\hat\tau_S)\) 恰好等于上述近似中的层内部分。于是大样本下，\(\operatorname{var}_{\mathrm{CRE}}(\hat\tau)-\operatorname{var}_{\mathrm{SRE}}(\hat\tau_S)\) 近似等于（见习题 5.4）
\[ \sum_{k=1}^{K} \frac{\pi_{[k]}}{n} \Bigl\{ \sqrt{\frac{n_0}{n_1}}\{\bar Y_{[k]}(1)-\bar Y(1)\} + \sqrt{\frac{n_1}{n_0}}\{\bar Y_{[k]}(0)-\bar Y(0)\} \Bigr\}^2 \ge 0. \tag{5.4} \]
(5.4) 为零，仅当对每个 \(k\) 括号内都为零。当协变量能预测潜在结果时，这些量通常不全为零——SRE 相对 CRE 便有了效率增益。只有协变量完全无预测力的极端情形，大样本增益才为零；那时有限样本里 SRE 甚至可能更差。这正与开篇 Box 的话彼此印证：能区组的就区组。
几点补充。第一，上面比的是抽样方差；比较估计方差，结论类似。第二，增大 \(K\) 可提高效率，但这依赖「层足够大」；实践中有权衡，不能任意增大 \(K\)。最极端是 \(n_{[k]1}=n_{[k]0}=1\)，即配对实验——第 7 章会谈。
5.4 CRE 中的事后分层
在带离散协变量 \(X\) 的 CRE 中，各层处理/对照人数随机；在 SRE 中则固定。可若条件于 \(\mathbf{n}=\{n_{[k]1},n_{[k]0}\}_{k=1}^K\) 做推断，CRE 就变成了 SRE。数学上，若 \(\mathbf{n}\) 各分量皆非零，则
\[ \operatorname{pr}_{\mathrm{CRE}}(Z=z\mid\mathbf{n}) = \Biggl(\prod_{k=1}^{K}\binom{n_{[k]}}{n_{[k]1}}\Biggr)^{-1}, \tag{5.5} \]
即 CRE 在给定 \(\mathbf{n}\) 后的条件分布，与 SRE 中 \(Z\) 的分布相同。证明留作习题 5.5。
因此，条件于 \(\mathbf{n}\)，可把带离散 \(X\) 的 CRE 当作 SRE 来分析：FRT 变成条件 FRT；Neyman 分析变成事后分层（post-stratification）：
\[ \hat\tau_{\mathrm{PS}}=\sum_{k=1}^{K}\pi_{[k]}\hat\tau_{[k]}, \]
形式与 \(\hat\tau_S\) 相同；条件于 \(\mathbf{n}\) 的方差，也与 SRE 下 \(\hat\tau_S\) 的方差相同。
Hennessy et al. (2016) 用模拟表明，条件 FRT 往往比无条件更强大；Miratrix et al. (2013) 从理论上说明，许多情形下事后分层比 \(\hat\tau\) 更有效——两者都要求 \(X\) 能预测结果。可模拟只覆盖有限数据生成过程，理论也假定各层足够大。条件 FRT 或事后分层不能走得太极端：\(K\) 越大，越容易出现某个 \(n_{[k]1}\) 或 \(n_{[k]0}\) 为零；过小或为零会大幅减少 FRT 的随机化次数，功效可能骤降。Neyman 一侧更刺眼：那时甚至无法定义 \(\hat\tau_{\mathrm{PS}}\) 与相应方差估计。
分层在设计阶段使用 \(X\)，事后分层在分析阶段使用 \(X\)——它们是使用 \(X\) 的一对对偶。渐近上，层较大时差别很小（Miratrix et al., 2013）。
5.4.1 Meinert et al. (1970) 的例子
数据来自 Meinert et al. (1970) 报告的一次 CRE（Rothman et al., 2008 也用过）。处理是 tolbutamide，对照是安慰剂。按年龄 \(<55\) 与 \(\ge 55\) 分层后，两层各自的估计、事后分层估计、以及忽略协变量的粗估计与标准误如下：
| 层 1 | 层 2 | 事后分层 | 粗估计 | |
|---|---|---|---|---|
| 估计 | \(-0.034\) | \(-0.036\) | \(-0.035\) | \(-0.045\) |
| 标准误 | \(0.031\) | \(0.060\) | \(0.032\) | \(0.033\) |
粗估计与事后分层并未给出本质不同的结论；但粗估计比两层各自的估计都更大，而事后分层落在两层估计之间——像把两层的声音，轻轻合在一起。
5.4.2 Chong et al. (2016) 的例子
Chong et al. (2016) 在秘鲁 Cajamarca 一所乡村中学，对 219 名学生做了 SRE（2009 学年）。他们先向村诊所提供铁剂补充，并培训工作人员向当面索取的青少年免费发放；再把学生随机分到三臂，观看三种视频：「足球运动员」鼓励补铁、「医生」鼓励补铁、以及完全不提铁的对照视频。实验按班级（1–5）分层。
结果之一是 2009 年第三、四季度平均成绩；重要背景协变量是基线贫血状态。本章只用「医生」对对照的子集。用 Neyman_SRE 按班级分层可得分层估计；再把班级与基线贫血交互成 \(5\times 2=10\) 层，做事后分层。比较如下：
| 估计 | 标准误 | \(t\) | \(p\) | |
|---|---|---|---|---|
| 分层 | 0.406 | 0.202 | 2.005 | 0.045 |
| 分层后再事后分层 | 0.463 | 0.190 | 2.434 | 0.015 |
事后分层给出小得多的 \(p\) 值。这个例子也说明：事后分层不只用于 CRE，也可用于已有分层、再加入额外离散协变量的 SRE。
5.5 实践问题
如何选择 \(X\) 来构造 SRE？理论上，\(X\) 应能预测潜在结果。有时实验者从预试验等背景知识里已足够清楚，选择就直截了当；有时背景不够清晰，便按后勤便利选择——例如研究区域指示、学生队列。
\(K\) 的选择是相关问题。理论上，若各层都够大，更多分层提高效率；但 \(K\) 极大时效率反而可能下降。模拟里常看见增大 \(K\) 的边际收益递减。坊间说法是：\(K=5\) 往往已够用（这个「魔法数字 5」会在第 11 章再次出现）。有人偏爱最极端的 SRE：\(K=n/2\)，即配对实验——第 7 章再谈。
若协变量是多维连续的，SRE 还能用吗？若有预试验，可先为 \(Y(0)\) 建模型，再把预测 \(\hat Y(0)\) 离散化当作 \(X\)。若没有预试验、也不想临时离散化，可用更一般的策略——再随机化（rerandomization），那是第 6 章的主题。
5.6 习题
5.1 CRE 中的协变量平衡
在 CRE 下证明 (5.2)。
5.2 常倾向得分的后果
证明 (5.3)。
5.3 常个体效应的后果
假定 \(\tau_i=\tau\) 对所有 \(i\)。考虑加权估计类 \(\hat\tau_w=\sum_{k=1}^{K}w_{[k]}\hat\tau_{[k]}\)，权重非负。求使 \(\hat\tau_w\) 对 \(\tau\) 无偏的权重条件；并在无偏估计中找出方差最小的权重。
5.4 比较 CRE 与 SRE
证明 (5.4)。
5.5 从 CRE 到 SRE
证明 (5.5)。
5.6 5.2.2 节的更多 FRT
用其他检验统计量扩展 5.2.2 节的分析。
5.7 Imbens and Rubin (2015) 中的 SRE 与 FRT
Imbens and Rubin (2015) 讨论了田纳西 STAR 实验（1985–1986）中幼儿园数据的一次 SRE：层对应学校，分析单元是教师/班级；处理 \(1\) 为小班（每师 13–17 人），\(0\) 为常规班（22–25 人）；结果是标准化平均数学成绩。原书给出了各层的 treatment 与 outcome 列表。请用 \(\hat\tau_S\)、\(W_S\) 与 \(\tilde W\) 做 FRT，并比较 \(p\) 值。
注： 本书用 \(Z\) 表示处理，Imbens and Rubin (2015) 用 \(W\)。
5.8 多中心试验
Gould (1998, Table 1) 报告了一项多中心试验的汇总数据（multicenter.csv）：29 个中心为层；患者随机分入对照、非那雄胺 1mg、非那雄胺 5mg；结果是总症状评分相对基线的变化。个体结果未公开，故无法做 FRT；但 Neyman 推断只需汇总统计量。请分别报告「1mg vs 对照」「5mg vs 对照」的点估计与方差估计。
5.9 数据再分析
再分析第 4.5.3 节的 LaLonde 数据，做 Fisher 与 Neyman 两套推断。原实验是 CRE；现假装它是 SRE：分别按种族、婚姻状态、高中文凭指示「分层」再分析，并与 CRE 下的结果比较。
5.10 推荐阅读
Miratrix et al. (2013) 为事后分层提供了扎实理论，并与分层比较。一个主要理论结果是：渐近上差别很小，尽管有限样本可以不同。