第 25 章　工具变量方法的应用：孟德尔随机化

Katan (1986) 担心那些观察性研究：它们提示低血清胆固醇与癌症风险有关联。可我们谈过，观察性研究会受未测混杂之苦。于是，把看见的关联读成因果，并不容易。Katan (1986) 盯的这个问题里，甚至可能是癌症早期反过来造成了低血清胆固醇。用标准流行病学研究把血清胆固醇对癌症的因果效应拆开，看起来很难。Katan (1986) 的论点是：载脂蛋白 E（Apolipoprotein E）基因与血清胆固醇有关联，却并不直接改癌症状态。因此，若低血清胆固醇真的导致癌症，我们应当在「带有会改血清胆固醇的基因型」与「不带」的人之间，看见癌症风险的差别。用我们因果推断的语言说：Katan (1986) 提议把载脂蛋白 E 基因当作 IV。

Katan (1986) 并没有做任何数据分析，只是提出一套概念上的设计：它既能对付未测混杂，也能对付反向因果（reverse causality）。打那以后，多亏现代的全基因组关联研究，出现了更复杂、更精细的研究。它们用遗传信息当流行病学研究里暴露的 IV，去估暴露对结果的因果效应。这些研究都受到孟德尔第二定律（Mendel’s second law）即自由组合定律（law of random assortment）的启发：一种性状的遗传，独立于其他性状的遗传。因此，用遗传信息当 IV 的方法，叫做孟德尔随机化（Mendelian Randomization, MR）。

25.1　背景与动机

从图上看，图 25.1 画出处理 \(D\)、结果 \(Y\)、未测混杂 \(U\)，以及遗传 IV \(G_1,\ldots,G_p\) 的因果图。许多 MR 研究里，遗传 IV 是单核苷酸多态性（single-nucleotide polymorphism, SNP）。因为存在多效性（pleiotropy），1 遗传 IV 可能对关心的结果有直接效应，所以图 25.1 也允许排除限制被打破。

     G₁ ········· α₁ ·········
         \γ₁                    ↘
     G₂ ──γ₂──→ D ────β────→ Y
         ↗              ↗
        :              /
     Gₚ ──γₚ──→      U
         ········· αₚ ·········

原书图 25.1：MR 的因果图。实线是基因对 \(D\) 的效应 \(\gamma\)，以及 \(D\) 对 \(Y\) 的效应 \(\beta\)；虚线是基因对 \(Y\) 的直接效应 \(\alpha\)（多效性）；\(U\) 同时指向 \(D\) 与 \(Y\)。

标准线性 IV 模型假定：IV 对结果没有直接效应。下面的定义 25.1 同时给出结构式与约简式。

定义 25.1（线性 IV 模型）　标准线性 IV 模型

\[ \begin{align} Y &= \beta_0+\beta D+\beta_u U+\varepsilon_Y, \tag{25.1} \\ D &= \gamma_0+\gamma_1 G_1+\cdots+\gamma_p G_p+\gamma_u U+\varepsilon_D, \tag{25.2} \end{align} \]

有约简式

\[ \begin{align} Y &= (\beta_0+\beta\gamma_0)+\beta\gamma_1 G_1+\cdots+\beta\gamma_p G_p+(\beta_u+\beta\gamma_u)U+\varepsilon_Y, \tag{25.3} \\ D &= \gamma_0+\gamma_1 G_1+\cdots+\gamma_p G_p+\gamma_u U+\varepsilon_D. \tag{25.4} \end{align} \]

定义 25.2 允许排除限制被打破。这时 \(G_1,\ldots,G_p\) 不是有效 IV。

定义 25.2（可能含无效 IV 的线性模型）　线性模型

\[ \begin{align} Y &= \beta_0+\beta D+\alpha_1 G_1+\cdots+\alpha_p G_p+\beta_u U+\varepsilon_Y, \tag{25.5} \\ D &= \gamma_0+\gamma_1 G_1+\cdots+\gamma_p G_p+\gamma_u U+\varepsilon_D, \tag{25.6} \end{align} \]

