第 18 章　有未测混杂时，平均因果效应的敏感性分析

Cornfield 型敏感性分析，对二值结果、风险比、再条件于观测协变量，最是顺手。Ding and VanderWeele (2016) 也给过平均因果效应的 Cornfield 型方法，可不够一般，用起来也不方便。下面换一条更直接的路：盯住潜在结果的条件期望。好处是，平均因果效应那些常用估计量，都能收进同一套敏感性分析框架里。想法早出现在 Robins (1999) 与 Scharfstein et al. (1999)。这一章采用 Lu and Ding (2023) 最近的表述。

这条路，跟给平均潜在结果推最坏情形界，其实是亲戚。我先把「界」这个更简单的想法过一遍，再把它伸到敏感性分析。

18.1　引子

回想观察性研究的典范设定：\(\{Z_i,X_i,Y_i(1),Y_i(0)\}_{i=1}^{n}\stackrel{\mathrm{IID}}{\sim}\{Z,X,Y(1),Y(0)\}\)，关心平均因果效应

\[ \tau=\mathrm{E}\{Y(1)-Y(0)\}. \]

它拆成

\[ \begin{aligned} \tau &= \bigl[ \mathrm{E}(Y\mid Z=1)\operatorname{pr}(Z=1) + \mathrm{E}\{Y(1)\mid Z=0\}\operatorname{pr}(Z=0) \bigr] \\ &\quad - \bigl[ \mathrm{E}\{Y(0)\mid Z=1\}\operatorname{pr}(Z=1) + \mathrm{E}(Y\mid Z=0)\operatorname{pr}(Z=0) \bigr]. \end{aligned} \]

于是真正难的，是去估那两个反事实均值

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

估它们，大体有两种极端策略。

第一种，第三部分已经用过：靠可忽略性。假定

\[ \mathrm{E}\{Y(1)\mid Z=1,X\} = \mathrm{E}\{Y(1)\mid Z=0,X\}, \]

\[ \mathrm{E}\{Y(0)\mid Z=1,X\} = \mathrm{E}\{Y(0)\mid Z=0,X\}, \]

反事实均值就能用可观测量写出来：

\[ \mathrm{E}\{Y(1)\mid Z=0\} = \mathrm{E}\bigl\{\mathrm{E}(Y\mid Z=1,X)\mid Z=0\bigr\}, \]

同理

\[ \mathrm{E}\{Y(0)\mid Z=1\} = \mathrm{E}\bigl\{\mathrm{E}(Y\mid Z=0,X)\mid Z=1\bigr\}. \]

第二种在下一节：除了结果夹在 \(\underline{y}\) 与 \(\overline{y}\) 之间，什么都不假定。二值结果取 \(\underline{y}=0\)、\(\overline{y}=1\)，这很自然。有了有界，两个反事实均值也夹在 \(\underline{y}\) 与 \(\overline{y}\) 之间，于是 \(\tau\) 有最坏情形界。表 18.1 画出基本想法；第 18.2 节把这条策略说细一些。

表 18.1　有界结果 \([\underline{y},\overline{y}]\) 的科学表，其中 \(\underline{y}\) 与 \(\overline{y}\) 是两个常数

\(Z\)	\(Y(1)\)	\(Y(0)\)	下界 \(Y(1)\)	上界 \(Y(1)\)	下界 \(Y(0)\)	上界 \(Y(0)\)
\(1\)	\(Y_1(1)\)	\(?\)	\(Y_1(1)\)	\(Y_1(1)\)	\(\underline{y}\)	\(\overline{y}\)
\(\vdots\)	\(\vdots\)	\(\vdots\)	\(\vdots\)	\(\vdots\)	\(\vdots\)	\(\vdots\)
\(1\)	\(Y_{n_1}(1)\)	\(?\)	\(Y_{n_1}(1)\)	\(Y_{n_1}(1)\)	\(\underline{y}\)	\(\overline{y}\)
\(0\)	\(?\)	\(Y_{n_1+1}(0)\)	\(\underline{y}\)	\(\overline{y}\)	\(Y_{n_1+1}(0)\)	\(Y_{n_1+1}(0)\)
\(\vdots\)	\(\vdots\)	\(\vdots\)	\(\vdots\)	\(\vdots\)	\(\vdots\)	\(\vdots\)
\(0\)	\(?\)	\(Y_n(0)\)	\(\underline{y}\)	\(\overline{y}\)	\(Y_n(0)\)	\(Y_n(0)\)

