Stata-双样本孟德尔随机化法 (MR)-mrrobust

发布时间：2021-11-13 阅读 3132

Stata连享会主页 || 视频 || 推文 || 知乎 || Bilibili 站

温馨提示： 定期清理浏览器缓存，可以获得最佳浏览体验。

New！ lianxh 命令发布了：
随时搜索推文、Stata 资源。安装：
. ssc install lianxh
详情参见帮助文件 (有惊喜)：
. help lianxh
连享会新命令：cnssc, ihelp, rdbalance, gitee, installpkg

✌ 课程详情： https://gitee.com/lianxh/Course

⛳ 课程主页： https://gitee.com/lianxh/Course

⛳ Stata 系列推文：

全部 | Stata入门 | Stata教程 | Stata资源 | Stata命令
计量专题 | 论文写作 | 数据分享 | 专题课程
结果输出 | Stata绘图 | 数据处理 | Stata程序
回归分析 | 面板数据 | 交乘项-调节 | IV-GMM
内生性-因果推断 | 倍分法DID | 断点回归RDD | PSM-Matching | 合成控制法
Probit-Logit | 时间序列 | 空间计量 | 分位数回归 | 生存分析 | SFA-DEA
文本分析-爬虫 | Python-R-Matlab | 机器学习
Markdown | 工具软件 | 其它

☝ PDF下载 - 推文合集

作者：彭甲超 (中国地质大学)
邮箱：pengjiachao@cug.edu.cn

编者按：本文摘译自下文，特此致谢！
Source：Spiller W, Davies N M, Palmer T M. Software application profile: mrrobust—a tool for performing two-sample summary Mendelian randomization analyses[J]. International journal of epidemiology, 2019, 48(3): 684. -PDF-

1. 简介

孟德尔随机化法 (Mendelian randomization，MR) 已经发展成为检验流行病学因果推断关系的一种有效方法。MR 方法通过使用遗传变异作为工具变量 (IV) 来推导暴露与结果变量之间的因果关系。要使得 MR 方法实现有效无偏估计，需要满足三条假设：

工具变量必须与暴露变量相关；
工具变量和结果变量不能有混杂因素；
工具变量不能影响结果变量，工具变量只能通过暴露变量与结果变量相关。

本文介绍了 Stata 的 mrrobust 软件包，可以实现双样本 MR analyses 分析过程。在 Stata 中，可以使用 IVREG2 实现两个样本等模块的单个级别 IV 分析。mrust 软件包解决了这一限制，提供了一套流行的双样本 MR 方法和灵敏度分析。接下来，本文将简要介绍几种 MR 估计方法。

考虑 MR 研究中关于遗传变异基因样本 $G_{1}, . . ., G_{J}$ ，存在连续暴露变量 $X$ 和结果变量 $Y$ 。所有混淆因素都包含在变量 $U$ 中。假设所有的遗传变体都是有效的 IV，并进一步假设变量之间的所有关系都是线性的，且不存在异质性 (Bowden 等，2016)：

X ∣ G_{j} = γ_{0} + γ_{j} G_{j} + ϵ_{X j}

Y ∣ G_{j} = Γ_{0} + Γ_{j} G_{j} + ϵ_{Y j}

假设所有的变异都与暴露变量有关，则对于所有的 $j$ 有 $γ_{j} \neq 0$ 。参考 Bowden 等 (2016)，设定 $Γ_{j} = β γ_{j}$ 。误差项 $ϵ_{X j}$ 和 $ϵ_{Y j}$ 为正态分布，存在相关性，包含混杂项 $U$ 和 $G_{j}$ 。

逆方程加权估计 (IVW) 估计：暴露变量 $X$ 对结果变量 $Y$ 的因果影响，可以使用第 $j$ 个变量作为遗传结果关联和基因暴露关联估计的比率来估计：

{\hat{β}}_{j} = \frac{{\hat{Γ}}_{j}}{{\hat{γ}}_{j}}