有约简式

\[ \begin{align} Y &= (\beta_0+\beta\gamma_0)+(\alpha_1+\beta\gamma_1)G_1+\cdots+(\alpha_p+\beta\gamma_p)G_p \\ &\qquad +(\beta_u+\beta\gamma_u)U+\varepsilon_Y, \tag{25.7} \\ D &= \gamma_0+\gamma_1 G_1+\cdots+\gamma_p G_p+\gamma_u U+\varepsilon_D. \tag{25.8} \end{align} \]

定义 25.1 与 25.2 和第 23 章的线性 IV 模型略有不同。它们把混杂 \(U\) 显式写了出来。可这点差别并不从根本上改后面的讨论。

在定义 25.1 里，有排除限制，于是

\[ \Gamma_j=\beta\gamma_j, \qquad (j=1,\ldots,p). \]

在定义 25.2 里，没有排除限制，于是

\[ \Gamma_j=\alpha_j+\beta\gamma_j, \qquad (j=1,\ldots,p). \]

若有个体数据，就可以在定义 25.1 的线性 IV 模型下，用经典 TSLS 去估 \(\beta\)。可大多数 MR 研究并没有个体数据，而是来自多项全基因组关联研究的汇总统计量（summary statistics）。一套典型的、基于汇总统计量的 MR 研究，有下面的数据结构。

假设 25.1（基于汇总统计量的 MR 研究）　(a) 我们有处理对遗传 IV 的回归系数 \(\hat\gamma_1,\ldots,\hat\gamma_p\)，以及标准误 \(\mathrm{se}_{D1},\ldots,\mathrm{se}_{Dp}\)。假定依概率有

\[ \hat\gamma_1\to\gamma_1,\ \ldots,\ \hat\gamma_p\to\gamma_p, \tag{25.9} \]

并忽略标准误里的不确定性。

  1. 我们有结果对遗传 IV 的回归系数 \(\hat\Gamma_1,\ldots,\hat\Gamma_p\)，以及标准误 \(\mathrm{se}_{Y1},\ldots,\mathrm{se}_{Yp}\)。假定依概率有

\[ \hat\Gamma_1\to\Gamma_1,\ \ldots,\ \hat\Gamma_p\to\Gamma_p, \tag{25.10} \]

并忽略标准误里的不确定性。

  1. 假定 \(\hat\gamma_1,\ldots,\hat\gamma_p,\hat\Gamma_1,\ldots,\hat\Gamma_p\) 联合正态，并且相互独立。

回归系数的渐近正态，可以由中心极限定理撑腰。大样本里，标准误是真标准误的准确估计。因此，真正微妙的假定，是回归系数的联合独立性。\(\hat\gamma_j\) 与 \(\hat\Gamma_j\) 之间的独立说得通，因为它们常常来自不同样本。

\(\hat\gamma_j\) 之间的独立，在 \(G_j\) 相互独立、并且 \(D\) 的真线性模型带同方差误差时，也可以说得通。见附录 B.4。可若线性模型的误差是异方差的，这条假定就别扭了。没有 \(G_j\) 的独立，也很难为系数的独立辩护。若回归系数来自非线性模型，就更难。\(\hat\Gamma_j\) 之间的独立，理由类似。

这一章盯的是：基于假设 25.1 这些汇总统计量，怎样对 \(\beta\) 做统计推断。

25.2　基于汇总统计量的 MR

25.2.1　固定效应估计量

在定义 25.1 下，\(\alpha_j=0\)，从而对所有 \(j\) 都有 \(\beta=\Gamma_j/\gamma_j\)。一条简单的路，是所谓的元分析（meta-analysis; Bowden et al., 2018）：把多个估计 \(\hat\beta_j=\hat\Gamma_j/\hat\gamma_j\) 合成同一个参数 \(\beta\)，也叫固定效应。用 delta 方法（见例 A.3），比值估计量 \(\hat\beta_j\) 的近似平方标准误是

