第 19 章　有未测混杂时，匹配观察性研究的 Rosenbaum 式 (p) 值

Rosenbaum (1987b) 给匹配观察性研究提过一套敏感性分析。它也适用于更一般的匹配（Rosenbaum, 2002b），可一对一匹配时，理论最干净。与第 17、18 章不同：Rosenbaum 这一套，最顺手的场合，是匹配观察性研究里，去检验「没有个体处理效应」这条尖锐原假设。

19.1　匹配数据的敏感性分析模型

考虑一项匹配观察性研究。用 \((i,j)\) 表示第 \(i\) 对里的第 \(j\) 个单元（\(i=1,\ldots,n\)；\(j=1,2\)）。若匹配是精确的，单元 \((i,1)\) 与 \((i,2)\) 有相同的协变量 \(X_i\)。假定 IID 抽样，并把定义 11.1 稍稍伸开，把倾向得分写成

\[ e_{ij}=\operatorname{pr}\bigl\{Z_{ij}=1\mid X_i,Y_{ij}(1),Y_{ij}(0)\bigr\}. \]

令 \(\mathbb{S}_i=\{Y_{i1}(1),Y_{i1}(0),Y_{i2}(1),Y_{i2}(0)\}\) 表示第 \(i\) 对里全部潜在结果。条件于事件 \(Z_{i1}+Z_{i2}=1\)，有

\[\begin{align*} \pi_{i1} &= \operatorname{pr}\bigl\{Z_{i1}=1\mid X_i,\mathbb{S}_i,Z_{i1}+Z_{i2}=1\bigr\} \\ &= \frac{ \operatorname{pr}\{Z_{i1}=1,Z_{i2}=0\mid X_i,\mathbb{S}_i\} }{ \operatorname{pr}\{Z_{i1}+Z_{i2}=1\mid X_i,\mathbb{S}_i\} } \\ &= \frac{ \operatorname{pr}\{Z_{i1}=1,Z_{i2}=0\mid X_i,\mathbb{S}_i\} }{ \operatorname{pr}\{Z_{i1}=1,Z_{i2}=0\mid X_i,\mathbb{S}_i\} + \operatorname{pr}\{Z_{i1}=0,Z_{i2}=1\mid X_i,\mathbb{S}_i\} } \\ &= \frac{e_{i1}(1-e_{i2})}{e_{i1}(1-e_{i2})+(1-e_{i1})e_{i2}}. \end{align*}\]

定义 \(o_{ij}=e_{ij}/(1-e_{ij})\) 为单元 \((i,j)\) 接受处理的优势，于是

\[ \pi_{i1}=\frac{o_{i1}}{o_{i1}+o_{i2}}. \]

有可忽略性时，\(e_{ij}\) 只是 \(X_i\) 的函数，因而 \(e_{i1}=e_{i2}\)、\(\pi_{i1}=1/2\)。再条件于协变量与潜在结果，处理分配机制就与处理、对照概率相等的 MPE 相同。这正是第 15.1 节分析匹配观察性研究的那条策略。

一般情形里，\(e_{ij}\) 也会是未观测潜在结果的函数，可以从 \(0\) 走到 \(1\)。Rosenbaum (1987b) 的敏感性分析模型，给优势比 \(o_{i1}/o_{i2}\) 加上界。

假设 19.1（Rosenbaum 的敏感性分析模型）　优势比有界：

\[ o_{i1}/o_{i2}\le\Gamma, \qquad o_{i2}/o_{i1}\le\Gamma, \qquad (i=1,\ldots,n), \]

其中 \(\Gamma\ge 1\) 事先指定。等价地，

\[ \frac{1}{1+\Gamma}\le\pi_{i1}\le\frac{\Gamma}{1+\Gamma}, \qquad (i=1,\ldots,n). \]

在假设 19.1 下，我们面对的是一场有偏的 MPE：各对里处理与对照的概率既不相等，也对与对之间可以不同。当 \(\Gamma=1\) 时，\(\pi_{i1}=1/2\)，回到标准 MPE。因此 \(\Gamma>1\) 度量的是：匹配时漏掉的变量，让局面离开理想 MPE 有多远。

19.2　Rosenbaum 模型下的最坏情形 \(p\) 值

考虑检验尖锐原假设

\[ H_{0\mathrm{f}}:\ Y_{ij}(1)=Y_{ij}(0) \quad\text{对 }i=1,\ldots,n\text{ 与 }j=1,2, \]