18.2　几乎不加假定时，平均因果效应的 Manski 型最坏情形界

假定结果夹在 \(\underline{y}\) 与 \(\overline{y}\) 之间。从分解

\[ \mathrm{E}\{Y(1)\} = \mathrm{E}\{Y(1)\mid Z=1\}\operatorname{pr}(Z=1) + \mathrm{E}\{Y(1)\mid Z=0\}\operatorname{pr}(Z=0) \]

出发，\(\mathrm{E}\{Y(1)\}\) 的下界是

\[ \mathrm{E}\{Y\mid Z=1\}\operatorname{pr}(Z=1) + \underline{y}\,\operatorname{pr}(Z=0), \]

上界是

\[ \mathrm{E}\{Y\mid Z=1\}\operatorname{pr}(Z=1) + \overline{y}\,\operatorname{pr}(Z=0). \]

同理，从分解

\[ \mathrm{E}\{Y(0)\} = \mathrm{E}\{Y(0)\mid Z=1\}\operatorname{pr}(Z=1) + \mathrm{E}\{Y(0)\mid Z=0\}\operatorname{pr}(Z=0) \]

出发，\(\mathrm{E}\{Y(0)\}\) 的下界是

\[ \underline{y}\,\operatorname{pr}(Z=1) + \mathrm{E}\{Y\mid Z=0\}\operatorname{pr}(Z=0), \]

上界是

\[ \overline{y}\,\operatorname{pr}(Z=1) + \mathrm{E}\{Y\mid Z=0\}\operatorname{pr}(Z=0). \]

把这些界合起来，平均因果效应 \(\tau=\mathrm{E}\{Y(1)\}-\mathrm{E}\{Y(0)\}\) 的下界是

\[ \mathrm{E}\{Y\mid Z=1\}\operatorname{pr}(Z=1) + \underline{y}\,\operatorname{pr}(Z=0) - \overline{y}\,\operatorname{pr}(Z=1) - \mathrm{E}\{Y\mid Z=0\}\operatorname{pr}(Z=0), \]

上界是

\[ \mathrm{E}\{Y\mid Z=1\}\operatorname{pr}(Z=1) + \overline{y}\,\operatorname{pr}(Z=0) - \underline{y}\,\operatorname{pr}(Z=1) - \mathrm{E}\{Y\mid Z=0\}\operatorname{pr}(Z=0). \]

区间长度是 \(\overline{y}-\underline{y}\)。信息量不大，可仍好过先验界 \([\underline{y}-\overline{y},\ \overline{y}-\underline{y}]\)，那条的长度是 \(2(\overline{y}-\underline{y})\)。若不再加假定，观测数据的分布并不能唯一钉住 \(\tau\)。这时我们说 \(\tau\) 是部分识别的，正式定义如下。

定义 18.1（部分识别）　若观测数据的分布，与参数 \(\theta\) 的多个取值都相容，就说 \(\theta\) 是部分识别的（partially identified）。

对照定义 10.1 与定义 18.1。若 \(\theta\) 被观测数据分布唯一确定，它就可识别；否则就只是部分可识别。因此：有可忽略性时 \(\tau\) 可识别；没有可忽略性时，它只是部分可识别。

Cochran (1953) 在有缺失的调查里用过最坏情形界，后来又放下了——结果往往太保守。同样，上面 \(\tau\) 的最坏情形界，实践里常常没什么意思，因为区间经常盖住 \(0\)。而且，结果一旦无界，这条策略就用不上。

Manski 把这个想法带到因果推断（Manski, 1990），也带到许多别的计量模型（Manski, 2003）。用尽量少的假定去给因果参数划界，本身就有力量；若再配上别的定性假定，就更有用。Manski (2003) 综述过许多策略。例如，我们也许相信处理不会伤害任何单元，于是有单调性（monotonicity）：\(Y(1)\ge Y(0)\)。这时 \(\tau\) 的下界变成 \(0\)，上界不变。另一类假定是 \(Z=I\{Y(1)\ge Y(0)\}\)，也就是：人选处理，依据的是两份潜在结果之差。它也能把 \(\tau\) 的界收窄。更细的讨论超出本书范围。