\[ \mathrm{se}_j^{2} = \bigl(\mathrm{se}_{Yj}^{2}+\hat\beta_j^{2}\mathrm{se}_{Dj}^{2}\bigr)/\hat\gamma_j^{2}. \]

因此，估 \(\beta\) 的最佳线性组合，是按方差倒数做的 Fisher 加权（Fisher weighting；见习题 A.6）：

\[ \hat\beta_{\mathrm{fisher0}} = \frac{\sum_{j=1}^{p}\hat\beta_j/\mathrm{se}_j^{2}}{\sum_{j=1}^{p}1/\mathrm{se}_j^{2}}, \]

其方差为 \(\bigl(\sum_{j=1}^{p}1/\mathrm{se}_j^{2}\bigr)^{-1}\)。若忽略 \(\hat\gamma_j\) 带来的、由 \(\mathrm{se}_{Dj}\) 度量的不确定性，估计量化成

\[ \hat\beta_{\mathrm{fisher1}} = \frac{\sum_{j=1}^{p}\hat\beta_j\hat\gamma_j^{2}/\mathrm{se}_{Yj}^{2}}{\sum_{j=1}^{p}\hat\gamma_j^{2}/\mathrm{se}_{Yj}^{2}} = \frac{\sum_{j=1}^{p}\hat\Gamma_j\hat\gamma_j/\mathrm{se}_{Yj}^{2}}{\sum_{j=1}^{p}\hat\gamma_j^{2}/\mathrm{se}_{Yj}^{2}}, \]

其方差为 \(\bigl(\sum_{j=1}^{p}\hat\gamma_j^{2}/\mathrm{se}_{Yj}^{2}\bigr)^{-1}\)。基于 \(\hat\beta_{\mathrm{fisher1}}\) 的推断并非最优，可实践里用得更广（Bowden et al., 2018）。\(\hat\beta_{\mathrm{fisher0}}\) 与 \(\hat\beta_{\mathrm{fisher1}}\) 都叫做固定效应估计量（fixed-effect estimator）。

盯那个并非最优、却更简单的 \(\hat\beta_{\mathrm{fisher1}}\)。在定义 25.2 下，可以说明依概率有

\[ \hat\beta_{\mathrm{fisher1}} \to \frac{\sum_{j=1}^{p}\Gamma_j\gamma_j/\mathrm{se}_{Yj}^{2}}{\sum_{j=1}^{p}\gamma_j^{2}/\mathrm{se}_{Yj}^{2}} = \beta + \frac{\sum_{j=1}^{p}\alpha_j\gamma_j/\mathrm{se}_{Yj}^{2}}{\sum_{j=1}^{p}\gamma_j^{2}/\mathrm{se}_{Yj}^{2}}. \]

若对所有 \(j\) 都有 \(\alpha_j=0\)，则 \(\hat\beta_{\mathrm{fisher1}}\) 相合。即便并非如此，只要 \(\alpha_j\) 与 \(\gamma_j\) 按 \(1/\mathrm{se}_{Yj}^{2}\) 加权的内积为零，它仍可能相合。这件事在遗传工具很多、并且排除限制的违背（由 \(\alpha_j\) 捕捉）是从均值为零的分布里独立抽出时，可以成立。

25.2.2　Egger 回归

从定义 25.1 出发。真参数满足

\[ \Gamma_j=\beta\gamma_j \qquad (j=1,\ldots,p); \]

换成估计，上面的恒等式只近似成立：

\[ \hat\Gamma_j\approx\beta\hat\gamma_j \qquad (j=1,\ldots,p). \]

