第 21 章　从实验角度看工具变量

工具变量方法（instrumental variable method）一直是计量经济学里一件利器。处理与结果之间没有无混杂时，它仍能识别因果效应。它靠的是多出来的一个变量，叫做工具变量（instrumental variable, IV），并要求这个变量满足某些条件。这些条件，第一次遇见时并不好消化。某种意义上，IV 方法像变魔术。这一章换一个不那么像魔术的角度看：从鼓励设计（encouragement design）出发。这又一次呼应 Dorn (1953) 的建议——计划观察性研究的人，应当总问自己一句：1

若这件事能用受控实验来做，你会怎么做？

IV 方法在实验里的对应物，是鼓励设计中的不依从（noncompliance; Zelen, 1979; Powers and Swinton, 1984; Holland, 1986）。

21.1　鼓励设计与不依从

考虑一场实验，单元记作 \(i=1,\ldots,n\)。令 \(Z_i\) 为分配的处理（treatment assigned），\(1\) 为处理、\(0\) 为对照。令 \(D_i\) 为接受的处理（treatment received），同样 \(1\) 为处理、\(0\) 为对照。若某些单元 \(Z_i\ne D_i\)，不依从就出现了。不依从很常见，尤其是以人为实验单元的鼓励设计：实验者没法强迫人去接受处理，只能鼓励他们去做。令 \(Y_i\) 为关心的结果。

先考虑 \(Z\) 的完全随机化，协变量 \(X\) 暂时放下。接受的处理有潜在值 \(\{D_i(1),D_i(0)\}\)，结果有潜在值 \(\{Y_i(1),Y_i(0)\}\)，都相对于处理分配的水平 \(1\) 与 \(0\)。观测值分别是 \(D_i=Z_i D_i(1)+(1-Z_i)D_i(0)\)、\(Y_i=Z_i Y_i(1)+(1-Z_i)Y_i(0)\)。

为记号简单，假定

\[ \{Z_i,D_i(1),D_i(0),Y_i(1),Y_i(0)\}_{i=1}^{n} \ \stackrel{\mathrm{IID}}{\sim}\ \{Z,D(1),D(0),Y(1),Y(0)\}, \]

不致混淆时省掉下标 \(i\)。

从 CRE 开始。

假设 21.1（随机化）　\(Z\perp\!\!\!\perp\{D(1),D(0),Y(1),Y(0)\}\)。

随机化允许识别 \(D\) 与 \(Y\) 上的平均因果效应：

\[ \tau_D=\mathrm{E}\{D(1)-D(0)\}=\mathrm{E}(D\mid Z=1)-\mathrm{E}(D\mid Z=0) \]

以及

\[ \tau_Y=\mathrm{E}\{Y(1)-Y(0)\}=\mathrm{E}(Y\mid Z=1)-\mathrm{E}(Y\mid Z=0). \]

可以用简单均值差估计量 \(\hat\tau_D\) 与 \(\hat\tau_Y\) 去估它们。

把 \(\hat\tau_Y\) 连同标准误一起报告，叫做意向性治疗（intention-to-treat, ITT）分析。它估的是处理分配对结果的效应，假设 21.1 的 CRE 为它撑腰。可它未必回答科学问题本身——接受的处理对结果的因果效应。

21.2　潜在依从状态与效应

21.2.1　非参数识别

跟着 Imbens and Angrist (1994) 与 Angrist et al. (1996)，我们按 \(U_i=\{D_i(1),D_i(0)\}\) 的联合潜在值把总体分层。\(D\) 是二值的，于是有四种组合：

\[ U_i = \begin{cases} \mathrm{a}, & D_i(1)=1\text{ 且 }D_i(0)=1;\\ \mathrm{c}, & D_i(1)=1\text{ 且 }D_i(0)=0;\\ \mathrm{d}, & D_i(1)=0\text{ 且 }D_i(0)=1;\\ \mathrm{n}, & D_i(1)=0\text{ 且 }D_i(0)=0, \end{cases} \]

其中 \(\mathrm{a}\) 是总是接受者（always taker），\(\mathrm{c}\) 是依从者（complier），\(\mathrm{d}\) 是违抗者（defier），\(\mathrm{n}\) 是从不接受者（never taker）。我们不能同时看见 \(D_i(1)\) 与 \(D_i(0)\)，因而 \(U_i\) 是单元 \(i\) 依从行为的潜在变量。

基于 \(U\)，用全概率公式把 \(Y\) 上的平均因果效应拆成四项：

\[\begin{align} \tau_Y &= \mathrm{E}\{Y(1)-Y(0)\mid U=\mathrm{a}\}\operatorname{pr}(U=\mathrm{a}) \nonumber\\ &\quad + \mathrm{E}\{Y(1)-Y(0)\mid U=\mathrm{c}\}\operatorname{pr}(U=\mathrm{c}) \nonumber\\ &\quad + \mathrm{E}\{Y(1)-Y(0)\mid U=\mathrm{d}\}\operatorname{pr}(U=\mathrm{d}) \nonumber\\ &\quad + \mathrm{E}\{Y(1)-Y(0)\mid U=\mathrm{n}\}\operatorname{pr}(U=\mathrm{n}). \tag{21.1} \end{align}\]

