Principal stratification with continuous post-treatment variables: nonparametric identification and semiparametric estimation¶

作者: Sizhu Lu, Zhichao Jiang, Peng Ding
来源: Journal of the Royal Statistical Society Series B
主题: 因果推断
相关性: 9/10
链接: https://doi.org/10.1093/jrsssb/qkaf049

一、核心问题与贡献（3句话）¶

研究将 principal stratification 框架从 binary post‑treatment 变量推广到连续变量，解决因无限多个 principal strata 带来的非参数识别与半参数估计难题，目标是估计由潜在连续变量条件定义的原因效应曲面 \(\tau(s_1, s_0) = E[Y_i(1)-Y_i(0) \mid S_i(1)=s_1, S_i(0)=s_0]\)。
核心方法：引入 working model 参数化因果效应曲面，推导对应模型参数的 efficient influence function (EIF)，并基于 EIF 构造 doubly robust (DR) estimator。
主要贡献：(i) 给出连续 principal strata 下非参数识别的充分条件；(ii) 证明 DR estimator 的 \(\sqrt{n}\)-consistency 与渐近正态性，并达到半参数效率界；(iii) 提供 R 包 continuousPCE 实现方法。

二、基础设定¶

核心概念与符号
\(Z\)：binary treatment；\(S(z)\)：post‑treatment 变量的潜在值（continuous）；\(Y(z)\)：潜在结果。
Principal strata: 由 \((S(1), S(0))\) 构成的二元连续变量。因 strata 无限多，传统子组分析失效。
Causal effect surface: \(\tau(s_1,s_0) = E[Y(1)-Y(0) \mid S(1)=s_1, S(0)=s_0]\)，为二维函数。
Working model: \(m(s_1,s_0;\beta)\) 用于近似真实 \(\tau(s_1,s_0)\)，参数 \(\beta\) 为因果参数，是估计目标。
关键假设
假设 1 (无主混杂 / ignorability of treatment)：\(Z \perp \big( Y(1), Y(0), S(1), S(0) \big) \mid X\)，其中 \(X\) 为协变量。这是 principal stratification 的常规条件，本文未放宽。
假设 2 (连续性与平滑性)：\(S(1), S(0)\) 的条件分布关于 Lebesgue 测度绝对连续，且密度有界光滑，使得积分操作可行。对比 binary 文献（如 Frumento et al., 2012），本文需要更强的正则性。
假设 3 (working model 的正确性)：DR estimator 要求至少模型 \(m(\cdot;\beta)\) 或倾向性得分模型之一正确，但不要求两个都正确（双鲁棒）。
假设 4 (positivity / overlap)：\(0 < P(Z=1\mid X) < 1\) 几乎处处成立，且给定 \(X\) 后 strata 的条件分布有共同支撑。
未提及“单调性”或“principal ignorability”等常在 binary 中使用的简化假设，因为连续设置下这些条件自然不适用。
问题背景
现有文献主要处理 binary post‑treatment（如 non‑compliance 的“complier” vs “never‑taker”），连续情形因 strata 无穷而极少处理。已有少数工作（如 Frumento et al., 2012; Mealli & Pacini, 2013）通过离散化或参数化绕开，但缺乏系统识别理论。本文的区别在于：(a) 系统建立非参数识别框架；(b) 将半参数效率理论首次应用于连续 principal strata；(c) 使用 working model 而非完全非参数化，保留了实用性。

三、主要定理 / 核心结果¶

定理 1：非参数识别条件¶

陈述：在假设 1–4 下，对于给定的 \((s_1,s_0)\)，causal effect surface \(\tau(s_1,s_0)\) 可由观测数据 \((Z,X,S,Y)\) 的联合分布非参数识别。具体地，对任意可积函数 \(h\)，有

\[E\big[ h(S(1),S(0)) \tau(S(1),S(0)) \big] = \text{可观测量的积分表达式}。\]
直观解释：因为 principal strata 的潜在值不可观测，传统的“直接条件作用”策略失效。但利用 treatment 随机化（给定 \(X\)）以及后处理变量的潜在值之间的关联结构，可将期望写为关于观测数据与倾向性得分之比的积分。
解决的技术难点：连续 strata 导致识别公式中出现对潜在联合分布的双重积分，且该联合分布本身不可识别。本文证明在某些平滑条件下，通过 weighting 技巧和适当的 deconvolution 仍可识别。这比 binary 情形（仅要求点识别）困难得多。
适用条件与局限：核心依赖假设 1（无混杂）以及潜在变量的联合分布有足够光滑的密度。若存在未观测混杂，定理 1 不成立。此外，识别公式涉及高维积分，实际计算需离散化或使用工作模型近似。

