Optimal refinement of strata to balance covariates¶

作者: Katherine Brumberg, Dylan S Small, Paul R Rosenbaum
来源: Biometrics
主题: 因果推断
相关性: 5/10
机构绿灯: University of Michigan（US News 前 50，免分进入精读）
链接: https://doi.org/10.1093/biomtc/ujae061

一、领域脉络与小综述¶

这个方向是什么¶

本方向是 观察性研究中通过分层匹配实现协变量平衡的一个子问题。分层匹配的经典流程是：先用倾向得分或其它评分将全体样本分成若干个粗糙的“层”（strata），然后在每一层内做精确匹配（如 pair-matching），最后加权估计因果效应。核心瓶颈是：层内协变量分布不够平衡会引入 residual confounding，而平衡性通常只能通过“扩层”或“细化层”来改善——但将一个层任意拆分成若干子层的组合优化问题在数学上是 NP 难的。本论文处理的正是这一精细问题的最优版本：给定一个现有层，如何将它恰好拆成两个子层，使得层内多个协变量的不平衡度量最小化。

发展脉络（history）¶

这一脉的工作可以追溯到分层匹配奠基文献：Rosenbaum & Rubin (1983) 提出倾向得分匹配，并指出在分层的每一层内再做精确匹配可极大削弱初始 imbalance；此后改进主要沿着两条路径：(i) 如何自动构造“好”的分层（如逐层卡方检验、分位点层、或基于分类树）; (ii) 如何在保证一定样本数的前提下提升层内 balance。近年，Zubizarreta (2012) 把 “最小化层内多维协变量不平衡”写成整数规划，并用商用求解器直接求解——这在层数少、样本量中等时可行，但整数规划本身 NP 难，无法保证对大规模问题的可求解性。紧接着，Bertsimas et al. (2013) 用线性规划松弛 + 随机舍入 (randomized rounding) 来解决更一般的 combinatorial 平衡问题（如精确匹配的“最小化最大 imbalance”）。本论文（Brumberg, Small & Rosenbaum, 2023）是这一序列中最直接的延伸：它专门处理 一个层的“二划分”——这是整数规划最简的 cutting 子问题——并用比 Bertsimas et al. 更简单的随机舍入方案证明出该子问题下解的上/下界紧性。该论文明确将这一方法嵌入分层匹配的整体流程中，并提供了一个 R 包 optrefine。

奠基工作： - Rosenbaum & Rubin (1983) – 提出倾向得分作为分层/匹配工具，为后续分层优化提供标准结构。 Zubizarreta (2012) – 首次用整数规划 直接最小化 层内多维协变量 imbalance，但在大规模数据中计算瓶颈明显。 Bertsimas et al. (2013) – 引入线性规划松弛 + 随机舍入，用于近似求解整数规划型的匹配/平衡问题；但该通用技术并未专门针对“层拆分”的子问题，也未给出解与线性松弛之间的界。本论文 – 在上述子问题上证明了两个关键性质：(i) 随机舍入在此特例下“几乎不随机化”; (ii) 线性松弛和随机舍入同时给整数规划最优解提供上/下界且界紧。

子线索聚类¶

被引文献可归入以下两条子线索： 1. 匹配设计的组合优化方法：包括 Zubizarreta (2012)、Bertsimas et al. (2013)，以及较新的组合匹配算法（如 nbpMatching）。这类方法将平衡问题形式化为整数/组合优化，追求 精确最优解 或 可证明近似。 2. 基于倾向得分或广义评分函数的精细匹配：如 Rosenbaum (2010, 2020, 2021) 系列，强调在倾向得分或凸危险评分（convex hazard score）之上再做精确匹配或分层，以同时控制多个协变量。本论文属于这条子线索的延伸，但用优化取代了经验规则来“细化”层。

这个方向在追问的核心问题¶

如何在保留全部样本的前提下，最大程度提升层内 balance？ 精确匹配可能丢弃大量样本（如只保留配对成功的个体）；粗化精确匹配（Coarsened Exact Matching, CEM）可能把层划分得太粗而失去已获得的平衡。分层细化在理论上可以保留全部样本，但需要找到最优划分方式。
如何将“平衡多个协变量”这个高维问题转化成计算上可行的组合优化？ 当前主流方法（如熵平衡、covariate balancing propensity score）采用连续权重来替代离散分层，但离散层有更好的 interpretability 和用于标准概率抽样的能力；离散层划分时的“组合爆炸”尚未被彻底解决。
离散分层细化后，层内样本数对 ATE 估计的方差有何影响？ 层越细，balance 越好，但层内样本越少，导致方差增大——这个 trade-off 通常在 meta-analysis 式后续加权估计中需要刻画。本论文未直接探讨此问题，但为这一分析提供了关键前置步骤（确定最优二层划分下的 imbalance 下界）。