18.3　平均因果效应的敏感性分析

第一种策略偏乐观：给定观测协变量后，潜在结果在处理组与对照组之间没有差别。第二种偏悲观：反事实均值，完全不从观测数据里推断。下面这条，走中间。

18.3.1　识别公式

定义

\[ \frac{\mathrm{E}\{Y(1)\mid Z=1,X\}}{\mathrm{E}\{Y(1)\mid Z=0,X\}} =\varepsilon_1(X), \qquad \frac{\mathrm{E}\{Y(0)\mid Z=1,X\}}{\mathrm{E}\{Y(0)\mid Z=0,X\}} =\varepsilon_0(X), \]

它们就是敏感性参数（sensitivity parameters）。为简单起见，可以再假定它们是不依赖 \(X\) 的常数。实践里，需要把它们固定，或让它们在事先指定的范围内变动。回想 \(\mu_1(X)=\mathrm{E}(Y\mid Z=1,X)\)、\(\mu_0(X)=\mathrm{E}(Y\mid Z=0,X)\)，分别是处理下与对照下、观测结果的条件均值。两个反事实均值，以及平均因果效应，可以如下识别。

定理 18.1　若 \(\varepsilon_1(X)\) 与 \(\varepsilon_0(X)\) 已知，则

\[ \mathrm{E}\{Y(1)\mid Z=0\} = \mathrm{E}\bigl\{\mu_1(X)/\varepsilon_1(X)\mid Z=0\bigr\}, \]

\[ \mathrm{E}\{Y(0)\mid Z=1\} = \mathrm{E}\bigl\{\mu_0(X)\varepsilon_0(X)\mid Z=1\bigr\}, \]

因而

\[\begin{align} \tau &= \mathrm{E}\bigl\{ZY+(1-Z)\mu_1(X)/\varepsilon_1(X)\bigr\} - \mathrm{E}\bigl\{Z\mu_0(X)\varepsilon_0(X)+(1-Z)Y\bigr\} \tag{18.1} \\ &= \mathrm{E}\bigl\{Z\mu_1(X)+(1-Z)\mu_1(X)/\varepsilon_1(X)\bigr\} - \mathrm{E}\bigl\{Z\mu_0(X)\varepsilon_0(X)+(1-Z)\mu_0(X)\bigr\}. \tag{18.2} \end{align}\]

定理 18.1 的证明留给习题 18.1。有了拟合的结果模型，(18.1) 与 (18.2) 便启发 \(\tau\) 的预测型与投影型估计量：

\[ \begin{aligned} \hat\tau^{\mathrm{pred}} &= \Biggl\{ n^{-1}\sum_{i=1}^{n}Z_i Y_i + n^{-1}\sum_{i=1}^{n}(1-Z_i)\hat\mu_1(X_i)/\varepsilon_1(X_i) \Biggr\} \\ &\quad - \Biggl\{ n^{-1}\sum_{i=1}^{n}Z_i\hat\mu_0(X_i)\varepsilon_0(X_i) + n^{-1}\sum_{i=1}^{n}(1-Z_i)Y_i \Biggr\}, \end{aligned} \]

以及

\[ \begin{aligned} \hat\tau^{\mathrm{proj}} &= \Biggl\{ n^{-1}\sum_{i=1}^{n}Z_i\hat\mu_1(X_i) + n^{-1}\sum_{i=1}^{n}(1-Z_i)\hat\mu_1(X_i)/\varepsilon_1(X_i) \Biggr\} \\ &\quad - \Biggl\{ n^{-1}\sum_{i=1}^{n}Z_i\hat\mu_0(X_i)\varepsilon_0(X_i) + n^{-1}\sum_{i=1}^{n}(1-Z_i)\hat\mu_0(X_i) \Biggr\}. \end{aligned} \]

「预测型」与「投影型」的用语来自抽样调查文献（Firth and Bennett, 1998; Ding and Li, 2018）；也见第 6.2.2.2 节。\(\hat\tau^{\mathrm{pred}}\) 与 \(\hat\tau^{\mathrm{proj}}\) 略有不同：前者在能看见观测结果时就用观测结果，后者把观测结果也换成拟合值。