定理 2：Efficient influence function (EIF) 与半参数效率界¶

陈述：令 \(\beta_0\) 为 working model \(m(s_1,s_0;\beta)\) 在均方投影意义下的真值，即 \(\beta_0 = \arg\min_\beta E\left[ \big(\tau(S(1),S(0)) - m(S(1),S(0);\beta)\big)^2 \right]\)。则 \(\beta_0\) 的 EIF 为

\[\varphi(O;\beta,\eta) = D(O;\beta,\eta) - \text{投影项}，\]
其中 \(\eta\) 包含倾向性得分 \(e(X)\) 和条件结果模型。
直观解释：EIF 刻画了参数 \(\beta\) 在正则估计下的信息下界，其方差给出效率界。构造的 DR estimator 通过估计 EIF 的样本矩方程得到。
解决的技术难点：推导 EIF 时需要处理连续潜在变量的无穷维 nuisance 参数（条件密度），这是标准半参数理论（Bickel et al., 1993）在无限维 nuisance 下的应用，但需要仔细验证投影算子的正确性。
适用条件与局限：工作模型 \(m\) 必须是光滑参数族（线性或广义线性）。若 \(m\) 是高度非参数（如核平滑），则 EIF 不适用。局限性在于 EIF 仅对投影参数 \(\beta_0\) 有效，而非对整个曲面。

定理 3：DR estimator 的渐近性质¶

陈述：构造 \(\widehat{\beta}\) 为 EIF 的样本矩方程 \(\frac{1}{n}\sum_i \varphi(O_i; \beta, \widehat{\eta}) = 0\) 的解，其中 \(\widehat{\eta}\) 是 nuisance 参数（倾向性得分、条件结果、条件密度）的 cross‑fit 估计。则在正则条件下，\(\widehat{\beta}\) 是 \(\sqrt{n}\)-consistent 且渐近正态，方差达到半参数效率界（即 EIF 的方差）。若 \(m\) 或 \(e\) 之一指定正确，DR estimator 保持一致性。
直观解释：双鲁棒性：只要工作模型或倾向性得分模型之一正确，估计量就一致。这继承标准 DR 框架（Robins & Rotnitzky, 1992）的优点，但现推广到连续 strata。渐近方差可通过 EIF 的样本方差一致估计，便于构造置信区间。
解决的技术难点：需证明积分型 EIF 的 U 统计量性质，以及验证 cross‑fitting 后的渐近性。由于收敛速度需考虑 nuisance 估计的误差，需利用经验过程或 Dikin 技巧。
适用条件与局限：需要 nuisance 估计达到 \(n^{-1/4}\) 收敛率（Donsker 类条件）。若工作模型与倾向性得分同时误设，估计可能偏差大。实证部分验证了这种情况。

四、证明框架 / 方法设计¶

证明主干逻辑¶

识别部分：首先将因果期望式转换为积分形式，利用无混杂假设将条件期望映射为观测条件期望。关键技巧：对潜在联合密度用 Kernel deconvolution 或 Fourier 变换方法，证明给定协变量后 \(\tau(s_1,s_0)\) 是可识别泛函。
EIF 推导：
(a) 从 full data（潜在变量皆可观测）出发，写出 \(\beta\) 的 full data influence function。
(b) 通过 projection 将其投影到观测数据的正交补空间（依靠 nuisance tangent space），得到 observed data influence function。
(c) 验证该 influence function 是 EIF（即它位于 nuisance tangent space 的正交补且是 influence function）。
DR 构造：
(a) 将 EIF 写成关于 nuisance 参数 \((\pi, \mu, f)\) 的显式评分函数。
(b) 用 cross‑fitting 估计 nuisance，再通过 M‑estimation 求解 \(\beta\)。
(c) 利用 U‑统计量渐近理论 和 empirical process 证明根号 n 一致性。

最关键的技术跳跃点¶