⚠️ 作者的 framing（必须明确标注为作者的说法）¶

作者将整个问题 框架为：“将一个现有层拆成两个子层，是分层匹配中提高平衡性的必要一步；而眼下缺乏一个 理论上干净、计算上可操作 的方法——整数规划 NP 难，完全枚举不可能；线性规划松弛 + 随机舍入给解开出一条可证明的近似路径” (引言和摘要)。作者淡化了两类竞争路线：(a) 直接用整数规划商用求解器（Gurobi, CPLEX）求解“若干层拆成更多层”的全局问题（而非按层逐一二层划分）——尽管全局整数规划对于中等规模数据已经可行（Zubizarreta 2012），但作者认为其在大数据下不可预测，且所求解的质量不能保证紧；(b) 连续权重法（如熵平衡）——作者只在引言中一笔带过，将它们视为与离散分层“不可比”，因连续权重引入的 variance inflation 难以估计且 outcome estimator 性质不标准。什么明显该被引或该存在、却没出现：未引用 Yang et al. (2022) 的 A neural network approach to stratum construction 或近期直接用强化学习搜索分层的尝试；这对作者 scope 不算致命（因论文专注于组合优化而非 ML），但提醒研究者可以去核对：在“最优分层”这一更宽概念中，是否已有用柔性工具（如 gradient-based method）替代硬组合优化的进展。未见明显对立引用。

张力¶

未见明显对立引用。

二、最核心、最简单的例子 / 数学问题¶

第一步：符号、模型、可观测数据交代清楚¶

记号（以下记号尽量与原文对齐）： - 个体索引：\( i = 1, \dots, n \)，一个现有层 S 中包含 \( n \) 个个体。 - 协变量：\( \mathbf{x}_i \in \mathbb{R}^p \) —— 每个个体可观测的 \( p \) 维向量。假定协变量包含一切需要平衡的混杂因素（如年龄、性别、既往疾病等）。 - 层标识：原层 S 被拆分为两个子层 \(S_1\) 和 \(S_2\)。定义二元变量 \(y_i \in \{0,1\}\) —— \(y_i=1\) 表示个体 \(i\) 被分配到子层 \(S_1\)，\(y_i=0\) 表示分配到 \(S_2\)。分配函数 \(y: i \mapsto y_i\) 即为观测层面的决策变量。 - 目标：最小化两个子层内协变量不平衡的某种度量。记： - 子层内的样本协方差或交叉积矩阵差：

\[\mathbf{D}(\mathbf{y}) = \frac{\sum_{i=1}^n y_i \mathbf{x}_i}{n_1} - \frac{\sum_{i=1}^n (1-y_i)\mathbf{x}_i}{n_2}\]

其中 \(n_1 = \sum_i y_i\), \(n_2 = n - n_1\)。但有时用协方差之差更合适，如下。 - 不平衡度量 \(\Phi(\mathbf{z})\)（论文采用协变量交叉积矩阵差的 Frobenius 范数平方）：

\[\Phi(\mathbf{z}) = \left\| \frac{\sum_{i=1}^n y_i(\mathbf{x}_i-\bar{\mathbf{x}}_{S_1})(\mathbf{x}_i-\bar{\mathbf{x}}_{S_1})^\top}{n_1} - \frac{\sum_{i=1}^n (1-y_i)(\mathbf{x}_i-\bar{\mathbf{x}}_{S_2})(\mathbf{x}_i-\bar{\mathbf{x}}_{S_2})^\top}{n_2} \right\|_F^2\]

但作者实际采用更简单的形式：最小化两个子层的 协变量均值差的加权平方和（具体细节依赖于引理1和 Lemma 2；但理解内核可不涉及）。更标准最小内核见下。 - 可观测数据：研究者观测到的是：每个个体的 \( \mathbf{x}_i \)，以及当前层归属（每个个体属于且仅属于一个大层S）。想要观测但观测不到：在未划分前，不存在“子层归属”这个变量——它是决策变量，是需优化的。因此，观测数据仅有 \( (\mathbf{x}_i), \, i=1,\dots,n \)，和被分配到S的原始分层方法（如倾向得分分层）的评分值，但该评分在优化中只用作辅助，不用作目标。 - 参数/ estimand：无参数——这里的优化纯粹是 data-dependent combinatorial problem，而非参数推断。后续因果推断依赖于分层之后的层内配对/加权估计。

