Stata连享会 主页 || 视频 || 推文 || 知乎 || Bilibili 站
温馨提示: 定期 清理浏览器缓存,可以获得最佳浏览体验。
New!
lianxh
命令发布了:
随时搜索推文、Stata 资源。安装:
. ssc install lianxh
详情参见帮助文件 (有惊喜):
. help lianxh
连享会新命令:cnssc
,ihelp
,rdbalance
,gitee
,installpkg
⛳ Stata 系列推文:
作者:彭甲超 (中国地质大学)
邮箱: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-
目录
孟德尔随机化法 (Mendelian randomization,MR) 已经发展成为检验流行病学因果推断关系的一种有效方法。MR 方法通过使用遗传变异作为工具变量 (IV) 来推导暴露与结果变量之间的因果关系。要使得 MR 方法实现有效无偏估计,需要满足三条假设:
本文介绍了 Stata 的 mrrobust
软件包,可以实现双样本 MR analyses 分析过程。在 Stata 中,可以使用 IVREG2
实现两个样本等模块的单个级别 IV 分析。mrust
软件包解决了这一限制,提供了一套流行的双样本 MR 方法和灵敏度分析。接下来,本文将简要介绍几种 MR 估计方法。
考虑 MR 研究中关于遗传变异基因样本
假设所有的变异都与暴露变量有关,则对于所有的
逆方程加权估计 (IVW) 估计:暴露变量
如果遗传方差满足 IV 假设,则
其中
Mr-Egger 回归估计:Bowden 等 (2015) 提出了一种对加总数据进行 MR 的稳健方法,即 Mr-Egger 回归。这种方法源于 Meta 分析文献中用于评估小型研究偏差 (发表偏倚) 的一种方法 (Egger 等,1997),该方法会对基因暴露系数
其中所有的
中位数估计:当所有的遗传变异都是有效的 IV 估计时,IVW 估计是一种有效的分析方法。但是,即使只有一个基因变异是无效的,它也是有偏见的。因此,IVW 预估可以说是有 0% 的细分水平。中位数估计比率提供了对因果效应的一致估计。
具体地说,设
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
表示暴露变量。
本节通过 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)
固定效应 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
与普通的 MR 方法相比,Egger 回归方法放宽了工具变量-无水平多向性假设。MR-Egger 相比于 IVW 方法,Egger 回归允许了非零的回归方程截距:
使用上述公式可以接受所有工具变量的水平多效性是不平衡或方向性的。在这种情形下,MR-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
------------------------------------------------------------------------------
. 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
------------------------------------------------------------------------------
. 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
------------------------------------------------------------------------------
. 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
------------------------------------------------------------------------------
. 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
------------------------------------------------------------------------------
. 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
------------------------------------------------------------------------------
. 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
. ssc install addplot, replace
. mreggerplot chdbeta chdse ldlcbeta ldlcse if sel1==1
. qui gr export mreggerplot.svg, width(600) replace
中位数估计方法的优势在于,只要存在一半的工具变量为有效 IV,即无水平多效性,与混杂无关联,且与暴露变量
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
------------------------------------------------------------------------------
. mrfunnel chdbeta chdse ldlcbeta ldlcse if sel1==1, xlrange(0 10)
. qui gr export mrfunnel.svg, width(600) replace
为了评估 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
通常,我们很容易混淆回归中的相关性和因果性,那么我们在实际操作中该如何解释本身具有强相关性事件之间的因果关系呢?MR 方法在解决此类问题上发挥重要作用,并被广泛应用于科研领域。MR的原理并不复杂,但是在实际操作中,还是需要注意几点:
Note:产生如下推文列表的 Stata 命令为:
lianxh IV, m
安装最新版lianxh
命令:
ssc install lianxh, replace
免费公开课
最新课程-直播课
专题 | 嘉宾 | 直播/回看视频 |
---|---|---|
⭐ 最新专题 | 文本分析、机器学习、效率专题、生存分析等 | |
研究设计 | 连玉君 | 我的特斯拉-实证研究设计,-幻灯片- |
面板模型 | 连玉君 | 动态面板模型,-幻灯片- |
面板模型 | 连玉君 | 直击面板数据模型 [免费公开课,2小时] |
⛳ 课程主页
⛳ 课程主页
关于我们
课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法
等
连享会小程序:扫一扫,看推文,看视频……
扫码加入连享会微信群,提问交流更方便
✏ 连享会-常见问题解答:
✨ https://gitee.com/lianxh/Course/wikis
New!
lianxh
和songbl
命令发布了:
随时搜索连享会推文、Stata 资源,安装命令如下:
. ssc install lianxh
使用详情参见帮助文件 (有惊喜):
. help lianxh