于是 \(\tau_Y\) 是四个潜在子组效应的加权平均。下面把这些潜在组看得更细。

假设 21.2 把 (21.1) 的第三项限制为零。

假设 21.2（单调性）　\(\operatorname{pr}(U_i=\mathrm{d})=0\)，或对所有 \(i\) 有 \(D_i(1)\ge D_i(0)\)。也就是说，没有违抗者。

当分到对照臂的人完全接触不到处理，即对所有单元 \(D_i(0)=0\) 时，这是单侧不依从（one-sided noncompliance），假设 21.2 自动成立。在随机化下，假设 21.2 有一条可检验的推论：

\[ \operatorname{pr}(D=1\mid Z=1)\ge\operatorname{pr}(D=1\mid Z=0). \tag{21.2} \]

可假设 21.2 比不等式 (21.2) 强得多。前者在个体层面限制 \(D_i(1)\) 与 \(D_i(0)\)，后者只在平均上限制。尽管如此，当可检验推论 (21.2) 成立时，我们并不能用观测数据去推翻假设 21.2。

假设 21.3 把 (21.1) 的第一项与最后一项限制为零：处理分配对结果起作用，只通过接受的处理。

假设 21.3（排除限制）　对总是接受者 \(U_i=\mathrm{a}\) 与从不接受者 \(U_i=\mathrm{n}\)，有 \(Y_i(1)=Y_i(0)\)。

等价的说法是：对所有 \(i\)，\(D_i(1)=D_i(0)\) 必须蕴含 \(Y_i(1)=Y_i(0)\)。它要求：处理分配只有在改变了接受的处理时，才改变结果。在双盲 RCT 里，2 这在生物学上说得通，因为结果只依赖真正接受的处理。也就是说，若分配没有改变接受的处理，它也就不改变结果。若分配对结果有不经过「接受的处理」的直接效应，它就会被打破。3 例如有些 RCT 并非双盲，分配可能沿着一些未知路径走到结果。

在假设 21.2 与 21.3 下，分解 (21.1) 只剩下第二项：

\[ \tau_Y=\mathrm{E}\{Y(1)-Y(0)\mid U=\mathrm{c}\}\operatorname{pr}(U=\mathrm{c}). \tag{21.3} \]

同理，可以把 \(D\) 上的平均因果效应也拆成四项：

\[\begin{align*} \tau_D &= \mathrm{E}\{D(1)-D(0)\mid U=\mathrm{a}\}\operatorname{pr}(U=\mathrm{a}) \\ &\quad + \mathrm{E}\{D(1)-D(0)\mid U=\mathrm{c}\}\operatorname{pr}(U=\mathrm{c}) \\ &\quad + \mathrm{E}\{D(1)-D(0)\mid U=\mathrm{d}\}\operatorname{pr}(U=\mathrm{d}) \\ &\quad + \mathrm{E}\{D(1)-D(0)\mid U=\mathrm{n}\}\operatorname{pr}(U=\mathrm{n}) \\ &= 0\times\operatorname{pr}(U=\mathrm{a}) + 1\times\operatorname{pr}(U=\mathrm{c}) + (-1)\times\operatorname{pr}(U=\mathrm{d}) + 0\times\operatorname{pr}(U=\mathrm{n}), \end{align*}\]

在假设 21.2 下化成

\[ \tau_D=\operatorname{pr}(U=\mathrm{c}). \tag{21.4} \]

由 (21.4)，依从者的比例 \(\pi_{\mathrm{c}}=\operatorname{pr}(U=\mathrm{c})\) 等于处理分配对 \(D\) 的平均因果效应，而这在 CRE 下可识别。虽然我们并不能根据观测数据认出每一个依从者，却可以根据 (21.4) 认出他们在整个人群里的比例。把 (21.3) 与 (21.4) 合起来，就有下面的结果。

定理 21.1　在假设 21.2–21.3 下，若 \(\tau_D\ne 0\)，则

\[ \mathrm{E}\{Y(1)-Y(0)\mid U=\mathrm{c}\} = \frac{\tau_Y}{\tau_D}. \]

跟着 Imbens and Angrist (1994) 与 Angrist et al. (1996)，我们定义一个新的因果效应。

定义 21.1（CACE 或 LATE）　定义

\[ \tau_{\mathrm{c}}=\mathrm{E}\{Y(1)-Y(0)\mid U=\mathrm{c}\} \]

为依从者平均因果效应（complier average causal effect, CACE），也叫局部平均处理效应（local average treatment effect, LATE）。它还有等价写法：

\[ \begin{aligned} \tau_{\mathrm{c}} &= \mathrm{E}\{Y(1)-Y(0)\mid D(1)=1,D(0)=0\} \\ &= \mathrm{E}\{Y(1)-Y(0)\mid D(1)>D(0)\}. \end{aligned} \]

CACE 是 \(Z\) 对 \(Y\) 在依从者——即 \(D(1)=1\)、\(D(0)=0\) 的那些人——上的平均因果效应；在单调性下，也就是 \(D(1)>D(0)\) 的那些人。根据定义 21.1，可以把定理 21.1 写成