更有意思的是，\(\tau\) 也可以用逆概率加权公式来识别。

定理 18.2　若 \(\varepsilon_1(X)\) 与 \(\varepsilon_0(X)\) 已知，则

\[ \mathrm{E}\{Y(1)\} = \mathrm{E}\left\{ w_1(X)\,\frac{Z}{e(X)}\,Y \right\}, \qquad \mathrm{E}\{Y(0)\} = \mathrm{E}\left\{ w_0(X)\,\frac{1-Z}{1-e(X)}\,Y \right\}, \]

其中

\[ w_1(X)=e(X)+\{1-e(X)\}/\varepsilon_1(X), \qquad w_0(X)=e(X)\varepsilon_0(X)+1-e(X). \]

定理 18.2 的证明留给习题 18.2。它在经典 IPW 公式上多乘了两个因子 \(w_1(X)\) 与 \(w_0(X)\)，二者同时依赖倾向得分与敏感性参数。有了拟合的倾向得分，定理 18.2 便启发 \(\tau\) 的如下估计量：

\[ \begin{aligned} \hat\tau^{\mathrm{ht}} &= n^{-1}\sum_{i=1}^{n} \frac{\bigl\{\hat e(X_i)\varepsilon_1(X_i)+1-\hat e(X_i)\bigr\}Z_i Y_i}{\varepsilon_1(X_i)\hat e(X_i)} \\ &\quad - n^{-1}\sum_{i=1}^{n} \frac{\bigl\{\hat e(X_i)\varepsilon_0(X_i)+1-\hat e(X_i)\bigr\}(1-Z_i)Y_i}{1-\hat e(X_i)} \end{aligned} \]

以及

\[ \begin{aligned} \hat\tau^{\mathrm{haj}} &= \frac{ \sum_{i=1}^{n} \dfrac{\bigl\{\hat e(X_i)\varepsilon_1(X_i)+1-\hat e(X_i)\bigr\}Z_i Y_i}{\varepsilon_1(X_i)\hat e(X_i)} }{ \sum_{i=1}^{n}Z_i/\hat e(X_i) } \\ &\quad - \frac{ \sum_{i=1}^{n} \dfrac{\bigl\{\hat e(X_i)\varepsilon_0(X_i)+1-\hat e(X_i)\bigr\}(1-Z_i)Y_i}{1-\hat e(X_i)} }{ \sum_{i=1}^{n}(1-Z_i)/\bigl(1-\hat e(X_i)\bigr) }. \end{aligned} \]

更有意思的是，同时用上拟合的倾向得分与结果模型，下面这个 \(\tau\) 的估计量是双重稳健的：

\[ \hat\tau^{\mathrm{dr}} = \hat\tau^{\mathrm{ht}} - n^{-1}\sum_{i=1}^{n} \bigl\{Z_i-\hat e(X_i)\bigr\} \left\{ \frac{\hat\mu_1(X_i)}{\hat e(X_i)\varepsilon_1(X_i)} + \frac{\hat\mu_0(X_i)\varepsilon_0(X_i)}{1-\hat e(X_i)} \right\}. \]

也就是说，在 \(\varepsilon_1(X_i)\) 与 \(\varepsilon_0(X_i)\) 已知时，只要倾向得分模型或结果模型对上其中一个，\(\hat\tau^{\mathrm{dr}}\) 就对 \(\tau\) 相合。方差可以用自助法来近似。技术细节见 Lu and Ding (2023)。

当 \(\varepsilon_1(X_i)=\varepsilon_0(X_i)=1\) 时，上面这些估计量退回第三部分的预测型估计量、IPW 估计量与双重稳健估计量。

18.3.2　例子

18.3.2.1　敏感性分析的 R 函数

下面的 R 函数计算敏感性分析的点估计。