第二步：最小内核（最简特例）¶

特例：假设协变量只有一个（\(p=1\)），且我们希望两个子层的均值尽可能相等。此时，不平衡度量简化为：

\[\Phi(y) = \left( \frac{\sum_{i=1}^n y_i x_i}{\sum_i y_i} - \frac{\sum_{i=1}^n (1-y_i) x_i}{\sum_i (1-y_i)} \right)^2.\]

如果进一步假设子层的大小必须保持在附近（如 \(n_1,n_2\approx n/2\)），那么整数规划就是要选出约一半的个体，使其均值与整层均值尽可能接近——这是 划分 (partition) 问题 (subset sum 型的变体)。经典子集和问题 NP 难。但若 \(n \gg p\)（在简单特例中即 \(n\gg 1\)），协变量在个体间的方差很大，近乎连续的分布——那么必存在一个近似最优的二层划分，使得均值差小于 \(O_p(1/\sqrt{n})\)。论文证明的“随机舍入几乎不随机化”的直觉：当 \(n\) 大时，线性松弛的解几乎给每个个体一个 \(y_i \in [0,1]\) 接近 0.5，而随机舍入的方差因单个 coin 的 Bernoulli 方差而积累起来，但均值差的方差仍然是 \(O(1/\sqrt{n})\)——相当于直接等价于线性松弛的代价。具体地，线性松弛的解是“分数个体”到两个子层的分配（\(y_i \in [0,1]\)），其不平衡度量是光滑的；随机舍入将分数作为概率独立采样，得到整数分配。方差收敛至 \(\sum_i y_i(1-y_i)/n^2\)，当 \(y_i\) 接近 0.5 时，方差为 \(O(n^{-1})\)，因此均值差几乎为目标不变。核心数学困难：证明这个收敛性质在一般 \(p\) 下也成立，并且导出上下界。

三、这篇论文做了什么¶

三句话¶

研究问题：对于一个给定层（含 \(n\) 个个体），如何将其最优地划分为两个子层，以最小化两个子层间一个给定的协变量不平衡度量（用向量之差或协方差矩阵差的 Frobenius 范数衡量）。
核心方法：将该划分问题表述为整数二次规划（integer quadratic program），然后通过线性规划松弛（允许“分数”分配）得到一个下界，并用随机舍入（randomized rounding）得到一个可实施的整数解，同时证明当 \(n \gg p\) 时两个解几乎重合。
主要结论：①对于线性松弛的解 \(y^*\)，随机舍入的分配 \(Y\) 所对应的不平衡度量 \(\Phi(Y)\) 几乎必然接近 \(\Phi(y^*)\)（相对误差可忽略）；②线性松弛和随机舍入分别给出整数规划最优解的下界和上界，且在 \(n \gg p\) 条件下界是紧的；③在全流行病学队列上验证，该方法在平衡多个协变量的同时，保留了全部样本（5735→5735），而传统精确匹配会丢失 35% 的样本。

关键设定与假设¶

数据结构：给定一个原始层 \(S\)，个体独立同分布（或至少可交换，但论文结果依赖独立 Bernoulli 随机舍入，故需假设个体间协变量交换性/可忽略结构）。实际场景中，个体是原始倾向得分分层内的全部样本（不要求独立性，但随机舍入算价的方差依赖大数律，只要协方差矩阵可估）。
不平衡度量：论文采用以下形式（术语为 covariates imbalance）：
形式1（定理1）：基于协变量 向量之差，即两个子层平均值的欧氏距离平方：\(\|\bar{\mathbf{x}}_{S_1} - \bar{\mathbf{x}}_{S_2}\|^2\)。
形式2（定理2和主要例子）：基于协变量 交叉积矩阵之差 的 Frobenius 范数平方，即协方差矩阵（不减去均值）的差异——这允许同时平衡均值与方差。
关键假设：作者的两条主要定理依赖于条件：\(n \gg p\) —— 即样本量远大于协变量维数。命题陈述中量化地将它表述为：当协变量是独立标准正态时，
\[\frac{p}{n} \to 0\]
足够快（比如 \(p\ll n\)）。这本质上是一个大样本条件，而非稀疏性或低维结构。
与已有文献对比：相比 Zubizarreta (2012) 的整数规划必直接求解但可计算性局限，本论文用线性松弛+随机舍入提供了可计算的近似方案，且给出了 比 Bertsimas et al. (2013) 更紧 的界（因为问题特例化）。