如果遗传方差满足 IV 假设，则 $Γ_{j} = β γ_{j}$ ，并且比率估计是渐近一致的。此外，如果遗传变异是不相关的，则可以使用 meta-analysis 中的公式将来自每个遗传变异的比率估计值组合成总体估计值：

{\hat{β}}_{I V W} = \frac{\sum_{j} {\hat{γ}}_{j}^{2} σ_{Y j}^{- 2} {\hat{β}}_{j}}{\sum_{j} {\hat{γ}}_{j}^{2} σ_{Y j}^{- 2}}

其中 $σ_{Y_{j}}$ 是变量 $j$ 的基因-结果关联估计的标准误差，这被称为逆方差加权 (IVW) 估计量。假设遗传变异是不相关的，IVW 估计与通常用于个体水平数据的两阶段最小二乘估计渐近相等。如果所有的遗传变异都满足 IV 假设，则 IVW 估计是因果效应的一致估计 (即，随着样本量的增加而收敛到真值)，因为它是个体比率估计的加权平均值。

Mr-Egger 回归估计：Bowden 等 (2015) 提出了一种对加总数据进行 MR 的稳健方法，即 Mr-Egger 回归。这种方法源于 Meta 分析文献中用于评估小型研究偏差 (发表偏倚) 的一种方法 (Egger 等，1997)，该方法会对基因暴露系数 ${\hat{γ}}_{j}$ 上的基因结果系数 ${\hat{Γ}}_{j}$ 进行加权线性回归：

{\hat{Γ}}_{j} = β_{0 E} + β_{E} {\hat{γ}}_{j}

其中所有的 ${\hat{γ}}_{j}$ 关联被设定为正，并且回归中的权重是基因-结果关联的逆方差 $(σ_{Y j}^{- 2})$ 。如果回归模型中没有截距项，则 MR-Egger 斜率估计 ${\hat{β}}_{E}$ 将等于 IVW 估计。

中位数估计：当所有的遗传变异都是有效的 IV 估计时，IVW 估计是一种有效的分析方法。但是，即使只有一个基因变异是无效的，它也是有偏见的。因此，IVW 预估可以说是有 0% 的细分水平。中位数估计比率提供了对因果效应的一致估计。

具体地说，设 ${\hat{β}}_{j}$ 表示第 $j$ 个有序比估计 (从小到大排列)。如果遗传变异总数为奇数 $(J = 2 k + 1)$ ，则简单的中位数估计为中比估计 ${\hat{β}}_{k + 1}$ ；如果为偶数 $(J = 2 k)$ ，则将中位数内插在两个中值估计 $\frac{1}{2} ({\hat{β}}_{k} + {\hat{β}}_{k + 1})$ 之间。具体内容及解释可参见 Bowden 等 (2016) 的图 2，当然中位数估计可以继续升级为加权中位数估计 (Weighted Median Estimator) 和惩罚加权中位数估计 (Penalized Weighted Median Estimator)。

2. 命令介绍

mrrobust 程序包含一组命令，用于使用基因型-表型和基因型-结果关联的汇总数据执行两样本 MR 分析。具体为基础命令 mr、安装依赖项 mrdeps、MR-Egger 和逆方差加权 (IVW) 估计 mregger、MR-Egger 模型的模拟外推算法 mreggersimex、散点图显示特定估计值拟合线和置信区间 mreggerplot、加总数据中值估计 mrmedian、个体数据中值估计 mrmedianobs、加总数据模态估计 mrmodal、模态估计密度图 mrmodalplot、Ratio (Wald) 估计 mratio、加总数据生成比率 (Wald) 估计 mrivests、森林图命令 mrdorest、IV 估计漏斗图 mrfunnel、多变量逆方差加权估计 mrmvivw、多变量 MR-Egger 回归估计 mrmvegger、Leave one out 分析估计 mrleaveoneout。

*命令安装
cnssc install lxhget, replace
lxhget mrrobust.pkg, install replace

