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

发布时间:2021-11-13 阅读 172

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 系列推文:

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 研究中关于遗传变异基因样本 G1,...,GJ,存在连续暴露变量 X 和结果变量 Y。所有混淆因素都包含在变量 U 中。假设所有的遗传变体都是有效的 IV,并进一步假设变量之间的所有关系都是线性的,且不存在异质性 (Bowden 等,2016):

假设所有的变异都与暴露变量有关,则对于所有的 j 有 γj0。参考 Bowden 等 (2016),设定 Γj=βγj。误差项 ϵXj 和 ϵYj 为正态分布,存在相关性,包含混杂项 U 和 Gj

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

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

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

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

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

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

具体地说,设 β^j 表示第 j 个有序比估计 (从小到大排列)。如果遗传变异总数为奇数 (J=2k+1),则简单的中位数估计为中比估计 β^k+1;如果为偶数 (J=2k),则将中位数内插在两个中值估计 12(β^k+β^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_gdvarname_gdse 表示结果变量;
  • varname_gpvarname_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 回归允许了非零的回归方程截距:

使用上述公式可以接受所有工具变量的水平多效性是不平衡或方向性的。在这种情形下,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
------------------------------------------------------------------------------
  • 报告 IGX2 统计量和异质性 Q 检验的 Egger 回归。异质性 Q 检验公式为:
. 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 命令中的 chdbetachdse 为结果变量,ldlcbetaldlcse 为暴露变量。
. 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) 零模态估计,即结果变量和暴露变量关联估计及其个体基因型的标准误差。mrmodaloptions 主要包括: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小时]
  • Note: 部分课程的资料,PPT 等可以前往 连享会-直播课 主页查看,下载。

课程主页

课程主页

关于我们

  • Stata连享会 由中山大学连玉君老师团队创办,定期分享实证分析经验。
  • 连享会-主页知乎专栏,700+ 推文,实证分析不再抓狂。直播间 有很多视频课程,可以随时观看。
  • 公众号关键词搜索/回复 功能已经上线。大家可以在公众号左下角点击键盘图标,输入简要关键词,以便快速呈现历史推文,获取工具软件和数据下载。常见关键词:课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法

连享会小程序:扫一扫,看推文,看视频……

扫码加入连享会微信群,提问交流更方便

✏ 连享会-常见问题解答:
https://gitee.com/lianxh/Course/wikis

New! lianxhsongbl 命令发布了:
随时搜索连享会推文、Stata 资源,安装命令如下:
. ssc install lianxh
使用详情参见帮助文件 (有惊喜):
. help lianxh