OS_est_sa = function(z, y, x, out.family = gaussian,
                     truncps = c(0, 1), e1 = 1, e0 = 1)
{
  ## fitted propensity score
  pscore = glm(z ~ x, family = binomial)$fitted.values
  pscore = pmax(truncps[1], pmin(truncps[2], pscore))
  ## fitted potential outcomes
  outcome1 = glm(y ~ x, weights = z,
                 family = out.family)$fitted.values
  outcome0 = glm(y ~ x, weights = (1 - z),
                 family = out.family)$fitted.values
  ## outcome regression estimator
  ace.reg = mean(z*y) + mean((1-z)*outcome1/e1) -
            mean(z*outcome0*e0) - mean((1-z)*y)
  ## IPW estimators
  w1 = pscore + (1 - pscore)/e1
  w0 = pscore*e0 + (1 - pscore)
  ace.ipw0 = mean(z*y*w1/pscore) -
             mean((1 - z)*y*w0/(1 - pscore))
  ace.ipw  = mean(z*y*w1/pscore)/mean(z/pscore) -
             mean((1 - z)*y*w0/(1 - pscore)) /
             mean((1 - z)/(1 - pscore))
  ## doubly robust estimator
  aug = outcome1/pscore/e1 + outcome0*e0/(1 - pscore)
  ace.dr = ace.ipw0 + mean((z - pscore)*aug)
  return(c(ace.reg, ace.ipw0, ace.ipw, ace.dr))
}

标准误的计算，我留给习题 18.3。

18.3.2.2　回到例 10.3

取

\[ \varepsilon_1(X)=\varepsilon_0(X) \in \{1/2,\ 1/1.7,\ 1/1.5,\ 1/1.3,\ 1,\ 1.3,\ 1.5,\ 1.7,\ 2\}, \]

用下面的 R 代码，得到一整面双重稳健估计：

> nhanes_bmi = read.csv("nhanes_bmi.csv")[, -1]
> z = nhanes_bmi$School_meal
> y = nhanes_bmi$BMI
> x = as.matrix(nhanes_bmi[, -c(1, 2)])
> x = scale(x)
>
> E1 = c(1/2, 1/1.7, 1/1.5, 1/1.3, 1, 1.3, 1.5, 1.7, 2)
> E0 = c(1/2, 1/1.7, 1/1.5, 1/1.3, 1, 1.3, 1.5, 1.7, 2)
> EST = outer(E1, E0)
> ll1 = length(E1)
> ll0 = length(E0)
> for(i in 1:ll1)
+   for(j in 1:ll0)
+     EST[i, j] = OS_est_sa(z, y, x, e1 = E1[i], e0 = E0[j])[4]

表 18.2 给出点估计。估计的符号，对大于 \(1\) 的敏感性参数并不敏感，对小于 \(1\) 的却相当敏感。若参加餐食计划的人倾向于有更高的 BMI（也就是 \(\varepsilon_1(X)>1\) 且 \(\varepsilon_0(X)>1\)），餐食计划对 BMI 的平均因果效应是负的。可若参加的人倾向于有更低的 BMI，这个结论就会相当敏感。

表 18.2　平均因果效应的敏感性分析（行是 \(\varepsilon_1\)，列是 \(\varepsilon_0\)）

	\(1/2\)	\(1/1.7\)	\(1/1.5\)	\(1/1.3\)	\(1\)	\(1.3\)	\(1.5\)	\(1.7\)	\(2\)
\(1/2\)	\(11.62\)	\(10.44\)	\(9.40\)	\(8.03\)	\(4.96\)	\(0.97\)	\(-1.69\)	\(-4.35\)	\(-8.34\)
\(1/1.7\)	\(9.22\)	\(8.05\)	\(7.00\)	\(5.64\)	\(2.57\)	\(-1.42\)	\(-4.08\)	\(-6.75\)	\(-10.74\)
\(1/1.5\)	\(7.63\)	\(6.45\)	\(5.41\)	\(4.05\)	\(0.97\)	\(-3.02\)	\(-5.68\)	\(-8.34\)	\(-12.33\)
\(1/1.3\)	\(6.03\)	\(4.86\)	\(3.81\)	\(2.45\)	\(-0.62\)	\(-4.61\)	\(-7.27\)	\(-9.94\)	\(-13.93\)
\(1\)	\(3.64\)	\(2.47\)	\(1.42\)	\(0.06\)	\(-3.01\)	\(-7.01\)	\(-9.67\)	\(-12.33\)	\(-16.32\)
\(1.3\)	\(1.80\)	\(0.63\)	\(-0.42\)	\(-1.78\)	\(-4.85\)	\(-8.85\)	\(-11.51\)	\(-14.17\)	\(-18.16\)
\(1.5\)	\(0.98\)	\(-0.19\)	\(-1.24\)	\(-2.60\)	\(-5.67\)	\(-9.66\)	\(-12.33\)	\(-14.99\)	\(-18.98\)
\(1.7\)	\(0.36\)	\(-0.82\)	\(-1.86\)	\(-3.23\)	\(-6.30\)	\(-10.29\)	\(-12.95\)	\(-15.61\)	\(-19.60\)
\(2\)	\(-0.35\)	\(-1.52\)	\(-2.57\)	\(-3.93\)	\(-7.00\)	\(-10.99\)	\(-13.65\)	\(-16.32\)	\(-20.31\)

