第 16 章 观察性研究因果效应中无混杂假定的困难
第三部分一直站在两条假定上：无混杂与重叠。两条都强，实践里都容易破。这一章先谈无混杂为什么难；第 17–19 章谈未测混杂时的各种敏感性分析；第 20 章再谈重叠。
16.1 因果图的一点基础
Pearl (1995) 把因果图（causal diagram）引进经验研究，当作因果推断的一件利器。Pearl (2000) 是因果图的教材。这里我只用它当一件直觉工具，把变量之间的因果关系画出来。
例如，若有因果图
X
↙ ↘
Z → Y并关心 \(Z\) 对 \(Y\) 的因果效应，就可以把它读成下面的数据生成过程：
\[ \begin{cases} X\sim F_X(x),\\ Z=f_Z(X,\varepsilon_Z),\\ Y(z)=f_Y(X,z,\varepsilon_{Y(z)}), \end{cases} \]
其中对 \(z=0,1\) 都有 \(\varepsilon_Z\perp\!\!\!\perp\varepsilon_{Y(z)}\)。上面说的是：协变量 \(X\) 从分布 \(F_X(x)\) 里来；处理分配是 \(X\) 的函数，再加随机误差 \(\varepsilon_Z\)；潜在结果 \(Y(z)\) 是 \(X\)、\(z\) 与随机误差 \(\varepsilon_{Y(z)}\) 的函数。从这些方程很容易读出 \(Z\perp\!\!\!\perp Y(z)\mid X\)，也就是无混杂假定成立。
若因果图变成
X U
↓ ↘ ↙ ↓
Z → Y就可以读成：
\[ \begin{cases} X\sim F_X(x),\\ U\sim F_U(u),\\ Z=f_Z(X,U,\varepsilon_Z),\\ Y(z)=f_Y(X,U,z,\varepsilon_{Y(z)}), \end{cases} \]
其中对 \(z=0,1\) 都有 \(\varepsilon_Z\perp\!\!\!\perp\varepsilon_{Y(z)}\)。从方程很容易读出：\(Z\perp\!\!\!\perp Y(z)\mid(X,U)\)，但 \(Z\not\perp\!\!\!\perp Y(z)\mid X\)。也就是说，无混杂在给定 \((X,U)\) 时成立，只给定 \(X\) 时不成立。这时 \(U\) 是未测混杂变量（unmeasured confounder）。
16.2 评估无混杂假定
无混杂假定
\[ Z\perp\!\!\!\perp Y(1)\mid X, \qquad Z\perp\!\!\!\perp Y(0)\mid X \]
蕴含
\[ \operatorname{pr}\{Y(1)\mid Z=1,X\} = \operatorname{pr}\{Y(1)\mid Z=0,X\}, \]
\[ \operatorname{pr}\{Y(0)\mid Z=1,X\} = \operatorname{pr}\{Y(0)\mid Z=0,X\}. \]
所以它基本上要求：反事实分布 \(\operatorname{pr}\{Y(1)\mid Z=0,X\}\) 等于观测分布 \(\operatorname{pr}\{Y(1)\mid Z=1,X\}\)，反事实分布 \(\operatorname{pr}\{Y(0)\mid Z=1,X\}\) 等于观测分布 \(\operatorname{pr}\{Y(0)\mid Z=0,X\}\)。反事实分布并不能从数据里直接认出，因而若不加额外假定，无混杂从根本上不可检验。下面谈两种评估无混杂的策略。这里「评估」比「检验」弱：前者是补充分析，给最初的分析撑腰或拆台；后者才是正式的统计检验。
16.2.1 用阴性结局
假定 \(Y^{\mathrm{n}}\) 是一个与 \(Y\) 相似的结局，理想情形下与 \(Y\) 分享同一套混杂结构。若我们相信 \(Z\perp\!\!\!\perp Y(z)\mid X\)，往往也会相信 \(Z\perp\!\!\!\perp Y^{\mathrm{n}}(z)\mid X\)。而且，我们事先知道 \(Z\) 对 \(Y^{\mathrm{n}}\) 的效应：
\[ \tau(Z\to Y^{\mathrm{n}})=\mathrm{E}\{Y^{\mathrm{n}}(1)-Y^{\mathrm{n}}(0)\}. \]
一个要紧的例子是 \(\tau(Z\to Y^{\mathrm{n}})=0\)。满足这些要求的因果图如下：
X → Yⁿ
↙ ↘
Z → Y例 16.1 Cornfield et al. (1959) 用观察性研究，关心吸烟对肺癌的因果作用。他们控制了许多重要背景变量，可仍可能有未测混杂把观测效应带偏。为把「这是因果」的证据加强一点，他们也报告了吸烟对车祸的效应——接近零，正是生物学事先预期的。于是，即便分析里没法排除未测混杂，基于阴性结局的这份补充分析，仍让吸烟导致肺癌的证据更硬一点。
例 16.2 Imbens and Rubin (2015) 建议把滞后结果当作阴性结局。多数情形里，有理由相信滞后结果与当前结果有相似的混杂结构。滞后结果发生在处理之前，对它的平均因果效应必须是 0。可这套建议要小心用：多数研究里，我们干脆把滞后结果当作观测到的混杂变量。
某种意义上，第 11 章的协变量平衡检查，是使用阴性对照的一个特例。与「把滞后结果当阴性对照」类似的麻烦是：那些协变量通常已经写进无混杂假定里了。因此，平衡检查失败，真正证伪的往往不是无混杂本身，而是倾向得分模型的设定。
例 16.3 对老年人的观察性研究显示：接种流感疫苗之后，下一季因肺炎/流感住院以及全因死亡的风险会显著下降——即便已经调整了测到的协变量。Jackson et al. (2006) 对这么大的效应起疑，于是对阴性结局做了补充分析。疫苗接种多在秋天开始，流感传播却常常要到冬天才明显。按生物学，疫苗的效应应当在流感季最显著。可 Jackson et al. (2006) 发现，流感季之前的效应更大——提示观测到的效应来自未测混杂。
Jackson et al. (2006) 看起来最有说服力：流感季前与流感季中、与流感相关的结局，应当有相似的混杂模式。Cornfield et al. (1959) 的额外证据弱一些：车祸与肺癌，就吸烟而言，因果机制差得很远。其实 Fisher (1957) 的批评是：吸烟与肺癌的关系，可能来自一个未观测的遗传因素（见第 17 章）。这种遗传因素也许会同时影响吸烟与肺癌，却不太像也会影响车祸。
Lipsitch et al. (2010) 是近年谈阴性结局的一篇。Rosenbaum (1989) 讨论过「已知效应」在因果推断里的角色。
16.2.2 用阴性暴露
阴性暴露是阴性结局的对偶。假定 \(Z^{\mathrm{n}}\) 是一个与 \(Z\) 相似的处理变量，并与 \(Z\) 分享同一套混杂结构。若我们相信 \(Z\perp\!\!\!\perp Y(z)\mid X\)，往往也会相信 \(Z^{\mathrm{n}}\perp\!\!\!\perp Y(z)\mid X\)。而且，我们事先知道 \(Z^{\mathrm{n}}\) 对 \(Y\) 的效应
\[ \tau(Z^{\mathrm{n}}\to Y)=\mathrm{E}\{Y(1^{\mathrm{n}})-Y(0^{\mathrm{n}})\}. \]
一个要紧的例子是 \(\tau(Z^{\mathrm{n}}\to Y)=0\)。满足这些要求的因果图如下：
Zⁿ ← X
↙ ↘
Z → Y例 16.4 Sanderson et al. (2017) 给了许多阴性暴露的例子：要判断宫内暴露对日后结局的效应，就把母亲孕期暴露与关心结局的关联，拿去和父亲暴露与同一结局的关联对照。他们综述了母亲与父亲吸烟对子代结局的研究，以及母亲与父亲 BMI 对子代日后 BMI、自闭症谱系障碍的研究。在这些例子里，我们预期母亲暴露与结局的关联，会大于父亲暴露与结局的关联。
16.2.3 小结
若不加额外假定，无混杂从根本上不可检验。观察性研究里的阴性结局与阴性对照，并不能证明或证伪无混杂；可把它们放进补充分析，能把「这是因果」的证据加强一点。可这类补充分析往往并不轻松：它需要更多数据，更要紧的是，需要对因果问题有更深的理解，才能找到说得通的阴性结局与阴性对照。
16.3 过度调整的麻烦
前面许多方法，都在无混杂假定
\[ Z\perp\!\!\!\perp\{Y(1),Y(0)\}\mid X \]
下去估因果效应。这是一条条件于 \(X\) 的假定。选对那一套能保证条件独立的 \(X\)，因此很要紧。Rosenbaum (2002b) 写道：「没有理由回避对处理前就已描述研究对象的变量做调整。」类似地，Rubin (2007) 写道：「通常，假定越是条件化的，就越容易被接受。」两人都主张：应当控制所有观测到的预处理协变量。VanderWeele and Shpitser (2011) 称之为预处理准则（pretreatment criterion）。Pearl 不同意这条建议，并给出下面两个反例。
16.3.1 M 偏倚
M 偏倚（M-bias）出现在下面这种带 M 结构的因果图里：
U₁ U₂
↓ ↘ ↙ ↓
Z X Y从这张图可以读出数据生成过程：
\[ \begin{cases} U_1\perp\!\!\!\perp U_2,\\ X=f_X(U_1,U_2,\varepsilon_X),\\ Z=f_Z(U_1,\varepsilon_Z),\\ Y=Y(z)=f_Y(U_2,\varepsilon_Y), \end{cases} \]
其中 \((\varepsilon_X,\varepsilon_Z,\varepsilon_Y)\) 是独立的随机误差。这张图里，\(X\) 被观测到，\(U_1\) 与 \(U_2\) 没有。若我们改 \(Z\) 的取值，\(Y\) 的取值完全不动。所以 \(Z\) 对 \(Y\) 的真因果效应必须是 0。从数据生成方程可以读出 \(Z\perp\!\!\!\perp Y\)，因而 \(Z\) 与 \(Y\) 的关联是 0，尤其
\[ \tau_{\mathrm{PF}}=\mathrm{E}(Y\mid Z=1)-\mathrm{E}(Y\mid Z=0)=0. \]
也就是说：不调整协变量 \(X\)，简单估计量对真参数无偏。
可一旦条件于 \(X\)，就有 \(U_1\not\perp\!\!\!\perp U_2\mid X\)（\(X\) 是 \(U_1\) 与 \(U_2\) 的对撞因子），从而 \(Z\not\perp\!\!\!\perp Y\mid X\)，并且一般地
\[ \int\bigl\{\mathrm{E}(Y\mid Z=1,X=x)-\mathrm{E}(Y\mid Z=0,X=x)\bigr\}f(x)\,dx \ne 0. \]
为拿一点直觉，考虑正态线性模型：1
\[ \begin{cases} X=aU_1+bU_2+\varepsilon_X,\\ Z=cU_1+\varepsilon_Z,\\ Y=Y(z)=dU_2+\varepsilon_Y, \end{cases} \]
其中 \((U_1,U_2,\varepsilon_X,\varepsilon_Z,\varepsilon_Y)\stackrel{\mathrm{IID}}{\sim}\mathcal{N}(0,1)\)。我们有
\[ \operatorname{cov}(Z,Y)=\operatorname{cov}(cU_1+\varepsilon_Z,\,dU_2+\varepsilon_Y)=0, \]
可按习题 1.3 的结果，给定 \(X\) 时 \(Z\) 与 \(Y\) 的偏相关系数是2
\[ \begin{aligned} \rho_{ZY\mid X} &= \frac{\rho_{ZY}-\rho_{ZX}\rho_{YX}} {\sqrt{1-\rho_{ZX}^2}\sqrt{1-\rho_{YX}^2}}\\ &\propto -\rho_{ZX}\rho_{YX}\\ &\propto -\operatorname{cov}(Z,X)\operatorname{cov}(Y,X)\\ &= -abcd, \end{aligned} \]
也就是从 \(Z\) 到 \(Y\) 那条路上系数的乘积。于是未调整估计量无偏，调整后的估计量却有正比于 \(abcd\) 的偏倚。
下面这个简单例子说明 M 偏倚。
## M bias with large sample size
n = 10^6
U1 = rnorm(n)
U2 = rnorm(n)
X = U1 + U2 + rnorm(n)
Y = U2 + rnorm(n)
## with a continuous treatment Z
Z = U1 + rnorm(n)
round(summary(lm(Y ~ Z))$coef[2, 1], 3)
## [1] 0
round(summary(lm(Y ~ Z + X))$coef[2, 1], 3)
## [1] -0.2
## with a binary treatment Z
Z = (Z >= 0)
round(summary(lm(Y ~ Z))$coef[2, 1], 3)
## [1] 0.002
round(summary(lm(Y ~ Z + X))$coef[2, 1], 3)
## [1] -0.4216.3.2 Z 偏倚
考虑下面的因果图：
U
↙b ↘c
X ─a→ Z ─τ→ Y数据生成过程为3
\[ \begin{cases} Z=aX+bU+\varepsilon_Z,\\ Y(z)=\tau z+cU+\varepsilon_Y, \end{cases} \]
其中 \((U,X,\varepsilon_Z,\varepsilon_Y)\) 独立同分布于 \(\mathcal{N}(0,1)\)。在这套生成过程里，\(X\perp\!\!\!\perp U\)，\(X\not\perp\!\!\!\perp Z\)，且 \(X\) 只通过 \(Z\) 影响 \(Y\)。
未调整估计量是
\[ \begin{aligned} \tau_{\mathrm{unadj}} &= \frac{\operatorname{cov}(Z,Y)}{\operatorname{var}(Z)}\\ &= \frac{\operatorname{cov}(Z,\,\tau Z+cU)}{\operatorname{var}(Z)}\\ &= \tau+\frac{c\,\operatorname{cov}(aX+bU,\,U)}{\operatorname{var}(Z)}\\ &= \tau+\frac{cb}{a^2+b^2+1}, \end{aligned} \]
偏倚为 \(bc/(a^2+b^2+1)\)。把 \(Y\) 对 \((Z,X)\) 做 OLS，得到的调整估计量满足
\[ \begin{cases} \mathrm{E}\bigl\{Z(Y-\tau_{\mathrm{adj}}Z-\alpha X)\bigr\}=0,\\ \mathrm{E}\bigl\{X(Y-\tau_{\mathrm{adj}}Z-\alpha X)\bigr\}=0. \end{cases} \]
解这个关于 \((\tau_{\mathrm{adj}},\alpha)\) 的线性方程组，得到
\[ \tau_{\mathrm{adj}}=\tau+\frac{bc}{b^2+1}, \tag{16.1} \]
偏倚为 \(bc/(b^2+1)\)。解方程组的细节留给习题 16.1。
于是未调整估计量的偏倚，比调整估计量更小。更有意思的是：\(X\) 与 \(Z\) 的关联越强（用 \(a\) 度量），调整估计量的偏倚反而越大。
数学推导并不特别难。可这种偏倚看起来相当神秘。直觉在这里。处理是 \(X\)、\(U\) 与其他随机误差的函数。一旦条件于 \(X\)，它差不多就只是 \(U\) 与其他随机误差的函数。因此，条件于 \(X\) 让 \(Z\) 更不随机，更要命的是，让未测混杂 \(U\) 在 \(Z\) 里扮演更重的角色。于是 \(U\) 带来的混杂偏倚，会被条件于 \(X\) 放大。这个理想化的例子说明：对某些协变量过度调整，是危险的。
Heckman and Navarro-Lozano (2004) 在模拟里观察到过这种现象。Wooldridge (2016，技术报告在 2006) 在线性模型里核实过。Pearl (2010a, 2011) 用因果图解释过。Ding et al. (2017b) 给了更一般的理论，以及一些直觉。这种偏倚叫做 Z 偏倚（Z-bias）：Pearl 原先的论文里，用符号 \(Z\) 表示我们这里的变量 \(X\)。可全书里，\(Z\) 都用来表示处理。本书第五部分会说：若某个变量满足本小节这种因果图，它就是工具变量。这也解释了为什么这种偏倚又名工具变量偏倚（instrumental variable bias）。
下面这个简单例子说明 Z 偏倚。
## Z bias with large sample size
n = 10^6
X = rnorm(n)
U = rnorm(n)
Z = X + U + rnorm(n)
Y = U + rnorm(n)
round(summary(lm(Y ~ Z))$coef[2, 1], 3)
## [1] 0.333
round(summary(lm(Y ~ Z + X))$coef[2, 1], 3)
## [1] 0.501
## stronger association between X and Z
Z = 2*X + U + rnorm(n)
round(summary(lm(Y ~ Z))$coef[2, 1], 3)
## [1] 0.167
round(summary(lm(Y ~ Z + X))$coef[2, 1], 3)
## [1] 0.501
## even stronger association between X and Z
Z = 10*X + U + rnorm(n)
round(summary(lm(Y ~ Z))$coef[2, 1], 3)
## [1] 0.01
round(summary(lm(Y ~ Z + X))$coef[2, 1], 3)
## [1] 0.516.3.3 观察性研究里该调整哪些协变量？
我们永远不知道真正的数据生成过程，而它可能相当复杂。可下面这张因果图，能把许多想法说清楚。它已经排除了 16.3.1 节那种 M 偏倚的可能。
X_R
X_Z X X_Y
↘ ↙ ↘ ↙
Z → Y
↘ ↙
X_I上面这些协变量各有各的脾气：
- \(X\) 同时影响处理与结果。条件于 \(X\) 能保证无混杂，所以应当控制 \(X\)。
- \(X_R\) 是纯噪声，既不影响处理也不影响结果。把它放进分析，不会把估计带偏，却会在有限样本里引入不必需的变异。
- \(X_Z\) 是工具变量：它只通过处理影响结果。在上面这张图里，把它放进分析不会带偏，却会增大估计的变异。可一旦有未测混杂，把它放进分析就会放大偏倚，如 16.3.2 节所示。
- \(X_Y\) 只影响结果、不影响处理。不条件于它，无混杂仍然成立。可它们能预测结果，放进分析往往能提高精度。
- \(X_I\) 受处理与结果影响。它是处理后变量，不是预处理协变量。若目标是推断处理对结果的效应，就不该把它放进来。处理后变量在因果推断里的麻烦，第六部分再谈。
若我们相信上面这张因果图，至少该调整 \(X\) 以去掉偏倚；更理想一点，再调整 \(X_Y\) 以减小方差。
16.4 习题
16.1 Z 偏倚公式的更多细节
核实 (16.1)。
16.2 Cochran 公式，或遗漏变量偏倚公式
David Cox 爵士把下面的结果叫做 Cochran 公式（Cochran, 1938; Cox, 2007）；计量学家叫做遗漏变量偏倚公式（Angrist and Pischke, 2008）。一个特例出现在 Fisher (1925)。它也是附录 B.3 里 Frisch–Waugh–Lovell 定理的对偶。
公式有两个版本。下面所有向量都是列向量。
- （总体版本） 假定 \((y_i,x_{1i},x_{2i})_{i=1}^{n}\) 独立同分布，其中 \(y_i\) 是标量，\(x_{1i}\) 是 \(K\) 维，\(x_{2i}\) 是 \(L\) 维。
我们有如下随机变量的 OLS 分解：
\[ y_i=\beta_1^\top x_{1i}+\beta_2^\top x_{2i}+\varepsilon_i, \tag{16.2} \]
\[ y_i=\gamma^\top x_{1i}+e_i, \tag{16.3} \]
\[ x_{2i}=\delta^\top x_{1i}+v_i. \tag{16.4} \]
(16.2) 叫长回归，(16.3) 叫短回归。(16.4) 里 \(\delta\) 是矩阵，因为它是向量对向量的回归。可以把 (16.4) 看成：把 \(x_{2i}\) 的每个分量分别对 \(x_{1i}\) 回归。
证明 \(\gamma=\beta_1+\delta\beta_2\)。
- （样本版本） 有一个 \(n\times 1\) 向量 \(Y\)，一个 \(n\times K\) 矩阵 \(X_1\)，一个 \(n\times L\) 矩阵 \(X_2\)。不假定任何随机性。下面全是纯线性代数。
我们可以得到如下 OLS 拟合：
\[ Y=X_1\hat\beta_1+X_2\hat\beta_2+\hat\varepsilon, \]
\[ Y=X_1\hat\gamma+\hat e, \]
\[ X_2=X_1\hat\delta+\hat v, \]
其中 \(\hat\varepsilon,\hat e,\hat v\) 是残差。最后一次 OLS 指的是：把 \(X_2\) 的每一列对 \(X_1\) 做 OLS，因而残差 \(\hat v\) 是 \(n\times L\) 矩阵。
证明 \(\hat\gamma=\hat\beta_1+\hat\delta\hat\beta_2\)。
注： 乘积项 \(\delta\beta_2\) 与 \(\hat\delta\hat\beta_2\) 常被分别称为总体层面与样本层面的遗漏变量偏倚。
16.3 推荐阅读
Imbens (2020) 综述并比较了潜在结果与因果图在因果推断里各自扮演的角色。