\[ \tau_{\mathrm{c}}=\frac{\tau_Y}{\tau_D}, \]

也就是说，CACE 或 LATE 等于 \(Y\) 上平均因果效应与 \(D\) 上平均因果效应之比。再在假设 21.1 下，可以把 CACE 识别出来。

推论 21.1　在假设 21.1–21.3 下，

\[ \tau_{\mathrm{c}} = \frac{\mathrm{E}(Y\mid Z=1)-\mathrm{E}(Y\mid Z=0)}{\mathrm{E}(D\mid Z=1)-\mathrm{E}(D\mid Z=0)}. \]

因此，在 CRE、单调性与排除限制下，可以把 CACE 非参数识别为：结果均值差除以接受的处理的均值差。

21.2.2　估计

由推论 21.1，可以用一个简单的比去估 \(\tau_{\mathrm{c}}\)：

\[ \hat\tau_{\mathrm{c}}=\frac{\hat\tau_Y}{\hat\tau_D}, \]

这叫 Wald 估计量（Wald, 1940），也叫 IV 估计量。上面的讨论里，\(Z\) 充当 \(D\) 的 IV。

方差估计可以靠下面的启发式（见例 A.3）：

\[ \hat\tau_{\mathrm{c}}-\tau_{\mathrm{c}} = (\hat\tau_Y-\tau_{\mathrm{c}}\hat\tau_D)/\hat\tau_D \approx (\hat\tau_Y-\tau_{\mathrm{c}}\hat\tau_D)/\tau_D = \hat\tau_A/\tau_D, \]

其中 \(\hat\tau_A\) 是调整结果 \(A_i=Y_i-\tau_{\mathrm{c}}D_i\) 的均值差。于是 \(\hat\tau_{\mathrm{c}}\) 的渐近方差，接近 \(\hat\tau_A\) 的方差再除以 \(\tau_D^2\)。方差估计按下面几步走：

  1. 得到调整结果 \(\hat A_i=Y_i-\hat\tau_{\mathrm{c}}D_i\)（\(i=1,\ldots,n\)）；
  2. 基于调整结果，得到 Neyman 型方差估计

\[ \hat V_{\hat A} = \frac{\hat S_{\hat A}^2(1)}{n_1} + \frac{\hat S_{\hat A}^2(0)}{n_0}, \]

其中 \(\hat S_{\hat A}^2(1)\) 与 \(\hat S_{\hat A}^2(0)\) 分别是处理组与对照组里 \(\hat A_i\) 的样本方差； 3. 最终方差估计为 \(\hat V_{\hat A}/\hat\tau_D^2\)。

上述方差估计量的理由见习题 21.2。也可以用自助法去近似 \(\hat\tau_{\mathrm{c}}\) 的方差。下面的函数计算点估计 \(\hat\tau_{\mathrm{c}}\)，以及基于 \(\hat V_{\hat A}\) 与自助法的标准误。

## IV point estimator
IV_Wald = function(Z, D, Y)
{
  tau_D = mean(D[Z==1]) - mean(D[Z==0])
  tau_Y = mean(Y[Z==1]) - mean(Y[Z==0])
  CACE  = tau_Y/tau_D
  c(tau_D, tau_Y, CACE)
}

## IV se via the delta method
IV_Wald_delta = function(Z, D, Y)
{
  est       = IV_Wald(Z, D, Y)
  AdjustedY = Y - D*est[3]
  VarAdj    = var(AdjustedY[Z==1])/sum(Z) +
              var(AdjustedY[Z==0])/sum(1 - Z)
  c(est[3], sqrt(VarAdj)/abs(est[1]))
}

## IV se via the bootstrap
IV_Wald_bootstrap = function(Z, D, Y, n.boot = 200)
{
  est      = IV_Wald(Z, D, Y)
  CACEboot = replicate(n.boot, {
    id.boot = sample(1:length(Z), replace = TRUE)
    IV_Wald(Z[id.boot], D[id.boot], Y[id.boot])[3]
  })
  c(est[3], sd(CACEboot))
}

在原假设 \(\tau_{\mathrm{c}}=0\) 下，可以用 \(\hat V_Y/\hat\tau_D^2\) 去近似方差，其中 \(\hat V_Y\) 是 \(Y\) 均值差的 Neyman 型方差估计。若真正的 \(\tau_{\mathrm{c}}\) 不是零，这个方差估计并不相合。因此它适合检验，不适合估计。可它仍能帮我们看清 ITT 估计量与 Wald 估计量的关系。ITT 估计量 \(\hat\tau_Y\) 的估计标准误是 \(\sqrt{\hat V_Y}\)。Wald 估计量 \(\hat\tau_Y/\hat\tau_D\) 本质上等于 ITT 估计量再乘 \(1/\hat\tau_D>1\)：点估计的绝对值变大，同时估计标准误也按同一个因子变大。基于方差估计 \(\hat V_Y\) 与 \(\hat V_Y/\hat\tau_D^2\)，\(\tau_Y\) 与 \(\tau_{\mathrm{c}}\) 的置信区间分别是