*命令语法
mr subcommand ... [aweight] [if] [in] [, options]
mregger varname_gd varname_gp [aweight] [if] [in] [, options]
mreggersimex varname_gd varname_gp [aweight] [if] [in] [, options]
mreggerplot varname_gd varname_gdse varname_gp varname_gpse [if] [in] [, options]
mrmedian varname_gd varname_gdse varname_gp varname_gpse [if] [in] [, options]
mrmedianobs depvar [varlist1] (varname_endog = varlist_ivs) [if] [in] [, options]
mrmodal varname_gd varname_gdse varname_gp varname_gpse [if] [in] [, options]
mrmodalplot varname_gd varname_gdse varname_gp varname_gpse [if] [in] [, options]
mrratio #gd #gdse #gp [#gpse #cov] [, options]
mrivests varname_gd varname_gdse varname_gp [varname_gpse varname_cov] [if] [in] [, options]
mrforest varname_gd varname_gdse varname_gp varname_gpse [varname_cov] [if] [in] [, options]
mrfunnel varname_gd varname_gdse varname_gp varname_gpse [if] [in] [, options]
mrmvivw/mvivw/mvmr varname_gd varname_gp1 [varname_gp2 ...] [aweight] [if] [in] [, options]
mrmvegger varname_gd varname_gp1 [varname_gp2 ...] [aweight] [if] [in] [, options]
mrleaveoneout varname_gd varname_gp [if] [in], gyse(varname) [options]

varname_gd 和 varname_gdse 表示结果变量；
varname_gp 和 varname_gpse 表示暴露变量。

3. Stata 实例

本节通过 Bowden 等 (2016) 示例数据展示主要命令。导入示例数据，获取与结果变量 chdbeta 和暴露变量 ldlcbeta 相关的样本，同时生成不同的子样本。

. cnssc install lxhuse, replace
. lxhuse dodata.dta, clear

. *Select observations (p-value with exposure < 10^-8)
. gen byte sel1 = (ldlcp2 < 1e-8)

3.1 IVW 估计应用

固定效应 IVW 假设每个工具变量都能得到相同的估计，也就是没有水平多效性的 IV。随机效应逆方差加权放宽了这一假设，允许每个 IV 变量有不同的效应估计，在平衡水平多效性的情况下，仍可以提供无偏估计。

固定效应逆方差加权模型和随机效应的估计值是相同的，但是由于考虑到工具变量之间的异质性，随机效应模型的方差相对较大。以下仅展示固定效率 IVW 的估计结果，并绘制森林图 (部分)：

. *IVW (with fixed effect standard errors)
. mr egger chdbeta ldlcbeta [aw=1/(chdse^2)] if sel1==1, ivw fe

                                                      Number of genotypes = 73
                                      Residual standard error constrained at 1
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
chdbeta      |
    ldlcbeta |      0.482      0.038    12.60   0.000        0.407       0.556
------------------------------------------------------------------------------

. ssc install metan, replace
. mrforest chdbeta chdse ldlcbeta ldlcse if sel1==1, ivid(rsid) ///
>          xlabel(-5,-4,-3,-2,-1,0,1,2,3,4,5)
. qui gr export mrforest.svg, width(600) replace

3.2 Egger 回归应用

与普通的 MR 方法相比，Egger 回归方法放宽了工具变量-无水平多向性假设。MR-Egger 相比于 IVW 方法，Egger 回归允许了非零的回归方程截距：

β_{Y_{j}} = α_{j} + θ β_{X_{j}}

使用上述公式可以接受所有工具变量的水平多效性是不平衡或方向性的。在这种情形下，MR-Egger 方法仍可以得到无偏因果效应估计，但必须满足水平多效性效应与工具变量对暴露变量 $X$ 的效应不相关。下文展示了多种 Egger 回归方式：

固定效应 IVW 的 Egger 回归