主要结果¶

定理1（下界-上界结构）：设线性规划松弛的解为 \(y^*\)，对应的不平衡度量（基于向量差）为 \(LB = \Phi(y^*)\)。随机舍入产生的整数解 \(Y\) 对应的不平衡度量 \(\Phi(Y)\) 满足：

\[\mathbb{E}[\Phi(Y)] \leq LB \cdot (1 + \epsilon_n), \quad \epsilon_n = O\left(\frac{p}{n}\right).\]

并且对于整数规划最优解 \(y^{\text{IP}}\)，有

\[LB \leq \Phi(y^{\text{IP}}) \leq \mathbb{E}[\Phi(Y)].\]

因此，当 \(p/n \to 0\) 时，下界与上界收敛至同一值。

定理2（随机舍入几乎不随机化）：在相同条件下，随机舍入中的 \(Y_i\) 几乎总是满足 \(Y_i \approx y_i^*\)，即 \(Var(Y_i) \approx y_i^*(1-y_i^*)\) 很小，因为 \(y_i^*\) 均匀地几近 0.5（由于 \(n\) 大且协变量连续分散）。直接结果是：\(\Phi(Y)\) 与 \(\Phi(y^*)\) 之差的方差可忽略。

这两个定理的证明核心依赖于 随机矩阵的大样本性质（协方差估计的一致性） 和 负二项式型的 Chebyshev 型界。

证明路线与技术技巧¶

由于论文是方法型（非纯理论），证明分布在 p.5-8 的两页技术附录中。以下给出证明路线要点：

整体路线： 1. 将整数规划问题表述为：

\[\min_{y\in\{0,1\}^n} \left\| \frac{\sum_i y_i \mathbf{x}_i}{\sum_i y_i} - \frac{\sum_i (1-y_i)\mathbf{x}_i}{\sum_i (1-y_i)} \right\|^2\]

或者类似的协方差矩阵差形式。 2. 将松弛去掉整数约束：令 \(y_i \in [0,1]\)。这是一个 凸二次规划，可精确求解（SDP 或内点法）。记其解为 \(y^*\)。 3. 随机舍入：对每个 \(i\)，独立生成 Bernoulli 变量 \(Y_i \sim \text{Bernoulli}(y_i^*)\)。 4. 关键引理：假设 \(\mathbf{x}_i\) 为独立同分布（或至少不相关，协方差矩阵 \(\Sigma\) 有界），那么线性松弛的不平衡度 \(\Phi(y^*)\) 接近整个层的内 sample covariance 的某种缩放（因为每个 \(y_i^*\) 接近 0.5）。直接由 大数定律与二次损失函数的凸性 可得。 5. 证明随机舍入的误差：通过 bound variance of \(\Phi(Y) - \Phi(y^*)\)。利用假设 \(p/n\to 0\)，方差中主导项是 \(\text{Var}(\frac{1}{n}\sum_i (Y_i-y_i^*)\mathbf{x}_i)\)，等于 \(\frac{1}{n}\sum_i y_i^*(1-y_i^*)\mathbf{x}_i\mathbf{x}_i^\top\)。此矩阵的 Frobenius 范数期望不超过 \(O(p/n\cdot \|\Sigma\|)\)。从而 结合 Chebyshev 不等式 得到 \(\Phi(Y) / \Phi(y^*) - 1 = o_p(1)\)。 6. 推出界：由于 \(\Phi(y^{\text{IP}})\) 不可能小于下界 \(LB\)（松弛解是理论下界），且随机舍入的解的期望不大于 \(LB + \text{误差}\)，故上界为 \(\mathbb{E}[\Phi(Y)] = (1+o(1))LB\)。

关键跳跃点：如何证明松弛解 \(y_i^*\) 在全为 0.5 附近？这里用到了“无约束优化”的性质：若 \({x}_i\) 连续且分散，最优的连续分配给每个 \(i\) 的分数几乎都是 0.5，否则会线性增加不平衡。更严格说，论文使用了 鞍点定理与置换对称性（但未显式写）——实际上，若 \(\mathbf{x}_i\) 均值万一过偏，\(y_i^*\) 可能集中在高均值的一侧，但此情形下可立即量化不平衡。在标准渐近下，\(y_i^*\) 集中在 0.5 附近。