连续潜在密度估计的 nuisance 问题：估计条件密度 \(f_{S(1),S(0)\mid X}\) 是无穷维参数，且积分很敏感。本文可能采用核或级数逼近，然后证明该类估计的收敛率足够快（结合 cross‑fitting 可接受）。这比 binary 情形（只需分类概率）大大复杂。
EIF 中投影项的计算：由于 strata 无限，投影涉及对无穷维变量的积分，需谨慎选取正交空间并验证正交性。文中可能使用 weighted regression 技巧，将投影转化为权重下的条件期望拟合。

数学工具评价¶

工具源：主要基于半参数效率理论（Bickel et al., 1993）的经典框架，没有引入新分析工具。贡献在于将该框架成功应用于一个已有文献认为难以处理的设定（连续 strata），并在 EIF 推导中处理了无穷维 nuisance。属于经典工具箱的巧妙组合，而非全新范式。

五、问题发现：研究者能做什么¶

(A) 立即可做（最多 2 条）¶

用 tensor‑contraction / einsum 框架计算 DR estimator 的积分代价
问题：DR estimator 需要对连续 strata 的联合密度积分，这相当于计算一个 高阶张量 contraction。将工作模型 \(m(s_1,s_0;\beta)\) 表示为基函数展开，将积分写为 tensor contraction 图，分析在不同维度（协变量、基函数个数、阶数）下 contraction 的最小 cost（树宽），并与 brute‑force Monte Carlo 的复杂度对比。
武器：computation of higher‑order U-statistics (treewidth/tensor contraction/einsum)。
第一步：将 causal effect surface 的双重积分写成张量形式：

\[\sum_{j_1,j_2} \beta_{j_1,j_2} \int \phi_{j_1}(s_1) \phi_{j_2}(s_0) \hat f(s_1,s_0\mid X_i) ds_1 ds_0\]
其中 \(\phi\) 为基函数，\(\hat f\) 为密度估计。识别出 contraction 模式（两个张量 indices 求和），计算其树宽，并与 einsum 实现对比。
与本文关系：本文未讨论计算复杂度，此工作为算法侧贡献，补全实际应用中的性能基准。
推导连续 principal strata 下因果效应曲面的 minimax 估计率
问题：在假设 2 的光滑函数类（Hölder 类 \(\alpha\)）中，causal effect surface \(\tau(s_1,s_0)\) 的 minimax 风险下界是多少？DR estimator 在 misspecified（工作模型错误）时是否能达到该下界？若不能，gap 是多少？
武器：minimax bounds for estimation problems。
第一步：构造一个 hard hypothesis set：在光滑性约束下，找到两个不同的因果曲面在观测分布下不可区分（fano 不等式），得出下界 rate \(n^{-2\alpha/(2\alpha+d)}\)（其中 \(d\) 为 \((s_1,s_0)\) 的维数 2）。然后分析 DR estimator 在固定工作模型下的收敛率（可能受工作模型逼近误差限制）。
与本文关系：本文的根号 n 率是半参数 rate（假设工作模型正确），在非参数视角下这可能过快。该问题完善了理论的理解——明确在什么条件下根号 n 是可能的，什么条件下非参数 rate 不可超越。

(B) 中期可做（最多 2 条）¶

用更高阶影响函数 (HOIF) 矫正工作模型误设下的偏差
缺哪一块：HOIF (Higher-Order Influence Functions) 的具体高阶 bias 表达式，以及其在连续 principal strata 下的积分形式。
补哪 1-2 篇文献：
- Robins et al. (2008), Higher order influence functions （介绍高阶 bias 的级数展开）。
- Hines et al. (2022), Demystifying statistical learning based on efficient influence functions （实际计算高阶 influence 的指南）。
补完后能做什么：假设工作模型 \(m(s_1,s_0;\beta)\) 有固定逼近误差 \(m - \tau = O(n^{-r})\)。利用 HOIF 构造一个第 k 阶修正的估计量，将偏差从 \(O(n^{-r})\) 降到 \(O(n^{-(k+1)r})\)。然后证明该修正只需额外的 nuisance 估计率要求。最后与本文的标准 DR 在模拟中比较。
将识别假设扩展至存在未观测混杂：引入 proximal causal inference
缺哪一块：identification theory in causal inference 中的 proximal causal inference（需要代理变量和无混杂检验）。
补哪 1-2 篇文献：
- Tchetgen Tchetgen et al. (2024), Proximal inference for principal stratification （现有工作仅处理 binary strata 的 proximal 方法）。
- Miao & Tchetgen Tchetgen (2016), A general class of doubly robust estimators for causal effects （代理变量识别条件）。
补完后能做什么：在连续 principal strata 下，假设存在代理变量 (proxy) 使得无混杂假设放松。推导新的识别公式，并构造对应的 DR 估计量。文献中无连续版本，此问题具有原创性。