. mregger chdbeta ldlcbeta [aw=1/(chdse^2)] if sel1==1, ivw fe

                                                      Number of genotypes = 73
                                      Residual standard error constrained at 1
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
chdbeta      |
    ldlcbeta |      0.482      0.038    12.60   0.000        0.407       0.556
------------------------------------------------------------------------------

使用无约束残差方差「乘法随机效应」的 Egger 回归

. mregger chdbeta ldlcbeta [aw=1/(chdse^2)] if sel1==1

                                                      Number of genotypes = 73
                                              Residual standard error =  1.548
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
chdbeta      |
       slope |      0.617      0.103     5.97   0.000        0.415       0.820
       _cons |     -0.009      0.005    -1.60   0.110       -0.020       0.002
------------------------------------------------------------------------------

报告 $I_{G X}^{2}$ 统计量和异质性 Q 检验的 Egger 回归。异质性 Q 检验公式为：

Q = \sum_{j} Q_{j} = \sum_{j} w_{j}^{'} {({\hat{β}}_{j} - \hat{β})}^{2}

. ssc install heterogi, replace
. mregger chdbeta ldlcbeta [aw=1/(chdse^2)] if sel1==1, ///
> >       gxse(ldlcse) heterogi

                                           Q_GX statistic (weighted) = 3454.26
                                         I^2_GX statistic (weighted) =  97.92%
                                                      Number of genotypes = 73
                                              Residual standard error =  1.548
                Ruecker's Q for heterogeneity; chi2(71) = 170.11 (p =  0.0000)
                             I-squared statistic = 58.3% (95% CI 45.8%, 67.8%)
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
chdbeta      |
       slope |      0.617      0.103     5.97   0.000        0.415       0.820
       _cons |     -0.009      0.005    -1.60   0.110       -0.020       0.002
------------------------------------------------------------------------------

使用 $t$ 分布检验的 Egger 回归

. mregger chdbeta ldlcbeta [aw=1/(chdse^2)] if sel1==1, tdist

                                                      Number of genotypes = 73
                                              Residual standard error =  1.548
------------------------------------------------------------------------------
             | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
chdbeta      |
       slope |      0.617      0.103     5.97   0.000        0.411       0.824
       _cons |     -0.009      0.005    -1.60   0.114       -0.020       0.002
------------------------------------------------------------------------------

使用径向公式估计的 Egger 回归

. mregger chdbeta ldlcbeta [aw=1/(chdse^2)] if sel1==1, radial

                                                      Number of genotypes = 73
                                              Residual standard error =  1.547
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
radialGD     |
    radialGP |      0.643      0.116     5.55   0.000        0.416       0.870
       _cons |     -0.574      0.355    -1.62   0.106       -1.269       0.121
------------------------------------------------------------------------------

使用径向公式和报告异质性 (Rucker's) Q 检验的 Egger 回归

. mregger chdbeta ldlcbeta [aw=1/(chdse^2)] if sel1==1, ///
>         radial heterogi

                                                      Number of genotypes = 73
                                              Residual standard error =  1.547
                Ruecker's Q for heterogeneity; chi2(71) = 169.98 (p =  0.0000)
                             I-squared statistic = 58.2% (95% CI 45.8%, 67.8%)
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
radialGD     |
    radialGP |      0.643      0.116     5.55   0.000        0.416       0.870
       _cons |     -0.574      0.355    -1.62   0.106       -1.269       0.121
------------------------------------------------------------------------------

Egger 回归的自举法应用及可视化

. net install grc1leg, from(http://www.stata.com/users/vwiggins)
. mreggersimex chdbeta ldlcbeta [aw=1/chdse^2] if sel1==1, ///
>              gxse(ldlcse) seed(12345) noboot

                                                      Number of genotypes = 73
                                                    Bootstrap replications = 0
                                                  Simulation replications = 50
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       slope |      0.626          .        .       .            .           .
       _cons |     -0.009          .        .       .            .           .
------------------------------------------------------------------------------
. qui gr export mreggersimex-plot.svg, width(600) replace

MR-Egger 回归散点图示

. ssc install addplot, replace
. mreggerplot chdbeta chdse ldlcbeta ldlcse if sel1==1
. qui gr export mreggerplot.svg, width(600) replace

3.3 中位数估计应用

中位数估计方法的优势在于，只要存在一半的工具变量为有效 IV，即无水平多效性，与混杂无关联，且与暴露变量 $X$ 强相关，则因果效应估计是无偏的。加权中位数估计则是通过将每个工具变量的贡献用工具变量与结果变量的相关性的逆方差加权，从而使得较强的工具变量对估计值的贡献更大。

加权中值估计应用。mrmedian 命令中的 chdbeta 和 chdse 为结果变量，ldlcbeta 和 ldlcse 为暴露变量。

. mrmedian chdbeta chdse ldlcbeta ldlcse if sel1==1, weighted

                                                      Number of genotypes = 73
                                                           Replications = 1000
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
        beta |      0.458      0.062     7.33   0.000        0.336       0.581
------------------------------------------------------------------------------

mrmodal 估计及可视化应用。mrmodal 实现了加总数据的 Hartwig 等 (2017) 零模态估计，即结果变量和暴露变量关联估计及其个体基因型的标准误差。mrmodal 的 options 主要包括：level(#) 表示置信区间，默认95%；nome 表示 NOME assumption；nosave 表示不保存 IV 估计的密度值；phi(#) 表示带宽值；reps(#) 表示获得标准误差的重复次数；seed(#) 表示用于获得标准误差的随机数生成器的种子；weighted 表示加权IV 估计。主要结果如下所示：

. ssc install kdens, replace
. mrmodal chdbeta chdse ldlcbeta ldlcse if sel1==1

                                                      Number of genotypes = 73
                                                           Replications = 1000
                                                                       Phi = 1
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
        beta |      0.492      0.130     3.79   0.000        0.237       0.746
------------------------------------------------------------------------------

. mrmodalplot chdbeta chdse ldlcbeta ldlcse if sel1==1, seed(12345)

                                                      Number of genotypes = 73
                                                           Replications = 1000
                                                                     Phi = .25
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
        beta |      0.420      0.227     1.85   0.064       -0.025       0.864
------------------------------------------------------------------------------
                                                      Number of genotypes = 73
                                                           Replications = 1000
                                                                      Phi = .5
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
        beta |      0.422      0.198     2.13   0.033        0.034       0.810
------------------------------------------------------------------------------
                                                      Number of genotypes = 73
                                                           Replications = 1000
                                                                       Phi = 1
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
        beta |      0.492      0.136     3.63   0.000        0.226       0.758
------------------------------------------------------------------------------

mrmodal加权估计应用

. mrmodal chdbeta chdse ldlcbeta ldlcse if sel1==1, weighted

                                                      Number of genotypes = 73
                                                           Replications = 1000
                                                                       Phi = 1
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
        beta |      0.479      0.067     7.13   0.000        0.347       0.611
------------------------------------------------------------------------------

mrmodal NOME 假设估计应用

. mrmodal chdbeta chdse ldlcbeta ldlcse if sel1==1, nome

                                                      Number of genotypes = 73
                                                           Replications = 1000
                                                                       Phi = 1
------------------------------------------------------------------------------
             | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
        beta |      0.492      0.129     3.81   0.000        0.239       0.745
------------------------------------------------------------------------------

3.4 漏斗图示例

. mrfunnel chdbeta chdse ldlcbeta ldlcse if sel1==1, xlrange(0 10)
. qui gr export mrfunnel.svg, width(600) replace

3.5 Leave-one-out 分析可视化应用

为了评估 MR 估计是否由一个具有特别大的水平多效性效应的单个工具变量 (异常值) 影响或偏倚，我们可以通过轮流删除一个工具变量来重新估计因果效应，然后判断删除一个工具变量是否对因果效应的估计影响很大，从而识别强影响工具变量。Leave-one-out 分析主要汇报工具变量的森林图，提供两种估计方法，如下图所示：

. gen byte sel2 = (ldlcp2 < 1e-25)
. mrleaveoneout chdbeta ldlcbeta if sel2==1, gyse(chdse) genotype(rsid) noprint
. qui gr export mrleaveoneout-plot-01.svg, width(600) replace

. mrleaveoneout chdbeta ldlcbeta hdlcbeta tgbeta if sel2==1, ///
>               method(mvmr) gyse(chdse) genotype(rsid) noprint
. qui gr export mrleaveoneout-plot-02.svg, width(600) replace

4. 结语

通常，我们很容易混淆回归中的相关性和因果性，那么我们在实际操作中该如何解释本身具有强相关性事件之间的因果关系呢？MR 方法在解决此类问题上发挥重要作用，并被广泛应用于科研领域。MR的原理并不复杂，但是在实际操作中，还是需要注意几点：

一是找到合适的基因变异作为工具变量，该工具变量必须与 “暴露” 变量是明确相关且不能与混杂因素相关，而 “暴露” 变量与结果变量之间的关联一定会受到潜在因素的影响；
二是工具变量只能通过一条途径影响结果变量，即工具变量通过 “暴露” 变量影响结果变量而不是其他；
三是 “暴露” 变量与结果变量之间的关联无法直接观察获得，因为无法直接测量 “暴露” 变量，但工具变量是可测量的，并且工具变量与 “暴露” 变量直接的关联是已知的或者可测量的，并独立于其他因素而存在。

5. 参考资料

Bowden J, Davey Smith G, Burgess S, 2015, Mendelian randomization with invalid instruments: effect estimation and bias detection through Egger regression, International Journal of Epidemiology, 44, 2: 512-525. -PDF-
Bowden J, Davey Smith G, Haycock PC, Burgess S, 2016, Consistent estimation in Mendelian randomization with some invalid instruments using a weighted median estimator, Genetic Epidemiology, 40, 4: 304-314. -PDF-
Hartwig F P, Davey Smith G, Bowden J, 2017, Robust inference in two-sample Mendelian randomisation via the zero modal pleiotropy assumption, International Journal of Epidemiology, 46, 6: 1985-1998. -PDF-
Hemani G et al., 2018, The MR-Base platform supports systematic causal inference across the human phenome, eLife, 7:e34408. -PDF-
Sanderson E, Davey Smith G, Windmeijer F, Bowden J, 2019, An examination of multivariable Mendelian randomization in the single-sample and two-sample summary data settings, International Journal of Epidemiology, 48, 3: 713-727. -PDF-
Sanderson E, Spiller W, Bowden J, 2020, Testing and correcting for weak and pleiotropic instruments in two-sample multivariable Mendelian randomisation, bioRxiv, 2020.04.02.021980. -PDF-

6. 相关推文

Note：产生如下推文列表的 Stata 命令为：
lianxh IV, m
安装最新版 lianxh 命令：
ssc install lianxh, replace

专题	嘉宾	直播/回看视频
⭐ 最新专题		文本分析、机器学习、效率专题、生存分析等
研究设计	连玉君	我的特斯拉-实证研究设计，-幻灯片-
面板模型	连玉君	动态面板模型，-幻灯片-
面板模型	连玉君	直击面板数据模型 [免费公开课，2小时]

Stata-双样本孟德尔随机化法 (MR)-mrrobust

1. 简介

2. 命令介绍

3. Stata 实例

3.1 IVW 估计应用

3.2 Egger 回归应用

3.3 中位数估计应用

3.4 漏斗图示例

3.5 Leave-one-out 分析可视化应用

4. 结语

5. 参考资料

6. 相关推文

相关课程

免费公开课

最新课程-直播课

⛳ 课程主页

⛳ 课程主页

关于我们

热门资讯