这看起来像经典的最小二乘问题：把 \(\{\hat\Gamma_j\}_{j=1}^{p}\) 对 \(\{\hat\gamma_j\}_{j=1}^{p}\) 做回归。我们可以把 \(\hat\Gamma_j\) 对 \(\hat\gamma_j\) 做 WLS，带或不带截距，权重可以是 \(w_j\)，用来估 \(\beta\)。下面的结果来自附录 B.5 复习过的 WLS 代数性质。

不带截距时，\(\hat\gamma_j\) 的系数是

\[ \hat\beta_{\mathrm{egger1}} = \frac{\sum_{j=1}^{p}\hat\gamma_j\hat\Gamma_j w_j}{\sum_{j=1}^{p}\hat\gamma_j^{2} w_j}, \]

若 \(w_j=1/\mathrm{se}_{Yj}^{2}\)，它就化成 \(\hat\beta_{\mathrm{fisher1}}\)。这次 WLS 叫做 Egger 回归（Egger regression）。它比第 25.2.1 节的固定效应估计量更一般。带截距时，\(\hat\gamma_j\) 的系数是

\[ \hat\beta_{\mathrm{egger0}} = \frac{\sum_{j=1}^{p}(\hat\gamma_j-\hat\gamma_w)(\hat\Gamma_j-\hat\Gamma_w)w_j}{\sum_{j=1}^{p}(\hat\gamma_j-\hat\gamma_w)^{2}w_j}, \]

其中 \(\hat\gamma_w=\sum_{j=1}^{p}\hat\gamma_j w_j/\sum_{j=1}^{p}w_j\)、\(\hat\Gamma_w=\sum_{j=1}^{p}\hat\Gamma_j w_j/\sum_{j=1}^{p}w_j\) 分别是 \(\hat\gamma_j\) 与 \(\hat\Gamma_j\) 的加权平均。

即便不假定定义 25.2 下所有 \(\alpha_j\) 都是零，依概率仍有

\[ \hat\beta_{\mathrm{egger0}} \to \frac{\sum_{j=1}^{p}(\gamma_j-\gamma_w)(\Gamma_j-\Gamma_w)w_j}{\sum_{j=1}^{p}(\gamma_j-\gamma_w)^{2}w_j} = \beta + \frac{\sum_{j=1}^{p}(\gamma_j-\gamma_w)(\alpha_j-\alpha_w)w_j}{\sum_{j=1}^{p}(\gamma_j-\gamma_w)^{2}w_j}, \]

其中 \(\gamma_w\)、\(\Gamma_w\) 与 \(\alpha_w\) 是真参数相应的加权平均。因此，只要 \(\alpha_j\) 对 \(\gamma_j\) 的 WLS 系数为零，\(\hat\beta_{\mathrm{egger0}}\) 就对 \(\beta\) 相合。这比「所有 \(j\) 都有 \(\alpha_j=0\)」更弱。若 \(\gamma_j\) 与 \(\alpha_j\) 是相互独立的随机变量的实现，这条更弱的假定就成立；它叫做工具强度独立于直接效应（Instrument Strength Independent of Direct Effect, InSIDE）假定（Bowden et al., 2015）。更有意思的是，Egger 回归的截距是

\[ \hat\alpha_{\mathrm{egger0}}=\hat\Gamma_w-\hat\beta_{\mathrm{egger0}}\hat\gamma_w, \]

在 InSIDE 假定下，它依概率收敛到

\[ \Gamma_w-\beta\gamma_w=\alpha_w. \]

于是截距估的是直接效应的加权平均。

25.3　一个例子

先用下面的函数做 Fisher 加权。

fisher.weight = function(est, se)
{
  n = sum(est/se^2)
  d = sum(1/se^2)
  res = c(n/d, sqrt(1/d))
  names(res) = c("est", "se")
  res
}

我用 mr.raps 包（Zhao et al., 2020）里的 bmi.sbp 数据，说明基于不同方差估计的 Fisher 加权。\(\hat\beta_{\mathrm{fisher0}}\) 与 \(\hat\beta_{\mathrm{fisher1}}\) 的结果相近。

