第 26 章 主分层
第二至第五部分盯的是处理对结果的因果效应，必要时调整观测到的预处理协变量。许多应用里还有处理后变量 \(M\)：它发生在处理之后，又与结果有关。一个重要的问题是：怎样恰当地用这个 \(M\)。我先举几个动机例子，再介绍 Frangakis and Rubin (2002) 基于潜在结果的表述。
26.1 动机例子
例 26.1（不依从） 随机化实验里有不依从时，可以用 \(M\) 表示接受的处理：它受处理分配 \(Z\) 影响，又影响结果 \(Y\)。这个例子里，\(M\) 与第 21 章的 \(D\) 是同一个东西。
例 26.2（因死亡截断） 针对重症患者的随机化实验里，有些患者可能在结果 \(Y\)（例如生活质量）被测量之前就去世了。这个例子里，处理后变量 \(M\) 是存活状态的二值指示。
例 26.3（失业） 就业培训项目里，单元被随机分到处理组与对照组，并报告就业状态 \(M\) 与工资 \(Y\)。这时处理后变量是就业状态的二值指示 \(M\)。
例 26.4（替代终点） 临床试验里，关心的结果（例如 30 年生存）需要又长又贵的随访。实践者转而在随访早期收集另一些容易测量的变量。这些变量叫做「替代终点」（surrogate endpoints）。一个具体例子来自 HIV 患者的临床试验：候选替代终点 \(M\) 是 CD4 细胞计数（回想：CD4 细胞是对抗感染的白细胞）。
例 26.1–26.4 有一个共同点：变量 \(M\) 发生在处理之后，又与结果有关。\(M\) 有可能处在从 \(Z\) 到 \(Y\) 的因果路径上。图 26.1(a) 画的就是这种机制。例 26.1 对应图 26.1(a)。\(M\) 也可能不在从 \(Z\) 到 \(Y\) 的因果路径上。图 26.1(b) 画的就是那种机制。例 26.2 与 26.3 对应图 26.1(b)。例 26.4 可以对应图 26.1(a) 或 (b)，取决于替代终点怎么选。
(a) \(M\) 在从 \(Z\) 到 \(Y\) 的因果路径上
U
↙ ↘
Z → M → Y
↘_____↗
(b) \(M\) 不在从 \(Z\) 到 \(Y\) 的因果路径上
U
↙ ↘
Z → M Y
↘_____↗原书图 26.1：带处理后变量 \(M\) 的因果图。\(Z\) 随机化，\(U\) 表示未测混杂。(a) 有 \(M\to Y\)；(b) 没有。两条图里 \(Z\) 都有一条弯到 \(Y\) 的直接箭头。
实践里，底层因果图可以比图 26.1 复杂得多。这一章跟着 Frangakis and Rubin (2002) 的表述走，并不假定底层因果图。
26.2 条件于处理后变量的麻烦
对付处理后变量 \(M\) 的一种天真做法，是条件于它的观测值，把它当成预处理协变量。可 \(M\) 与 \(X\) 根本不同：一般说来 \(M\) 受处理影响，\(X\) 不受。还有一条「经验法则」：评估处理对结果的平均因果效应时，分析者不该条件于任何处理后变量（Cochran, 1957; Rosenbaum, 1984）。基于潜在结果，Frangakis and Rubin (2002) 给过下面这份有洞察的解释。
为简单起见，这一章盯 CRE。
假设 26.1（带中间变量的 CRE） 我们有
\[ Z\perp\!\!\!\perp\{M(1),M(0),Y(1),Y(0),X\}. \]
条件于 \(M=m\)，我们比较
\[ \operatorname{pr}(Y\mid Z=1,M=m) \tag{26.1} \]
与
\[ \operatorname{pr}(Y\mid Z=0,M=m). \tag{26.2} \]
这个比较看起来直觉：它度量的是，给定处理后变量同一个取值，处理组与对照组结果分布的差别。若 \(M\) 是预处理协变量，这个比较给出合理的子组效应。可若 \(M\) 是处理后变量，这个比较的解释就有麻烦。在假设 26.1 下，可以把 (26.1) 与 (26.2) 里的概率改写成
\[ \begin{aligned} \operatorname{pr}(Y\mid Z=1,M=m) &= \operatorname{pr}\{Y(1)\mid Z=1,M(1)=m\} \\ &= \operatorname{pr}\{Y(1)\mid M(1)=m\} \end{aligned} \]
以及
\[ \begin{aligned} \operatorname{pr}(Y\mid Z=0,M=m) &= \operatorname{pr}\{Y(0)\mid Z=0,M(0)=m\} \\ &= \operatorname{pr}\{Y(0)\mid M(0)=m\}. \end{aligned} \]
因此，在 CRE 下，比较 (26.1) 与 (26.2)，等价于比较不同子集上 \(Y(1)\) 与 \(Y(0)\) 的分布：若 \(Z\) 影响 \(M\)，则 \(M(1)=m\) 的那些单元，与 \(M(0)=m\) 的那些单元，并不是同一群人。于是，条件于 \(M=m\) 的比较，一般并没有因果解释，除非 \(M(1)=M(0)\)。1
回到例 26.1。比较 \(\operatorname{pr}(Y\mid Z=1,M=1)\) 与 \(\operatorname{pr}(Y\mid Z=0,M=1)\)，在单调性 \(M(1)\ge M(0)\) 下，等价于：把依从者与总是接受者的处理潜在结果，去跟总是接受者的对照潜在结果比。习题 22.8 第 3 问已经指出这种分析的毛病。
回到例 26.2。若处理改善了存活状态，处理就能比对照多救一些更弱的患者。这时，\(M(1)=1\) 的那些人比 \(M(0)=1\) 的那些人更弱，于是天真的比较给出偏向对照的有偏结果。
26.3 条件于处理后变量的潜在值
Frangakis and Rubin (2002) 提议：条件于处理后变量的联合潜在值 \(U=\{M(1),M(0)\}\)，再比较
\[ \operatorname{pr}\{Y(1)\mid M(1)=m_1,M(0)=m_0\} \]
与
\[ \operatorname{pr}\{Y(0)\mid M(1)=m_1,M(0)=m_0\}, \]
其中 \((m_1,m_0)\) 取某些值。这是在同一子集——\(M(1)=m_1\) 且 \(M(0)=m_0\) 的那些单元——上，比较处理与对照下的潜在结果。Frangakis and Rubin (2002) 把这套策略叫做主分层（principal stratification），把 \(\{M(1),M(0)\}\) 看成预处理协变量。基于这个想法，可以定义
\[ \tau(m_1,m_0)=\mathrm{E}\{Y(1)-Y(0)\mid M(1)=m_1,M(0)=m_0\} \]
为子组 \(M(1)=m_1\)、\(M(0)=m_0\) 上的主层平均因果效应（principal stratification average causal effect）。\(M\) 为二值时，有四个子组
\[ \begin{cases} \tau(1,1)=\mathrm{E}\{Y(1)-Y(0)\mid M(1)=1,M(0)=1\},\\ \tau(1,0)=\mathrm{E}\{Y(1)-Y(0)\mid M(1)=1,M(0)=0\},\\ \tau(0,1)=\mathrm{E}\{Y(1)-Y(0)\mid M(1)=0,M(0)=1\},\\ \tau(0,0)=\mathrm{E}\{Y(1)-Y(0)\mid M(1)=0,M(0)=0\}. \end{cases} \tag{26.3} \]
因为 \(\{M(1),M(0)\}\) 不受处理影响，它是一个协变量，于是 \(\tau(m_1,m_0)\) 是子组因果效应。对 \(M(1)=M(0)\) 的那些子组，处理并不改中间变量，因此 \(\tau(1,1)\) 与 \(\tau(0,0)\) 度量的是分离效应（dissociative effects）。对其余 \(m_1\ne m_0\) 的子组，主层平均因果效应 \(\tau(m_1,m_0)\) 度量的是关联效应（associative effects）。这套术语来自 Frangakis and Rubin (2002)，并不假定 \(M\) 处在从 \(Z\) 到 \(Y\) 的因果路径上。若我们有图 26.1(a)，可以把分离效应读成 \(Z\) 对 \(Y\)、不经过 \(M\) 的直接效应；可关联效应并不能简单地读成 \(Z\) 对 \(Y\) 的直接或间接效应。
例 26.1（不依从） 有不依从时，(26.3) 由总是接受者、依从者、违抗者与从不接受者上的平均因果效应组成（Imbens and Angrist, 1994; Angrist et al., 1996）。
例 26.2（因死亡截断） 结果只有在患者存活时才定义得好，因此 (26.3) 里三个子组因果效应没有意义，唯一定义得好的子组效应是
\[ \tau(1,1)=\mathrm{E}\{Y(1)-Y(0)\mid M(1)=1,M(0)=1\}. \tag{26.4} \]
它叫做存活者平均因果效应（survivor average causal effect; Rubin, 2006a）。它是那些无论处理状态如何都存活的单元上，处理对结果的平均因果效应。
例 26.3（失业） 失业问题与因死亡截断同构：工资只有在单元先就业了才定义得好。因此唯一定义得好的子组效应仍是 (26.4)，也就是就业者平均因果效应。先前，Heckman (1979) 提过一个模型，现在叫做 Heckman 选择模型（Heckman Selection Model），把失业者的工资看成缺失值，用来处理工资建模里的失业。2 可 Zhang and Rubin (2003) 与 Zhang et al. (2009) 认为：在潜在结果框架下，\(\tau(1,1)\) 是更有意义的量。
例 26.4（替代终点） 直觉上，我们想通过处理对替代终点的效应，去评估处理对结果的效应。因此，好的替代终点应当满足两条：第一，若处理不影响替代终点，它也不应影响结果；第二，若处理影响替代终点，它也应影响结果。第一条被 Frangakis and Rubin (2002) 叫做「因果必要性」（causal necessity），第二条被 Gilbert and Hudgens (2008) 叫做「因果充分性」（causal sufficiency）。基于 (26.3)，对二值替代终点，因果必要性要求 \(\tau(1,1)\) 与 \(\tau(0,0)\) 为零，因果充分性要求 \(\tau(1,0)\) 与 \(\tau(0,1)\) 不为零。
\[ M_i=1(X_i^{\mathrm{T}}\beta+u_i\ge 0). \]
第二，潜在对数工资由线性模型决定：
\[ Y_i^{*}=W_i^{\mathrm{T}}\gamma+v_i, \]
并且 \(Y_i^{*}\) 只有在 \(M_i=1\) 时才作为 \(Y_i\) 被观测到。在他的两阶段模型里，协变量 \(X_i\) 与 \(W_i\) 可以不同，误差 \((u_i,v_i)\) 是相关的二元正态。
26.4 统计推断及其困难
在例 26.1 里，若有随机化、单调性与排除限制，就可以识别依从者平均因果效应 \(\tau(1,0)\)。这是第 21 章推出的关键结果。
可在另一些例子里，我们加不上排除限制。例如，例 26.2 与 26.3 里主要关心的参数是 \(\tau(1,1)\)；例 26.4 里 \(\tau(1,1)\) 与 \(\tau(0,0)\) 都关心。没有排除限制，识别主层平均因果效应就非常难。有时连单调性都加不上，于是连潜在层的比例都认不出来。
26.4.1 特例：因死亡截断且结果为二值
我用处理、存活状态、结果都是二值的简单设定，说明基于主分层做统计推断的想法，尤其是它的困难。
在假设 26.1 之外，再加上单调性。
假设 26.2（单调性） \(M(1)\ge M(0)\)。
定理 22.1 表明：在假设 26.1 与 26.2 下，三个潜在层的比例可由
\[ \begin{aligned} \pi_{(1,1)} &= \operatorname{pr}(M=1\mid Z=0), \\ \pi_{(0,0)} &= \operatorname{pr}(M=0\mid Z=1), \\ \pi_{(1,0)} &= \operatorname{pr}(M=1\mid Z=1)-\operatorname{pr}(M=1\mid Z=0) \end{aligned} \]
识别。我们的目标是识别存活者平均因果效应 \(\tau(1,1)\)。先容易认出 \(\mathrm{E}\{Y(0)\mid M(1)=1,M(0)=1\}\)，因为观测组 \((Z=0,M=1)\) 里只有存活者：
\[ \mathrm{E}\{Y(0)\mid M(1)=1,M(0)=1\}=\mathrm{E}(Y\mid Z=0,M=1). \]
关键是认出 \(\mathrm{E}\{Y(1)\mid M(1)=1,M(0)=1\}\)。观测组 \((Z=1,M=1)\) 是两层 \((1,1)\) 与 \((1,0)\) 的混合物，因此
\[\begin{align} \mathrm{E}(Y\mid Z=1,M=1) &= \frac{\pi_{(1,1)}}{\pi_{(1,1)}+\pi_{(1,0)}} \mathrm{E}\{Y(1)\mid M(1)=1,M(0)=1\} \nonumber\\ &\qquad + \frac{\pi_{(1,0)}}{\pi_{(1,1)}+\pi_{(1,0)}} \mathrm{E}\{Y(1)\mid M(1)=1,M(0)=0\}. \tag{26.5} \end{align}\]
我们有两个未知参数，(26.5) 里却只有一条方程。因此不能从 (26.5) 唯一确定 \(\mathrm{E}\{Y(1)\mid M(1)=1,M(0)=1\}\)。可 (26.5) 对关心的量仍含有信息。也就是说，按定义 18.1，\(\mathrm{E}\{Y(1)\mid M(1)=1,M(0)=1\}\) 是部分识别的。
结果 \(Y\) 为二值时，我们知道 \(\mathrm{E}\{Y(1)\mid M(1)=1,M(0)=0\}\) 夹在 \(0\) 与 \(1\) 之间，于是 \(\mathrm{E}\{Y(1)\mid M(1)=1,M(0)=1\}\) 夹在下面两条方程的解之间：
\[ \mathrm{E}(Y\mid Z=1,M=1) = \frac{\pi_{(1,1)}}{\pi_{(1,1)}+\pi_{(1,0)}} \mathrm{E}\{Y(1)\mid M(1)=1,M(0)=1\} + \frac{\pi_{(1,0)}}{\pi_{(1,1)}+\pi_{(1,0)}} \]
以及
\[ \mathrm{E}(Y\mid Z=1,M=1) = \frac{\pi_{(1,1)}}{\pi_{(1,1)}+\pi_{(1,0)}} \mathrm{E}\{Y(1)\mid M(1)=1,M(0)=1\}. \]
因此，\(\mathrm{E}\{Y(1)\mid M(1)=1,M(0)=1\}\) 的下界是
\[ \frac{\{\pi_{(1,1)}+\pi_{(1,0)}\}\mathrm{E}(Y\mid Z=1,M=1)-\pi_{(1,0)}}{\pi_{(1,1)}}, \]
上界是
\[ \frac{\{\pi_{(1,1)}+\pi_{(1,0)}\}\mathrm{E}(Y\mid Z=1,M=1)}{\pi_{(1,1)}}. \]
然后可以推出 \(\tau(1,1)\) 的界，汇总在下面的定理 26.1。
定理 26.1 在假设 26.1 与 26.2 下，若 \(Y\) 为二值，则
\[\begin{align*} &\frac{\{\pi_{(1,1)}+\pi_{(1,0)}\}\mathrm{E}(Y\mid Z=1,M=1)-\pi_{(1,0)}}{\pi_{(1,1)}} -\mathrm{E}(Y\mid Z=0,M=1) \\ &\qquad \le \tau(1,1) \\ &\qquad \le \frac{\{\pi_{(1,1)}+\pi_{(1,0)}\}\mathrm{E}(Y\mid Z=1,M=1)}{\pi_{(1,1)}} -\mathrm{E}(Y\mid Z=0,M=1). \end{align*}\]
可以用 Imbens and Manski (2004) 的置信区间去包 \(\tau(1,1)\)，分两步：第一，得到估计的下界与上界 \([\hat l,\hat u]\)，以及估计标准误 \((\mathrm{se}_l,\mathrm{se}_u)\)；第二，把置信区间构造成 \([\hat l-z_{1-\alpha}\mathrm{se}_l,\ \hat u+z_{1-\alpha}\mathrm{se}_u]\)，其中 \(z_{1-\alpha}\) 是标准正态的 \(1-\alpha\) 分位数。Imbens and Manski (2004) 置信区间的有效性依赖一些正则条件。大多数因死亡截断的问题里，这些条件成立，因为下界与上界差得相当开，并且离极端值 \(-1\) 与 \(1\) 都有距离。技术细节我略过。
总结一下：这是一道很难的题——即便样本量无穷，也不能仅凭观测数据把参数认出来。我们可以给 \(\tau(1,1)\) 推大样本界，可基于这些界的统计推断并不是标准的。若没有单调性，大样本界的形式更复杂；见 Zhang and Rubin (2003) 与 Jiang et al. (2016, 附录 A)。
表 26.1 因死亡截断的数据，* 表示死亡患者的结果
处理 \(Z=1\)
| \(Y=1\) | \(Y=0\) | 合计 | |
|---|---|---|---|
| \(M=1\) | \(54\) | \(268\) | \(322\) |
| \(M=0\) | * | * | \(109\) |
对照 \(Z=0\)
| \(Y=1\) | \(Y=0\) | 合计 | |
|---|---|---|---|
| \(M=1\) | \(59\) | \(218\) | \(277\) |
| \(M=0\) | * | * | \(152\) |
26.4.2 一个应用
我用 Yang and Small (2016) 里的数据，来自急性呼吸窘迫综合征网络（Acute Respiratory Distress Syndrome Network）研究，含 \(861\) 名肺损伤与急性呼吸窘迫综合征患者。患者被随机分到较低潮气量或传统潮气量的机械通气。结果是二值指示：第 28 天患者能否无需辅助地自主呼吸。表 26.1 汇总观测数据。
先得到潜在层的点估计：
\[ \begin{aligned} \hat\pi_{(1,1)} &= \frac{277}{277+152}=0.646, \\ \hat\pi_{(0,0)} &= \frac{109}{109+322}=0.253, \\ \hat\pi_{(1,0)} &= 1-0.646-0.253=0.101. \end{aligned} \]
存活患者上结果的样本均值是
\[ \begin{aligned} \hat{\mathrm{E}}(Y\mid Z=1,M=1) &= \frac{54}{322}=0.168, \\ \hat{\mathrm{E}}(Y\mid Z=0,M=1) &= \frac{59}{277}=0.213. \end{aligned} \]
\(\mathrm{E}\{Y(1)\mid M(1)=1,M(0)=1\}\) 的界的估计是
\[ \left[ \frac{(0.646+0.101)\times 0.168-0.101}{0.646}, \ \frac{(0.646+0.101)\times 0.168}{0.646} \right] = [0.037,0.194], \]
于是 \(\tau(1,1)\) 的界是
\[ [0.037-0.213,\ 0.194-0.213]=[-0.176,-0.019]. \]
再用自助法把抽样不确定性收进来，基于 Imbens and Manski (2004) 方法的置信区间是 \([-0.267,0.039]\)，盖住了 \(0\)。
26.4.3 推广
Zhang and Rubin (2003) 开了大样本界这条文献。Imai (2008a) 与 Lee (2009) 是两篇后续。Cheng and Small (2006) 推过多处理臂时的界。Yang and Small (2016) 用次要结果、Yang and Ding (2018a) 用更细的生存信息，去收紧存活者平均因果效应的界。
26.5 主层得分方法
没有额外假定时，我们一般只能给主层内的因果效应划界，而不能把它们认出来。要非参数识别 \(\tau(m_1,m_0)\)，就得再加假定。这些假定怎么选，并没有共识。它们不可检验，是否说得通，取决于具体应用。
有一条研究路线，立在主层可忽略性（principal ignorability）假定下的主层得分（principal score）上，与无混杂观察性研究里的因果推断平行。为简单起见，我盯强单调性的情形。
26.5.1 强单调性下的主层得分方法
假设 26.3（强单调性） \(M(0)=0\)。
与可忽略性类似，我们再假定主层可忽略性。
假设 26.4（主层可忽略性） 我们有
\[ \mathrm{E}\{Y(0)\mid M(1)=1,X\}=\mathrm{E}\{Y(0)\mid M(1)=0,X\}. \]
假设 26.4 蕴含 \(\mathrm{E}\{Y(0)\mid M(1),X\}=\mathrm{E}\{Y(0)\mid X\}\)，或者等价地，
\[ \mathrm{E}\{Y(0)M(1)\mid X\}=\mathrm{E}\{Y(0)\mid X\}\mathrm{E}\{M(1)\mid X\}, \tag{26.6} \]
也就是：给定 \(X\) 时，\(Y(0)\) 与 \(M(1)\) 均值独立（或不相关）。假设 26.1、26.3 与 26.4 保证主层内因果效应的非参数识别，汇总在下面的定理 26.2。
定理 26.2 在假设 26.1、26.3 与 26.4 下，
- \(M(1)\) 的条件概率与边际概率 \(\pi(X)=\operatorname{pr}\{M(1)=1\mid X\}\) 与 \(\pi=\operatorname{pr}\{M(1)=1\}\)，可由
\[ \pi(X)=\operatorname{pr}(M=1\mid Z=1,X) \]
以及
\[ \pi=\operatorname{pr}(M=1\mid Z=1) \]
分别识别；
- 主层平均因果效应可由
\[ \tau(1,0)=\mathrm{E}(Y\mid Z=1,M=1)-\mathrm{E}\{\pi(X)Y\mid Z=0\}/\pi \]
以及
\[ \tau(0,0)=\mathrm{E}(Y\mid Z=1,M=0)-\mathrm{E}\{(1-\pi(X))Y\mid Z=0\}/(1-\pi) \]
识别。
条件概率 \(\pi(X)=\operatorname{pr}\{M(1)=1\mid X\}\) 叫做主层得分。定理 26.2 说：\(\tau(1,0)\) 与 \(\tau(0,0)\) 可以由带适当权重的均值差识别，权重依赖主层得分。
定理 26.2 的证明 我只证 \(\tau(1,0)\) 那条。我们有
\[ \mathrm{E}(Y\mid Z=1,M=1)=\mathrm{E}\{Y(1)\mid Z=1,M(1)=1\}=\mathrm{E}\{Y(1)\mid M(1)=1\}. \]
此外，
\[\begin{align*} \mathrm{E}\{\pi(X)Y\mid Z=0\}/\pi &= \mathrm{E}\{\pi(X)Y(0)\mid Z=0\}/\pi \\ &= \mathrm{E}\{\pi(X)Y(0)\}/\pi && \text{（假设 26.1）} \\ &= \mathrm{E}\bigl[\mathrm{E}\{\pi(X)Y(0)\mid X\}\bigr]/\pi && \text{（塔性质）} \\ &= \mathrm{E}\bigl[\pi(X)\mathrm{E}\{Y(0)\mid X\}\bigr]/\pi \\ &= \mathrm{E}\bigl[\mathrm{E}\{M(1)\mid X\}\mathrm{E}\{Y(0)\mid X\}\bigr]/\pi \\ &= \mathrm{E}\bigl[\mathrm{E}\{M(1)Y(0)\mid X\}\bigr]/\pi && \text{（由 (26.6)）} \\ &= \mathrm{E}\{M(1)Y(0)\}/\pi \\ &= \mathrm{E}\{Y(0)\mid M(1)=1\}. \end{align*}\]
\(\tau(0,0)\) 那条的证明留给习题 26.1。\(\square\)
定理 26.2 提示下面这对 \(\tau(1,0)\) 与 \(\tau(0,0)\) 的简单估计量：
- 只用处理组的数据，把 \(M\) 对 \(X\) 做逻辑回归，得到 \(\hat\pi(X_i)\)；
- 用 \(\hat\pi=\sum_{i=1}^{n}Z_i M_i/\sum_{i=1}^{n}Z_i\) 去估 \(\pi\)；
- 得到矩估计：
\[ \hat\tau(1,0) = \frac{\sum_{i=1}^{n}Z_i M_i Y_i}{\sum_{i=1}^{n}Z_i M_i} - \frac{\sum_{i=1}^{n}(1-Z_i)\hat\pi(X_i)Y_i}{\hat\pi\sum_{i=1}^{n}(1-Z_i)} \]
以及
\[ \hat\tau(0,0) = \frac{\sum_{i=1}^{n}Z_i(1-M_i)Y_i}{\sum_{i=1}^{n}Z_i(1-M_i)} - \frac{\sum_{i=1}^{n}(1-Z_i)\{1-\hat\pi(X_i)\}Y_i}{(1-\hat\pi)\sum_{i=1}^{n}(1-Z_i)}; \]
- 用自助法去近似 \(\hat\tau(1,0)\) 与 \(\hat\tau(0,0)\) 的方差。
26.5.2 一个例子
下面的函数可以计算 \(\hat\tau(1,0)\) 与 \(\hat\tau(0,0)\) 的点估计。
psw = function(Z, M, Y, X) {
## probabilities of 10 and 00
pi.10 = mean(M[Z==1])
pi.00 = 1 - pi.10
## conditional probabilities of 10 and 00
ps.10 = glm(M ~ X, family = binomial,
weights = Z)$fitted.values
ps.00 = 1 - ps.10
## PCEs 10 and 00
tau.10 = mean(Y[Z==1 & M==1]) -
mean(Y[Z==0]*ps.10[Z==0])/pi.10
tau.00 = mean(Y[Z==1 & M==0]) -
mean(Y[Z==0]*ps.00[Z==0])/pi.00
c(tau.10, tau.00)
}下面的函数可以计算点估计以及自助标准误。
psw.boot = function(Z, M, Y, X, n.boot = 500){
## point estimates
point.est = psw(Z, M, Y, X)
## bootstrap standard errors
n = length(Z)
boot.est = replicate(n.boot, {
id.boot = sample(1:n, n, replace = TRUE)
psw(Z[id.boot], M[id.boot], Y[id.boot], X[id.boot, ])
})
boot.se = apply(boot.est, 1, sd)
## results
res = rbind(point.est, boot.se)
rownames(res) = c("est", "se")
colnames(res) = c("tau10", "tau00")
return(res)
}回到第 21.5 节用过的数据。先前 IV 分析假定了排除限制。现在我们放下排除限制，改假定主层可忽略性。结果如下。
> jobsdata = read.csv("jobsdata.csv")
> getX = lm(treat ~ sex + age + marital
+ + nonwhite + educ + income,
+ data = jobsdata)
> X = model.matrix(getX)[, -1]
> Z = jobsdata$treat
> M = jobsdata$comply
> Y = jobsdata$job_seek
> table(Z, M)| \(M=0\) | \(M=1\) | |
|---|---|---|
| \(Z=0\) | \(299\) | \(0\) |
| \(Z=1\) | \(228\) | \(372\) |
> psw.boot(Z, M, Y, X)| tau10 | tau00 | |
|---|---|---|
| est | \(0.169\) | \(-0.099\) |
| se | \(0.104\) | \(0.156\) |
点估计 \(\hat\tau(1,0)\) 与基于 IV 分析的那个差得不多。点估计 \(\hat\tau(0,0)\) 靠近零，标准误却很大。这个例子里，排除限制看起来是说得通的假定。
26.5.3 推广
Follmann (2000)、Hill et al. (2002)、Jo and Stuart (2009)、Jo et al. (2011) 与 Stuart and Jo (2015) 开了用主层得分去识别主层内因果效应这条文献。Ding and Lu (2017) 给这套策略提供了理论基础。他们证明了定理 26.2，以及单调性下更一般的版本；见习题 26.2。
26.6 其他方法
要在没有排除限制时估主层平均因果效应，Zhang et al. (2009) 提议用正态混合模型。可基于正态混合模型的推断可以相当脆。一条策略是：在某些限制下，用额外信息去改善推断（Ding et al., 2011; Mealli and Pacini, 2013; Mattei et al., 2013; Jiang et al., 2016）。
概念上，主分层框架对一般的 \(M\) 都适用。多值的 \(M\) 会生出许多潜在主层，连续的 \(M\) 会生出无穷多潜在主层。那些情形里，连主层的概率都不容易认出，更不用说主层平均因果效应。Jiang and Ding (2021) 综述过一些有用的策略。
26.7 习题
26.1 补完定理 26.2 的证明
证明定理 26.2 里 \(\tau(0,0)\) 那条。
26.2 单调性下的主层得分方法
本题把定理 26.2 伸出去：假设 26.3 换成假设 26.2，假设 26.4 换成下面的假设 26.5。
假设 26.5（主层可忽略性） 我们有
\[ \mathrm{E}\{Y(1)\mid M(1)=1,M(0)=0,X\} = \mathrm{E}\{Y(1)\mid M(1)=1,M(0)=1,X\} \]
以及
\[ \mathrm{E}\{Y(0)\mid M(1)=1,M(0)=0,X\} = \mathrm{E}\{Y(0)\mid M(1)=0,M(0)=0,X\}. \]
定理 26.3 在假设 26.1、26.2 与 26.5 下，
- 条件主层得分与边际主层得分可由
\[ \begin{aligned} \pi_{(0,0)}(X) &= \operatorname{pr}(M=0\mid Z=1,X), \\ \pi_{(1,1)}(X) &= \operatorname{pr}(M=1\mid Z=0,X), \\ \pi_{(1,0)}(X) &= \operatorname{pr}(M=1\mid Z=1,X)-\operatorname{pr}(M=1\mid Z=0,X) \end{aligned} \]
以及
\[ \begin{aligned} \pi_{(0,0)} &= \operatorname{pr}(M=0\mid Z=1), \\ \pi_{(1,1)} &= \operatorname{pr}(M=1\mid Z=0), \\ \pi_{(1,0)} &= \operatorname{pr}(M=1\mid Z=1)-\operatorname{pr}(M=1\mid Z=0) \end{aligned} \]
分别识别；
- 主层平均因果效应可由
\[ \begin{aligned} \tau(1,0) &= \mathrm{E}\bigl\{w_{1,(1,0)}(X)Y\mid Z=1,M=1\bigr\} - \mathrm{E}\bigl\{w_{0,(1,0)}(X)Y\mid Z=0,M=0\bigr\}, \\ \tau(0,0) &= \mathrm{E}(Y\mid Z=1,M=0) - \mathrm{E}\bigl\{w_{0,(0,0)}(X)Y\mid Z=0,M=0\bigr\}, \\ \tau(1,1) &= \mathrm{E}\bigl\{w_{1,(1,1)}(X)Y\mid Z=1,M=1\bigr\} - \mathrm{E}(Y\mid Z=0,M=1) \end{aligned} \]
识别，其中
\[ \begin{aligned} w_{1,(1,0)}(X) &= \frac{\pi_{(1,0)}(X)}{\pi_{(1,0)}(X)+\pi_{(1,1)}(X)} \bigg/ \frac{\pi_{(1,0)}}{\pi_{(1,0)}+\pi_{(1,1)}}, \\ w_{0,(1,0)}(X) &= \frac{\pi_{(1,0)}(X)}{\pi_{(1,0)}(X)+\pi_{(0,0)}(X)} \bigg/ \frac{\pi_{(1,0)}}{\pi_{(1,0)}+\pi_{(0,0)}}, \\ w_{0,(0,0)}(X) &= \frac{\pi_{(0,0)}(X)}{\pi_{(1,0)}(X)+\pi_{(0,0)}(X)} \bigg/ \frac{\pi_{(0,0)}}{\pi_{(1,0)}+\pi_{(0,0)}}, \\ w_{1,(1,1)}(X) &= \frac{\pi_{(1,1)}(X)}{\pi_{(1,0)}(X)+\pi_{(1,1)}(X)} \bigg/ \frac{\pi_{(1,1)}}{\pi_{(1,0)}+\pi_{(1,1)}}. \end{aligned} \]
注： 基于定理 26.3，可以构造加权估计量。定理 26.3 是 Ding and Lu (2017) 的命题 2，原文还给了更细的估计细节。
26.3 观察性研究里的主层得分方法
本题把定理 26.2 伸出去：假设 26.1 换成下面的可忽略性。
假设 26.6 \(Z\perp\!\!\!\perp\{M(1),M(0),Y(1),Y(0)\}\mid X\)。
回想倾向得分的定义 \(e(X)=\operatorname{pr}(Z=1\mid X)\)。我们有下面的识别结果。
定理 26.4 在假设 26.6、26.3 与 26.4 下，
- \(M(1)\) 的条件概率与边际概率 \(\pi(X)=\operatorname{pr}\{M(1)=1\mid X\}\) 与 \(\pi=\operatorname{pr}\{M(1)=1\}\)，可由
\[ \pi(X)=\operatorname{pr}(M=1\mid Z=1,X) \]
以及
\[ \pi=\mathrm{E}\{\operatorname{pr}(M=1\mid Z=1,X)\} \]
分别识别；
- 主层平均因果效应可由
\[ \tau(1,0) = \mathrm{E}\left\{\frac{M}{\pi}\frac{Z}{e(X)}Y\right\} - \mathrm{E}\left\{\frac{\pi(X)}{\pi}\frac{1-Z}{1-e(X)}Y\right\} \]
以及
\[ \tau(0,0) = \mathrm{E}\left\{\frac{1-M}{1-\pi}\frac{Z}{e(X)}Y\right\} - \mathrm{E}\left\{\frac{1-\pi(X)}{1-\pi}\frac{1-Z}{1-e(X)}Y\right\} \]
识别。
证明定理 26.4。
26.4 一般的主层得分方法
在假设 26.6、26.2 与 26.5 下，把定理 26.3 与 26.4 伸出去。
注： 见 Jiang et al. (2022)。
26.5 推荐阅读
Frangakis and Rubin (2002) 提出主分层框架。Zhang and Rubin (2003) 推过存活者平均因果效应的大样本界。Jiang and Ding (2021) 综述过识别主层内因果效应的各种策略。Jiang et al. (2022) 给观察性研究里这套策略做过统一讨论，并提议主层内因果效应的多重稳健估计量。