\[ \hat\tau_Y\pm z_{1-\alpha/2}\sqrt{\hat V_Y} \]

以及

\[ \hat\tau_Y/\hat\tau_D\pm z_{1-\alpha/2}\sqrt{\hat V_Y/\hat\tau_D^2} = \bigl(\hat\tau_Y\pm z_{1-\alpha/2}\sqrt{\hat V_Y}\bigr)/\hat\tau_D, \]

其中 \(z_{1-\alpha/2}\) 是标准正态的 \(1-\alpha/2\) 上侧分位数。这两个区间给出相同的定性结论：要么都盖住零，要么都不盖。某种意义上，IV 分析与 \(Y\) 的 ITT 分析给出同样的定性信息，只是手续更绕。

21.3　协变量

21.3.1　CRE 里的协变量调整

现在考虑带协变量的完全随机化实验，并假定

\[ Z\perp\!\!\!\perp\{D(1),D(0),Y(1),Y(0),X\}. \]

有了协变量 \(X\)，可以对 \(D\) 与 \(Y\) 分别做 Lin (2013) 估计 \(\hat\tau_{D,L}\) 与 \(\hat\tau_{Y,L}\)，再得到 \(\hat\tau_{\mathrm{c},L}=\hat\tau_{Y,L}/\hat\tau_{D,L}\)。\(\hat\tau_{\mathrm{c},L}\) 的渐近方差可以用自助法来近似。下面的函数计算点估计 \(\hat\tau_{\mathrm{c},L}\) 以及基于自助法的标准误。

## covariate adjustment in IV analysis
IV_Lin = function(Z, D, Y, X)
{
  X     = scale(as.matrix(X))
  tau_D = lm(D ~ Z + X + Z*X)$coef[2]
  tau_Y = lm(Y ~ Z + X + Z*X)$coef[2]
  CACE  = tau_Y/tau_D
  c(tau_D, tau_Y, CACE)
}

## IV_adj se via the bootstrap
IV_Lin_bootstrap = function(Z, D, Y, X, n.boot = 200)
{
  X        = scale(as.matrix(X))
  est      = IV_Lin(Z, D, Y, X)
  CACEboot = replicate(n.boot, {
    id.boot = sample(1:length(Z), replace = TRUE)
    IV_Lin(Z[id.boot], D[id.boot], Y[id.boot], X[id.boot, ])[3]
  })
  c(est[3], sd(CACEboot))
}

21.3.2　条件随机化或无混杂观察性研究里的协变量

若随机化是条件成立的，即

\[ Z\perp\!\!\!\perp\{D(1),D(0),Y(1),Y(0)\}\mid X, \]

就必须调整协变量，以免带偏。分析其实也直截了当：第三部分已经有许多估计量，分别去估 \(Z\) 对 \(D\)、\(Z\) 对 \(Y\) 的效应。把它们放进比的公式 \(\hat\tau_{\mathrm{c}}=\hat\tau_Y/\hat\tau_D\)，再用自助法近似渐近方差即可。相应的估计量与方差估计我没有在正文里实现，留给习题 21.8。

21.4　弱工具变量

21.4.1　一点模拟

即便 \(\tau_D>0\)，\(\hat\tau_D\) 仍有正概率为零，因而 \(\hat\tau_{\mathrm{c}}\) 的方差是无穷（见习题 21.1）。前面正态近似给出的，并不是 \(\hat\tau_{\mathrm{c}}\) 自己的方差，而是它渐近分布的方差。这是一个细的技术点。当 \(\tau_D\) 靠近 \(0\)——也就是弱工具变量（weak IV）——比估计量 \(\hat\tau_{\mathrm{c}}=\hat\tau_Y/\hat\tau_D\) 的有限样本性质就很差。这时 \(\hat\tau_{\mathrm{c}}\) 有有限样本偏倚，渐近分布也不是正态，相应的 Wald 型置信区间覆盖也很糟。4 结果 \(Y\) 若是二值的，我们知道 \(\tau_Y\) 必在 \(-1\) 与 \(1\) 之间，可 \(\hat\tau_{\mathrm{c}}\) 并没有保证仍落在 \(-1\) 与 \(1\) 之间。

原书图 21.1(a)(b) 给出不同 \(\pi_{\mathrm{c}}\) 下 \(\hat\tau_{\mathrm{c}}\) 与 \(\hat\tau_{\mathrm{c},L}\) 的模拟直方图。数据生成过程的细节留给 R 代码。从图上看，\(\pi_{\mathrm{c}}\) 等于 \(0.2\) 与 \(0.1\) 时，两个估计量的分布都离正态很远。基于渐近正态的统计推断，这时靠不住。

原书图 21.1：强 IV（\(\pi_{\mathrm{c}}=0.5\)）时直方图还像钟形；弱 IV 时分布又尖又野，横轴可以拉到几十甚至上千。加协变量调整之后略好一些，可 \(\pi_{\mathrm{c}}=0.1\) 时仍远远谈不上正态。

21.4.2　对弱 IV 稳健的一套手续