18.4　习题

18.1　定理 18.1 的证明
证明定理 18.1。

18.2　定理 18.2 的证明
证明定理 18.2。

18.3　敏感性分析里的标准误
第 18.3.2 节只给了点估计。请给出相应的自助法标准误。

18.4　处理组上平均因果效应 \(\tau_{\mathrm{T}}\) 的敏感性分析
本题把第 13 章伸到有未测混杂的情形，去估

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

\(\mathrm{E}(Y\mid Z=1)\) 很容易用样本矩来估，\(\hat\mu_{\mathrm{t}1}=\sum_{i=1}^{n}Z_i Y_i\big/\sum_{i=1}^{n}Z_i\)。唯一的反事实项是 \(\mathrm{E}\{Y(0)\mid Z=1\}\)。因此只需要敏感性参数 \(\varepsilon_0(X)\)。\(\varepsilon_0(X)\) 已知时，有下面两条识别公式。

定理 18.3　若 \(\varepsilon_0(X)\) 已知，则

\[ \mathrm{E}\{Y(0)\mid Z=1\} = \mathrm{E}\bigl\{Z\mu_0(X)\varepsilon_0(X)\bigr\}/e = \mathrm{E}\left\{ e(X)\varepsilon_0(X)\,\frac{1-Z}{1-e(X)}\,Y \right\}/e, \]

其中 \(e=\operatorname{pr}(Z=1)\)。

证明定理 18.3。

注： 定理 18.3 启发用 \(\hat\tau_{\mathrm{T}}^*=\hat\mu_{\mathrm{t}1}-\hat\mu_{\mathrm{t}0}^*\) 去估 \(\tau_{\mathrm{T}}\)，其中

\[ \hat\mu_{\mathrm{t}0}^{\mathrm{reg}} = n_1^{-1}\sum_{i=1}^{n}Z_i\varepsilon_0(X_i)\hat\mu_0(X_i), \]

\[ \hat\mu_{\mathrm{t}0}^{\mathrm{ht}} = n_1^{-1}\sum_{i=1}^{n}\varepsilon_0(X_i)\hat o(X_i)(1-Z_i)Y_i, \]

\[ \hat\mu_{\mathrm{t}0}^{\mathrm{haj}} = \sum_{i=1}^{n}\varepsilon_0(X_i)\hat o(X_i)(1-Z_i)Y_i \bigg/ \sum_{i=1}^{n}\hat o(X_i)(1-Z_i), \]

而 \(\hat o(X_i)=\hat e(X_i)/\{1-\hat e(X_i)\}\) 是给定观测协变量后、处理的估计条件优势。还可以构造双重稳健估计量 \(\hat\tau_{\mathrm{T}}^{\mathrm{dr}}=\hat\mu_{\mathrm{t}1}-\hat\mu_{\mathrm{t}0}^{\mathrm{dr}}\)，其中

\[ \hat\mu_{\mathrm{t}0}^{\mathrm{dr}} = \hat\mu_{\mathrm{t}0}^{\mathrm{ht}} - n_1^{-1}\sum_{i=1}^{n} \varepsilon_0(X_i)\, \frac{\hat e(X_i)-Z_i}{1-\hat e(X_i)}\, \hat\mu_0(X_i). \]

更多细节见 Lu and Ding (2023)。

18.5　R 代码
实现习题 18.4 里的估计量。分析第 18.3.2 节用过的数据。

18.6　推荐阅读
Rosenbaum and Rubin (1983a) 与 Imbens (2003) 是敏感性分析的两篇经典论文，不过手续更复杂。