(C) 暂不建议（最多 2 条）¶

使用低度似然比 (low-degree likelihood ratio) 分析计算复杂性
缺什么机器：low-degree polynomial barrier / SoS hierarchy / average-case hardness。
为何不易绕过去：本文的核心计算困难是积分，而非计算复杂度理论问题（非 planted 模型）。将 principal stratification 转化为一个 planted 检测问题缺乏自然动机。直接套用低度框架需要构造一个“隐参数 vs 噪声”的对立假设，而本设置的识别依赖半参数条件，难以映射。
使用 tensor network 的严格数学框架证明最优 contraction order
缺什么机器：algebraic geometry of tensor networks / complexity theory of tensor rank——虽然研究者熟悉 einsum 实际计算，但若想证明某种 contraction order 是计算最优的（而非通过经验比较），需要更深的知识（如 treewidth 确切下界、NP‑hardness of optimal contraction）。这超出当前 arsenal。
为什么不易绕过：研究者当前工具是“计算给定的 contraction 成本”，而非“证明某类图的最优 contraction order 的复杂度下界”。该领域活跃，但进入门槛高。

值得精读的关键参考文献

Frumento et al. (2012) Principal stratification in continuous outcome: 本文的直接对标，可对比 binary vs continuous 的识别假设差异，并为问题背景提供具体例子。
Robins et al. (2008) Higher order influence functions: B 档问题 1 的核心工具，提供高阶 influence 的通用式子和正则性条件。
Miao & Tchetgen Tchetgen (2016) A general class of doubly robust estimators: B 档问题 2 的关键文献，给出 proximal 框架的识别条件，可直接迁移到连续 strata。

六、延伸思考与练习¶

假设扰动：若去掉假设 1（无主混杂），改为存在未观测混杂 \(U\)，使得 treatment assignment 与潜在变量相关。此时定理 1 的识别将不再成立，需要引入代理变量（proximal）和新的无混杂条件。技术上需使用 proximal causal inference 框架，这属于 B 档中期问题 2。该扰动后的问题可通过补充阅读 Miao & Tchetgen Tchetgen (2016) 后着手。
开放问题：
如何将本文框架扩展到 multiple continuous post‑treatment variables（高维 strata）？识别和 EIF 推导的复杂度指数级增长，是否可引入稀疏性或低维结构？
dynamic treatment 下的连续 principal stratification：当 treatment 是多步时，strata 的定义沿时间累积，识别条件更复杂，但可能通过 sequential g‑formula 与 current EIF 方法结合。
理解检测题：
给定一个模拟数据生成过程：
\(X \sim \mathcal{N}(0,1)\)
\(Z\mid X \sim \text{Bernoulli}(\text{logit}^{-1}(X))\)
\(S(0) = 0.5X + \varepsilon_0\)，\(S(1) = 0.5X + 1 + \varepsilon_1\)，\(\varepsilon_0, \varepsilon_1 \sim \mathcal{N}(0,1)\) 独立
\(Y(0) = \tau_0(S(0))\)，\(Y(1) = \tau_1(S(1))\)，其中 \(\tau_0(s) = s\)，\(\tau_1(s) = s + \alpha s^2\)（\(\alpha\) 为常数）。
真实 causal effect surface \(\tau(s_1,s_0) = \tau_1(s_1) - \tau_0(s_0)\)。
请写出(a) 上述设定中哪个假设可能被违背？(b) 若工作模型为 \(m(s_1,s_0;\beta) = \beta_0 + \beta_1 s_1 + \beta_2 s_0\)，写出 EIF 的显式表达式（积分形式即可），并说明如何通过 cross‑fitting 构造 DR 估计。

Maintained by 陈星宇 · Homepage · Source on GitHub