依据的是对内差 \(\hat\tau_i=(2Z_{i1}-1)(Y_{i1}-Y_{i2})\)（\(i=1,\ldots,n\)）。在 \(H_{0\mathrm{f}}\) 下，\(\lvert\hat\tau_i\rvert\) 是固定的；若 \(\hat\tau_i\ne 0\)，则 \(S_i=I(\hat\tau_i>0)\) 是随机的。考虑下面这一族检验统计量：

\[ T=\sum_{i=1}^{n}S_i q_i, \]

其中 \(q_i\ge 0\) 是 \((\lvert\hat\tau_1\rvert,\ldots,\lvert\hat\tau_n\rvert)\) 的函数。特例包括符号统计量、配对 \(t\) 统计量（差一个常数平移），以及 Wilcoxon 符号秩统计量：

\[ T=\sum_{i=1}^{n}S_i, \qquad T=\sum_{i=1}^{n}S_i\lvert\hat\tau_i\rvert, \qquad T=\sum_{i=1}^{n}S_i R_i, \]

其中 \((R_1,\ldots,R_n)\) 是 \((\lvert\hat\tau_1\rvert,\ldots,\lvert\hat\tau_n\rvert)\) 的秩。

一般的 \(\Gamma\) 下，假设 19.1 并不把各 \(\pi_{i1}\) 钉死，因而 \(T\) 在原假设下的精确分布可以很麻烦。好在我们并不需要知道精确分布，只要知道：在假设 19.1 下，能给出最大 \(p\) 值的那份最坏情形分布。它对应

\[ S_i\ \stackrel{\mathrm{IID}}{\sim}\ \mathrm{Bernoulli}\!\left(\frac{\Gamma}{1+\Gamma}\right). \]

相应的 \(T\) 有均值

\[ \mathrm{E}_{\Gamma}(T) = \frac{\Gamma}{1+\Gamma}\sum_{i=1}^{n}q_i \]

与方差

\[ \operatorname{var}_{\Gamma}(T) = \frac{\Gamma}{(1+\Gamma)^2}\sum_{i=1}^{n}q_i^2, \]

以及正态近似

\[ \frac{ T-\dfrac{\Gamma}{1+\Gamma}\displaystyle\sum_{i=1}^{n}q_i }{ \sqrt{ \dfrac{\Gamma}{(1+\Gamma)^2}\displaystyle\sum_{i=1}^{n}q_i^2 } } \to\mathcal{N}(0,1) \quad\text{依分布}. \]

实践里，可以把一串 \(p\) 值当成 \(\Gamma\) 的函数报出来。

19.3　例子

19.3.1　回到 LaLonde 数据

我们在匹配后的 LaLonde 数据上做 Rosenbaum 式敏感性分析。用 Matching 包可以构造匹配后的数据集。

library("Matching")
library("sensitivitymv")
library("sensitivitymw")
dat <- read.table("cps1re74.csv", header = TRUE)
dat$u74 <- as.numeric(dat$re74 == 0)
dat$u75 <- as.numeric(dat$re75 == 0)
y = dat$re78
z = dat$treat
x = as.matrix(dat[, c("age", "educ", "black",
                      "hispan", "married", "nodegree",
                      "re74", "re75", "u74", "u75")])
matchest = Match(Y = y, Tr = z, X = x)
ytreated = y[matchest$index.treated]
ycontrol = y[matchest$index.control]
datamatched = cbind(ytreated, ycontrol)

我们用检验统计量 \(T=\sum_{i=1}^{n}S_i\lvert\hat\tau_i\rvert\)。在理想 MPE、\(\Gamma=1\) 时，可以模拟 \(T\) 的分布，得到 \(p\) 值 \(0.002\)，见原书图 19.1 第一幅。\(\Gamma\) 稍稍升到 \(1.1\)，\(T\) 的最坏情形分布往右移，\(p\) 值升到 \(0.011\)。再升到 \(1.3\)，分布继续右移，\(p\) 值超过 \(0.05\)。原书图 19.2 给出 \(\hat\tau_i\) 的直方图，以及 \(p\) 值随 \(\Gamma\) 的变化；\(\Gamma=1.233\) 度量的是：在 \(0.05\) 水平上仍能拒绝原假设的、最大那份混杂。

也可以用 sensitivitymw 包里的 senmw，直接得到一串对着 \(\Gamma\) 的 \(p\) 值。

Gamma  = seq(1, 1.4, 0.001)
Pvalue = Gamma
for(i in 1:length(Gamma))
{
  Pvalue[i] = senmw(datamatched, gamma = Gamma[i],
                    method = "t")$pval
}