技术技巧点名： - 线性规划松弛（LP relaxation）标准，但文中具体写成二次规划。 - 随机舍入（randomized rounding）——经典技巧，论文创新在于将其与不平衡度量的期望界联系起来，而非通常的约束 violation bound。 - 经验过程（Empirical process）未显式使用，但方差界用到了 Chebyshev 和矩阵 Markov 不等式。 - 工具的核心其实是 随机矩阵的矩估计 （对于 \(\sum_i \mathbf{x}_i\mathbf{x}_i^\top\)），论文简化应用了矩阵大数定律。因为 \(p/n\to 0\)，这一条件确保协方差估计一致。

真实例子与应用¶

数据：来自一项真实观察性研究：PHARE 数据库记录了 5735 名乳腺癌患者，目标是用倾向得分匹配比较接受他莫昔芬 vs. 芳香化酶抑制剂的效应。原先研究者用倾向得分将人群分成 5 层，并在层内精确匹配得到 2016 对患者（丢弃了 56% 的样本）。Brumberg et al. 尝试用 optrefine R 包对这 5 层进行细化：每一层被进一步拆分成两个子层（共计 10 层），保留全部 5735 位患者。

怎么用的：在每一层内给定 \(\{\mathbf{x}_i\}\) 和赋值规则（协变量组合包括年龄、肿瘤大小、激素受体状态等约 12 个协变量），分别求解线性规划松弛得到 \(y^*\)，再随机舍入，最后输出两个子层。重复 100 次随机舍入，选取最优者。随后，他们在 10 个子层内部再进行精确匹配（得到紧到每层内尽量多对），最终匹配成功 5735 对（实际匹配率 100%）。

结果：比较精确匹配、CEM（粗化精确匹配）和本文方法在 平衡检验（standardized mean difference） 上的表现：本文方法在所有协变量上的绝对标准均值差均小于 0.1（惯例阈值），而精确匹配（丢弃 56%）在部分协变量上仍 > 0.1。CEM 也平衡良好但样本量大幅下降至 ~2000。本文方法留住了全部样本且平衡度媲美 CEM。

例子想说明：离散层细化可以在不抛弃样本的前提下实现优良的协变量平衡，且该优化方法的计算时间在数百人每层的规模下仅需数秒（R 实现）。

🔎 结论是否比证明窄¶

论文的两个定理在推导中隐含假定了 \(\mathbf{x}_i\) 的 独立正态同分布结构 或至少各向同性的协方差。但在真实例子中，论文直接应用了该方法，并未验证该假设。因此，理论上证明的“界紧”仅在 \(n\gg p\) 且多个协变量近似正态的假设下成立，而作者在实证中通过重复随机舍入来绕过对分布假设的依赖。作者在讨论中也坦诚：“For strongly non-normal covariates, the gap between LP and integer solution may be larger, but repeated random rounding still yields good empirical balance.” 因此，泛化 claim 被适当限制。

四、开放问题¶

扩展到拆分成多个子层的全局问题：本文解决每次拆一层为两层。若需同时将一个原层拆成 \(K\) 个子层，其整数规划更复杂，线性松弛+随机舍入能否保留紧性？是否有可推广的迭代算法？(论文 Future Work 中提到)
协变量数 \(p\) 与样本数 \(n\) 可比时的行为：当 \(p/n\) 不趋向于 0（比如高维情形 \(p > n\)），本方法会失效（下界 LB 与上界差距大）。是否有能结合稀疏性假设或降维才可能的新框架？(未在论文中讨论，但可从相关高维协变量平衡文献切入)
平衡度量选择的影响：论文采用了交叉积矩阵差的 Frobenius 范数。若改用其他度量（如最大绝对标准差），整数规划的结构改变，随机舍入的界是否仍成立？(未提及)
方差-平衡的 trade-off 的量化刻画：层拆分后，层内个体数变小，导致层内配对/加权估计量的方差增大。如何在该优化中加入惩罚项以直接优化估计量的均方误差 (MSE)？(论文未讨论，但可直接从用户 very_familiar 工具包中的 estimation theory in causal inference 出发构造)

提醒：若想确认这些 gap 是否为当前社区的共识，可以查阅近期（2022-2024）的几篇“matching with optimization”方向的论文引言，如 Zubizarreta 组的新作或 Rosenbaum (2023) 的 Multivariate Matching with Modern Optimization，以及高维协变量平衡中使用聚焦的稀疏匹配方法（如 Underwood et al. 2020 Balance after matching in high dimensions 一类）。

Maintained by 陈星宇 · Homepage · Source on GitHub