弱 IV 怎么办？从检验的角度看，有一个简单办法。因为 \(\tau_{\mathrm{c}}=\tau_Y/\tau_D\)，下面两条原假设在 \(\tau_D>0\) 时等价：

\[ H_0:\ \tau_{\mathrm{c}}=0 \quad\Longleftrightarrow\quad H_0':\ \tau_Y=0. \]

于是只需检验 \(H_0'\)，也就是 \(Z\) 对 \(Y\) 的平均因果效应为零。这呼应了第 21.2.2 节里 ITT 分析与 IV 分析的关系。

从估计的角度看，点估计的有限样本性质虽差，我们仍可以盯住置信区间。因为 \(\tau_{\mathrm{c}}=\tau_Y/\tau_D\)，这很像经典统计里的 Fieller–Creasy 问题。下面讨论一套给 \(\tau_{\mathrm{c}}\) 构造置信区间的策略，动机来自 Fieller (1954)；见附录 A.4.2。由置信区间与假设检验的对偶（见附录 A.2.5），可以通过反转一串原假设

\[ H_0(b):\ \tau_{\mathrm{c}}=b \]

来构造 \(\tau_{\mathrm{c}}\) 的置信集。给定真值 \(\tau_{\mathrm{c}}=b\)，有

\[ \tau_Y-b\tau_D=0. \]

于是 \(H_0(b)\) 等价于：调整结果 \(A_i(b)=Y_i-b D_i\) 上的平均因果效应为零，

\[ H_0(b):\ \tau_{A(b)}=0. \]

令 \(\hat\tau_A(b)\) 为 \(\tau_{A(b)}\) 的某个泛型点估计，相应方差估计为 \(\hat V_A(b)\)。点和方差估计都随设定而变。不带协变量的 CRE 里，\(\hat\tau_A(b)\) 是 \(A_i(b)\) 的均值差，\(\hat V_A(b)\) 是 Neyman 型方差估计。带协变量的 CRE 里，\(\hat\tau_A(b)\) 是对结果 \(A_i(b)\) 的 Lin (2013) 估计，\(\hat V_A(b)\) 是把 \(Y_i-b D_i\) 对 \((Z_i,X_i,Z_i X_i)\) 做 OLS 时的 EHW 方差，再加上第 9.1 节讨论过的校正项。5 在无混杂的观察性研究里，可以用第三部分那些策略，去估 \(A_i(b)\) 上的平均因果效应及其方差。

有了 \(\hat\tau_A(b)\) 与 \(\hat V_A(b)\)，就可以对 \(H_0(b)\) 做 Wald 型检验。把检验反转，得到 \(\tau_{\mathrm{c}}\) 的置信集：

\[ \left\{ b:\ \frac{\hat\tau_A(b)^2}{\hat V_A(b)} \le z_{1-\alpha/2}^2 \right\}, \]

其中 \(z_{1-\alpha/2}\) 是标准正态的 \(1-\alpha/2\) 上侧分位数。这接近计量里的 Anderson–Rubin 型置信区间（Anderson and Rubin, 1950）。因为它与 Fieller (1954) 也有联系，我把它叫做 Fieller–Anderson–Rubin（FAR） 置信区间。IV 强时，这些弱 IV 置信区间退回渐近置信区间；IV 弱时，它们还有额外保证。实践里我建议用它们。

例 21.1　为了对 FAR 置信区间有一点直觉，看不带协变量的 CRE 这个简单情形。置信区间里的二次不等式化成

\[ \begin{aligned} (\hat\tau_Y-b\hat\tau_D)^2 &\le z_{1-\alpha/2}^2 \Bigl[ n_1^{-1}\bigl\{\hat S_Y^2(1)+b^2\hat S_D^2(1)-2b\hat S_{YD}(1)\bigr\} \\ &\qquad + n_0^{-1}\bigl\{\hat S_Y^2(0)+b^2\hat S_D^2(0)-2b\hat S_{YD}(0)\bigr\} \Bigr], \end{aligned} \]

其中 \(\{\hat S_Y^2(1),\hat S_D^2(1),\hat S_{YD}(1)\}\) 与 \(\{\hat S_Y^2(0),\hat S_D^2(0),\hat S_{YD}(0)\}\) 分别是处理组与对照组里 \(Y\) 与 \(D\) 的样本方差与协方差。置信集可以是一个闭区间、两个不相连的区间、空集，或整条实线。细的讨论留给习题 21.4。

21.4.3　实施与模拟

下面的函数可以计算一串作为 \(b\) 的函数的 \(p\) 值。FARci 不用协变量；FARciX 用协变量，并依赖第 9.2 节定义的 linestimator。

FARci = function(Z, D, Y, Lower, Upper, grid)
{
  CIrange = seq(Lower, Upper, grid)
  Pvalue  = sapply(CIrange, function(t){
    Y_t    = Y - t*D
    Tauadj = mean(Y_t[Z==1]) - mean(Y_t[Z==0])
    VarAdj = var(Y_t[Z==1])/sum(Z) +
             var(Y_t[Z==0])/sum(1 - Z)
    Tstat  = Tauadj/sqrt(VarAdj)
    (1 - pnorm(abs(Tstat)))*2
  })
  return(list(CIrange = CIrange, Pvalue = Pvalue))
}

FARciX = function(Z, D, Y, X, Lower, Upper, grid)
{
  CIrange = seq(Lower, Upper, grid)
  X       = scale(X)
  Pvalue  = sapply(CIrange, function(t){
    Y_t   = Y - t*D
    linest = linestimator(Z, Y_t, X)
    Tstat  = linest[1]/linest[3]
    (1 - pnorm(abs(Tstat)))*2
  })
  return(list(CIrange = CIrange, Pvalue = Pvalue))
}

原书图 21.2 画出不同 \(\pi_{\mathrm{c}}\) 的模拟数据里，\(p\) 值作为 \(b\) 的函数。图 21.2(a) 与 21.2(b) 来自同一数据生成过程的两次实现，细节留给 R 代码。图 21.2(a) 里，FAR 置信集都是闭区间；图 21.2(b) 里，\(\pi_{\mathrm{c}}=0.2\) 与 \(\pi_{\mathrm{c}}=0.1\) 时置信集不再是闭区间。\(\pi_{\mathrm{c}}=0.5\) 时，两次实现里置信集的形状都还稳。

原书图 21.2：左列是第一次实现，右列是第二次。横轴是 \(\tau_{\mathrm{c}}\)，纵轴是 \(p\) 值，虚线在 \(0.05\)。强 IV 时曲线像一座山，在有限区间里穿过 \(0.05\)；弱 IV 的第二次实现里，曲线下去之后不再回到虚线下方，置信集就会往一边伸开。

21.5　应用

mediation 包里有一套数据 jobs，来自求职干预研究 JOBS II：这是一项随机化田野实验，关心就业培训干预对失业工人是否有效。变量 treat 指示参与者是否被随机选进 JOBS II 培训项目，comply 指示是否真正参加了项目。一个关心的结果是 job_seek，度量求职自我效能，取值 \(1\) 到 \(5\)。协变量包括 sex、age、marital、nonwhite、educ 与 income。

> jobsdata = read.csv("jobsdata.csv")
> Z = jobsdata$treat
> D = jobsdata$comply
> Y = jobsdata$job_seek
> getX = lm(treat ~ sex + age + marital +
+           nonwhite + educ + income,
+           data = jobsdata)
> X = model.matrix(getX)[, -1]

可以用 \(\hat\tau_{\mathrm{c}}\) 去估 \(\tau_{\mathrm{c}}\)，标准误来自 \(\hat V_{\hat A}\) 或自助法。再做协变量调整得到 \(\hat\tau_{\mathrm{c},L}\)，标准误用自助法。结果如下。点估计与标准误在几种方法之间相当稳。

> ## without covariates
> res = rbind(IV_Wald_delta(Z, D, Y),
+             IV_Wald_bootstrap(Z, D, Y, n.boot = 10^3))
> ## with covariates
> res = rbind(res,
+             IV_Lin_bootstrap(Z, D, Y, X, n.boot = 10^3))
> res = cbind(res, res[, 1] - 1.96*res[, 2],
+             res[, 1] + 1.96*res[, 2])
> row.names(res) = c("delta", "bootstrap", "with covariates")
> colnames(res) = c("est", "se", "lower CI", "upper CI")
> round(res, 3)
est se lower CI upper CI
delta \(0.109\) \(0.081\) \(-0.050\) \(0.268\)
bootstrap \(0.109\) \(0.083\) \(-0.054\) \(0.271\)
with covariates \(0.118\) \(0.082\) \(-0.042\) \(0.278\)

也可以用反转检验来构造 FAR 置信集。它们与上面的区间相近。

lower CI upper CI
without covariates \(-0.050\) \(0.267\)
with covariates \(-0.047\) \(0.282\)

原书图 21.3 画出一串检验的 \(p\) 值：上方面板不调整协变量，下面调整。两条曲线都很像一座山，大约在 \(-0.05\) 与 \(0.27\) 处穿过 \(0.05\)。

21.6　怎样读 CACE

潜在结果 \(\{D(1),D(0),Y(1),Y(0)\}\) 的记号，是相对于「分配的处理 \(Z\)」这场假想干预。因而 \(\tau_{\mathrm{c}}\) 是分配的处理对结果、在依从者上的平均因果效应。好在对依从者有 \(D=Z\)，于是也可以把 \(\tau_{\mathrm{c}}\) 读成：接受的处理对结果、在依从者上的平均因果效应。科学问题，它只回答了一部分。

有些论文用另一套记号。例如 Angrist et al. (1996) 用 \(Y_i(z,d)\) 表示单元 \(i\) 在一场 \(2\times 2\) 析因实验里的潜在结果，6 因素是分配的处理 \(z\) 与接受的处理 \(d\)。Angrist (2022, 第 3.1 节) 评论过这套记号的思想史。有了它，排除限制可以写成下面的形式。

假设 21.4（排除限制）　对所有 \(i\) 有 \(Y_i(z,d)=Y_i(d)\)，也就是说，潜在结果只是 \(d\) 的函数。

看下面的因果图，假设 21.4 排除了从 \(Z\) 到 \(Y\) 的直接箭头。这时 \(Z\) 就是 \(D\) 的 IV。

          U
         ↙ ↘
    Z → D → Y

在假设 21.4 下，加长的记号 \(Y_i(z,d)\) 收成 \(Y_i(d)\)，这也解释了「排除限制」这个名字。因此对 \(d=0,1\) 有 \(Y_i(1,d)=Y_i(0,d)\)；再配上假设 21.2，就有

\[ \begin{aligned} Y_i(z=1)-Y_i(z=0) &= Y_i\bigl(1,D_i(1)\bigr)-Y_i\bigl(0,D_i(0)\bigr) \\ &= \begin{cases} 0, & U_i=\mathrm{a},\\ 0, & U_i=\mathrm{n},\\ Y_i(d=1)-Y_i(d=0), & U_i=\mathrm{c}. \end{cases} \end{aligned} \]

上面我特意标明潜在结果是相对于 \(z\)、\(d\)，还是两者，以免混在一起。此前 \(\tau_Y\) 的分解仍然成立，于是有 Imbens and Angrist (1994) 与 Angrist et al. (1996) 的下面这条结果。

回想 \(D\) 上的平均因果效应 \(\tau_D=\mathrm{E}\{D(1)-D(0)\}\)，把 \(Y\) 上的平均因果效应定义为 \(\tau_Y=\mathrm{E}\{Y(D(1))-Y(D(0))\}\)，并把依从者平均因果效应定义为

\[ \tau_{\mathrm{c}}=\mathrm{E}\{Y(d=1)-Y(d=0)\mid U=\mathrm{c}\}. \]

定理 21.2　在假设 21.2–21.4 下，

\[ Y\bigl(D(1)\bigr)-Y\bigl(D(0)\bigr) = \bigl\{D(1)-D(0)\bigr\}\times\bigl\{Y(d=1)-Y(d=0)\bigr\} \]

且 \(\tau_{\mathrm{c}}=\tau_Y/\tau_D\)。

证明与定理 21.1 几乎相同，只是记号改了一下。留给习题 21.3。从记号 \(Y_i(d)\) 出发，把 \(\tau_{\mathrm{c}}\) 读成「接受的处理对结果、在依从者上的平均因果效应」，更顺口。

21.7　习题

21.1　Wald 估计量的方差
证明 \(\operatorname{var}(\hat\tau_{\mathrm{c}})=\infty\)。

21.2　Wald 估计量的渐近方差及其估计
考虑 \(n\to\infty\) 的大样本机制。先证明 \(\sqrt{n}(\hat\tau_{\mathrm{c}}-\tau_{\mathrm{c}})\to\mathcal{N}(0,V)\) 依分布，并找出 \(V\)。再证明 \(\hat V_{\hat A}/V\to 1\) 依概率。

21.3　Imbens and Angrist (1994) 与 Angrist et al. (1996) 主定理的证明
证明定理 21.2。

21.4　关于 FAR 置信集再多说几句
例 21.1 里的置信集可以是一个闭区间、两个不相连的区间、空集，或整条实线。请给出每一种情形的精确条件。

21.5　FAR 置信集的更多模拟
原书图 21.2 给出的是不使用协变量的 FAR 置信集。请在 CRE 里做平行模拟：FAR 置信集改为调整协变量。

21.6　二值 IV 与有序的接受处理
Angrist and Imbens (1995) 讨论过更一般的设定：二值 IV \(Z\)，有序的接受处理 \(D\in\{0,1,\ldots,J\}\)，以及结果 \(Y\)。接受的处理相对于二值 IV 有潜在结果 \(D(1)\) 与 \(D(0)\)，结果相对于二值 IV 与有序接受处理有潜在结果 \(Y(z,d)\)。把第 21.6 节的讨论以及相应的 IV 假定伸开如下。

假设 21.5　(1) 随机化：\(Z\perp\!\!\!\perp\{D(z),Y(z,d):z=0,1;\ d=0,1,\ldots,J\}\)；(2) 单调性：\(D(1)\ge D(0)\)；(3) 排除限制：对所有 \(z=0,1\) 与 \(d=0,1,\ldots,J\) 有 \(Y(z,d)=Y(d)\)。

他们证明了下面的定理 21.3。

定理 21.3　在假设 21.5 下，

\[ \frac{\mathrm{E}(Y\mid Z=1)-\mathrm{E}(Y\mid Z=0)}{\mathrm{E}(D\mid Z=1)-\mathrm{E}(D\mid Z=0)} = \sum_{j=1}^{J} w_j\,\mathrm{E}\bigl\{Y(j)-Y(j-1)\mid D(1)\ge j>D(0)\bigr\}, \]

其中

\[ w_j = \frac{\operatorname{pr}\{D(1)\ge j>D(0)\}}{\sum_{j'=1}^{J}\operatorname{pr}\{D(1)\ge j'>D(0)\}}. \]

证明定理 21.3。

注： 当 \(J=1\) 时，定理 21.3 退回定理 21.2。一般的 \(J\) 下，它说的是：标准 IV 公式识别的，是某些潜在子组效应的加权平均。权重与潜在组 \(D(1)\ge j>D(0)\) 的概率成正比，潜在子组效应 \(\mathrm{E}\{Y(j)-Y(j-1)\mid D(1)\ge j>D(0)\}\) 比较的是相邻的接受处理水平。可这份加权平均并不好解释，因为这些潜在组是重叠的。

证明可以很烦。一个技巧是把处理分配 \(z\) 下接受的处理与结果写成

\[ D(z)=\sum_{j=0}^{J}j\,I\{D(z)=j\}, \qquad Y\bigl(D(z)\bigr)=\sum_{j=0}^{J}Y(j)\,I\{D(z)=j\}, \]

从而

\[ D(1)-D(0) = \sum_{j=0}^{J} j\bigl[I\{D(1)=j\}-I\{D(0)=j\}\bigr] \]

以及

\[ Y\bigl(D(1)\bigr)-Y\bigl(D(0)\bigr) = \sum_{j=0}^{J} Y(j)\bigl[I\{D(1)=j\}-I\{D(0)=j\}\bigr]. \]

然后用下面的阿贝尔引理（Abel’s lemma），也称分部求和：

\[ \sum_{j=0}^{J}f_j(g_{j+1}-g_j) = f_J g_{J+1}-f_0 g_0-\sum_{j=1}^{J}g_j(f_j-f_{j-1}), \]

对适当指定的序列 \((f_j)\) 与 \((g_j)\)。

21.7　数据分析：流感疫苗鼓励设计（McDonald et al., 1992）
fludata.txt 来自 McDonald et al. (1992) 的随机化鼓励设计，Hirano et al. (2000) 也重新分析过。它包含下列变量：

变量 含义
assign 是否被鼓励接种流感疫苗（二值）
receive 是否接种了流感疫苗（二值）
outcome 是否发生流感相关住院（二值）
age 患者年龄
sex 患者性别
race 患者种族
copd 慢性阻塞性肺病
dm 糖尿病
heartd 心脏病
renal 肾病
liverd 肝病

请在调整与不调整协变量两种情形下分析这套数据。

21.8　条件于协变量的 IV 估计
用 R 实现第 21.3.2 节提到的估计量与相应方差估计。

注： 这道题对习题 21.9 有用。

21.9　数据分析：卡罗林斯卡数据
回到习题 12.5。Rubin (2008) 用卡罗林斯卡数据当 IV 方法的例子。在 karolinska.txt 里，患者是否在高容量医院诊断，可以看成是否在高容量医院治疗的 IV。这在给定其他观测协变量时说得通。更多细节见 Rubin (2008) 的分析。

请重新分析这套数据，假定 IV 在给定观测协变量后是随机分配的。

21.10　数据分析：一项就业培训项目
文件 jobtraining.rtf 描述了数据文件 X.csv 与 Y.csv。

X.csv 含预处理协变量。抽样权重变量 wgt 也可以看成一个协变量。许多先前的分析做了这个简化，尽管在调查数据的统计分析里，这件事一直有争议。请做有协变量与无协变量两种分析。

Y.csv 含抽样权重、分配的处理、接受的处理，以及许多处理后变量。因此，这套数据里有许多结果，取决于你关心的问题。数据也有许多麻烦。第一，有些结果缺失。第二，失业的人没有工资。第三，结果随时间重复观测。分析时请说明：你选了哪些关心的问题，用了哪些估计量。

注： Schochet et al. (2008) 分析过原始数据。Frumento et al. (2012) 基于后面第 26 章的框架，给过更精细的分析。

21.11　推荐阅读
Angrist et al. (1996) 在计量的 IV 视角与基于潜在结果的统计因果推断之间搭了一座桥，并用一个应用说明它有用。

IV 的另一些早期文献包括 Permutt and Hebel (1989)、Sommer and Zeger (1991)、Baker and Lindeman (1994)，以及 Cuzick et al. (1997)。

  1. 第 10 章也引过这句话。↩︎

  2. 一般而言，把实验盲起来更好，以免安慰剂效应、患者预期等带来各种偏倚。双盲 RCT 里，医生与患者都不知道处理；单盲 RCT 里，患者不知道、医生知道。有时双盲甚至单盲都做不到，那些试验叫做开放试验。↩︎

  3. 这里的「直接效应」用得并不正式。更细的讨论见后面第 27、28 章。↩︎

  4. 理论里常常假定 \(\tau_D\) 是 \(n^{-1/2}\) 这一阶。这个机制下，依从者的比例随 \(n\) 趋于零。IV 方法识别的，便是比例缩向零的一个子组平均因果效应。这是为理论分析而造的机制，实践里很难为它辩护。下面的讨论并不假定它。↩︎

  5. 若采用有限总体视角，就不需要这项校正。↩︎

  6. 「析因实验」这个名字来自实验设计文献（Dasgupta et al., 2015）。实验者把多个因素随机分配给每个单元。析因实验里的处理有多个水平。↩︎