原书图 19.1 是三幅直方图：\(\Gamma=1,1.1,1.3\) 时 \(T\) 的最坏情形分布，竖线标出观测到的 \(T\)；分布一右移，观测值就显得不那么极端。原书图 19.2 左边是 \(\hat\tau_i\) 的密度直方图，多数靠近 \(0\)、右边拖着长尾；右边是 \(p\) 值随 \(\Gamma\) 上升的曲线，在 \(1.233\) 处穿过 \(0.05\)。

19.3.2　Rosenbaum 两个 R 包里的例子

erpcp 数据来自 R 包 sensitivitymw。它含 \(n=39\) 对匹配：一名焊工、一名对照，按观测协变量年龄与吸烟匹配。结果是 DNA 洗脱率（DNA elution rate）。原书图 19.3(a) 给出对内结果差的直方图，以及基于配对 \(t\) 统计量、对着 \(\Gamma\) 的 \(p\) 值。下面的 R 代码生成图 19.3(a)。

par(mfrow = c(1, 2), mai = c(0.8, 0.8, 0.3, 0.3))
data(erpcp)
hist(erpcp[, 1] - erpcp[, 2], main = "erpcp",
     xlab = expression(hat(tau)[i]),
     freq = FALSE)
Gamma  = seq(1, 5, 0.005)
Pvalue = Gamma
for(i in 1:length(Gamma))
{
  Pvalue[i] = senmw(erpcp, gamma = Gamma[i], method = "t")$pval
}
gammastar = Gamma[which(Pvalue >= 0.05)[1]]
gammastar
plot(Pvalue ~ Gamma, type = "l",
     xlab = expression(Gamma),
     ylab = "p-value")
abline(h = 0.05, lty = 2)
abline(v = gammastar, lty = 2)

lead250 数据也来自 sensitivitymw。它含 \(n=250\) 对匹配：一名每日吸烟者、一名不吸烟对照，按 NHANES 里的观测协变量性别、年龄、种族、教育水平与家庭收入匹配。结果是血铅水平，单位 \(\mu\mathrm{g}/\mathrm{l}\)。原书图 19.3(b) 给出对内结果差的直方图，以及基于配对 \(t\) 统计量、对着 \(\Gamma\) 的 \(p\) 值。下面的 R 代码生成图 19.3(b)。

par(mfrow = c(1, 2), mai = c(0.8, 0.8, 0.3, 0.3))
data(lead250)
hist(lead250[, 1] - lead250[, 2],
     main = "lead250",
     xlab = expression(hat(tau)[i]),
     freq = FALSE)
Gamma  = seq(1, 2.5, 0.001)
Pvalue = Gamma
for(i in 1:length(Gamma))
{
  Pvalue[i] = senmw(lead250, gamma = Gamma[i], method = "t")$pval
}
gammastar = Gamma[which(Pvalue >= 0.05)[1]]
gammastar
plot(Pvalue ~ Gamma, type = "l",
     xlab = expression(Gamma), ylab = "p-value")
abline(h = 0.05, lty = 2)
abline(v = gammastar, lty = 2)

原书图 19.3 两例结构相同：左边是 \(\hat\tau_i\) 的直方图，右边是 \(p\) 值随 \(\Gamma\) 的曲线，虚线标出 \(0.05\) 水平与临界 \(\Gamma\)。焊工那例，曲线穿过 \(0.05\) 时 \(\Gamma\) 大约在 \(3\) 与 \(4\) 之间；血铅那例大约在 \(2\) 附近。

19.4　习题

19.1　假设 19.1 的一个模型
假设 19.1 有下面的等价形式。

假设 19.2　倾向得分满足模型

\[ \log\frac{\pi_{ij}}{1-\pi_{ij}} = g(X_i)+\gamma U_{ij}, \qquad (i=1,\ldots,n;\ j=1,2), \]

其中 \(g(\cdot)\) 是未知函数，\(U_{ij}\in[0,1]\) 是单元 \((i,j)\) 的一个有界未观测协变量。

证明：若假设 19.2 成立，则假设 19.1 必成立，且 \(\Gamma=e^{\gamma}\)。

注： Rosenbaum (2002b, 第 108 页) 也说明：若假设 19.1 成立，则对某些 \(U_{ij}\)，假设 19.2 也成立。假设 19.2 也许更好解释，因为 \(\gamma\) 度量的是 \(U\) 对处理的条件优势比的对数；逻辑回归系数的解释见附录 B.6。

19.2　Rosenbaum 方法的应用
用基于匹配的 Rosenbaum 方法，重新分析例 10.3。

19.3　推荐阅读
Rosenbaum (2015) 给他两个做匹配观察性研究敏感性分析的 R 包写过教程。