> bmisbp = read.csv("mr_bmisbp.csv")
> bmisbp$iv = with(bmisbp, beta.outcome/beta.exposure)
> bmisbp$se.iv = with(bmisbp, se.outcome/beta.exposure)
> bmisbp$se.iv1 = with(bmisbp,
+   sqrt(se.outcome^2 + iv^2*se.exposure^2)/beta.exposure)
> fisher.weight(bmisbp$iv, bmisbp$se.iv)
       est         se
0.31727680 0.05388827
> fisher.weight(bmisbp$iv, bmisbp$se.iv1)
       est         se
0.31576007 0.05893783

带或不带截距的 Egger 回归，点估计也很相近。可标准误与 Fisher 加权差得相当远。

> mr.egger = lm(beta.outcome ~ 0 + beta.exposure,
+               data = bmisbp,
+               weights = 1/se.outcome^2)
> summary(mr.egger)$coef
Estimate Std. Error \(t\) value \(\operatorname{Pr}(>\lvert t\rvert)\)
beta.exposure \(0.317\) \(0.111\) \(2.87\) \(0.0047\) 
> mr.egger.w = lm(beta.outcome ~ beta.exposure,
+                 data = bmisbp,
+                 weights = 1/se.outcome^2)
> summary(mr.egger.w)$coef
Estimate Std. Error \(t\) value \(\operatorname{Pr}(>\lvert t\rvert)\)
(Intercept) \(0.000\) \(0.002\) \(0.05\) \(0.957\)
beta.exposure \(0.317\) \(0.111\) \(2.86\) \(0.0048\)

原书图 25.2 画出原始数据，以及带截距的 Egger 回归拟合直线。

原书图 25.2：散点大小与方差倒数成比例，并叠上 Egger 回归直线。横轴是 \(\hat\gamma\)，纵轴是 \(\hat\Gamma\)。点多挤在原点附近，拟合直线略往上倾斜，截距几乎贴着零。

25.4　对基于 MR 的分析的几条批评

MR 是 IV 想法的一种应用。它依赖很强的假定。我从概念、生物学与技术三个角度给三组批评。

概念上，大多数基于 MR 的研究里，处理从潜在结果的角度看定义并不清楚。例如，处理常常被定义成胆固醇水平或 BMI。它们是复合变量，可以对应许多复杂、并不唯一的假想实验。对这些处理，SUTVA 常常不成立。回想第 2 章的讨论。

生物学上，IV 分析那些根本假定未必成立。孟德尔第二定律保证不同性状的遗传相互独立。可它并不保证候选 IV 与处理、结果之间的隐藏混杂独立。这些 IV 可能对混杂有直接效应。也可能有一些未测基因，同时影响 IV 与混杂。孟德尔第二定律也不保证排除限制。IV 完全可能沿着处理之外的其他因果路径走到结果。

技术上，MR 的统计假定相当强。线性 IV 模型本身就是一条强建模假定。\(\hat\gamma_j\) 之间、\(\hat\Gamma_j\) 之间的独立也很强。数据收集过程里的另一些麻烦，还会让 IV 假定更难解释。例如，处理与结果常常带测量误差；全基因组关联研究又常常基于病例对照设计，样本是条件于结果的（见附录 B.6.3）。

VanderWeele et al. (2014) 是一篇很好的综述，讨论 MR 在方法上的挑战。

25.5　习题

25.1　数据分析
分析 R 包 mr.raps 里的 bmi.bmi 数据。更多细节见该包以及 Zhao et al. (2020, 第 7.2 节)。

25.2　推荐阅读
Davey Smith and Ebrahim (2003) 综述过 MR 的潜力与局限。

第六部分　带处理后变量的因果机制

  1. 多效性指：一个基因影响两个或多个表面上不相关的表型性状。↩︎