👇 连享会 · 推文导航 | www.lianxh.cn
🍎 Stata：Stata基础 | Stata绘图 | Stata程序 | Stata新命令 📘 论文：数据处理 | 结果输出 | 论文写作 | 数据分享 💹 计量：回归分析 | 交乘项-调节 | IV-GMM | 时间序列 | 面板数据 | 空间计量 | Probit-Logit | 分位数回归 ⛳ 专题：SFA-DEA | 生存分析 | 爬虫 | 机器学习 | 文本分析 🔃 因果：DID | RDD | 因果推断 | 合成控制法 | PSM-Matching 🔨 工具：工具软件 | Markdown | Python-R-Stata 🎧 课程：最新专题 | 计量专题 | 关于连享会
🍓 课程推荐：2024 空间计量专题
主讲老师：范巧 (兰州大学)
课程时间：2024 年 10 月 2-4 日 (三天)
课程咨询：王老师 18903405450（微信）
课程特色 · 2024空间计量：
👉 一、从“零基础”到“高水平”的课程设计
兼顾基础知识、主流模型与前沿模型 既考虑软件安装、程序编写以及空间权重矩阵设计等 基础知识 讲授，更强调时空面板地理加权回归模型、贝叶斯空间计量模型、矩阵指数模型、空间计量交互模型与空间面板似不相关回归模型等 前沿模型 的传授。
👉 二、“保姆级”的空间计量代码
编写与校准所有模型的MATLAB代码，简化实操环节 模型的估计与检验等 仅按照提供的Excel数据版式 搜集与整理原始数据，即可一次性出结果并作图。
👉 三、“最多上新” 的内容体系
新增 矩阵指数模型、短面板空间似不相关模型、空间计量交互模型、贝叶斯空间计量模型等 新增 前沿应用案例，包括空间计量与索洛余值法、随机前沿分析与数据包络分析等的互嵌研究，阐释基于空间计量的产业空间结构优化评价方法。 新增 Dagum空间基尼系数、核密度估计、空间马尔科夫链与空间收敛性等内容，阐释现实研究中对空间收敛性的应用“谬误”。
作者：彭甲超 (中国地质大学)
邮箱：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. 简介
2. 命令介绍
3. Stata 实例
3.1 IVW 估计应用
3.2 Egger 回归应用
3.3 中位数估计应用
3.4 漏斗图示例
3.5 Leave-one-out 分析可视化应用
4. 结语
5. 参考资料
6. 相关推文
温馨提示： 文中链接在微信中无法生效。请点击底部「阅读原文」。或直接长按/扫描如下二维码，直达原文：
1. 简介
孟德尔随机化法 (Mendelian randomization，MR) 已经发展成为检验流行病学因果推断关系的一种有效方法。MR 方法通过使用遗传变异作为工具变量 (IV) 来推导暴露与结果变量之间的因果关系。要使得 MR 方法实现有效无偏估计，需要满足三条假设：
工具变量必须与暴露变量相关； 工具变量和结果变量不能有混杂因素； 工具变量不能影响结果变量，工具变量只能通过暴露变量与结果变量相关。
本文介绍了 Stata 的 mrrobust 软件包，可以实现双样本 MR analyses 分析过程。在 Stata 中，可以使用 IVREG2 实现两个样本等模块的单个级别 IV 分析。mrust 软件包解决了这一限制，提供了一套流行的双样本 MR 方法和灵敏度分析。接下来，本文将简要介绍几种 MR 估计方法。
考虑 MR 研究中关于遗传变异基因样本 ，存在连续暴露变量 和结果变量 。所有混淆因素都包含在变量 中。假设所有的遗传变体都是有效的 IV，并进一步假设变量之间的所有关系都是线性的，且不存在异质性 (Bowden 等，2016)：
假设所有的变异都与暴露变量有关，则对于所有的 有 。参考 Bowden 等 (2016)，设定 。误差项 和 为正态分布，存在相关性，包含混杂项 和 。
逆方程加权估计 (IVW) 估计：暴露变量 对结果变量 的因果影响，可以使用第 个变量作为遗传结果关联和基因暴露关联估计的比率来估计：
如果遗传方差满足 IV 假设，则 ，并且比率估计是渐近一致的。此外，如果遗传变异是不相关的，则可以使用 meta-analysis 中的公式将来自每个遗传变异的比率估计值组合成总体估计值：
其中 是变量 的基因-结果关联估计的标准误差，这被称为逆方差加权 (IVW) 估计量。假设遗传变异是不相关的，IVW 估计与通常用于个体水平数据的两阶段最小二乘估计渐近相等。如果所有的遗传变异都满足 IV 假设，则 IVW 估计是因果效应的一致估计 (即，随着样本量的增加而收敛到真值)，因为它是个体比率估计的加权平均值。
Mr-Egger 回归估计：Bowden 等 (2015) 提出了一种对加总数据进行 MR 的稳健方法，即 Mr-Egger 回归。这种方法源于 Meta 分析文献中用于评估小型研究偏差 (发表偏倚) 的一种方法 (Egger 等，1997)，该方法会对基因暴露系数 上的基因结果系数 进行加权线性回归：
其中所有的 关联被设定为正，并且回归中的权重是基因-结果关联的逆方差 。如果回归模型中没有截距项，则 MR-Egger 斜率估计 将等于 IVW 估计。
中位数估计：当所有的遗传变异都是有效的 IV 估计时，IVW 估计是一种有效的分析方法。但是，即使只有一个基因变异是无效的，它也是有偏见的。因此，IVW 预估可以说是有 0% 的细分水平。中位数估计比率提供了对因果效应的一致估计。
具体地说，设 表示第 个有序比估计 (从小到大排列)。如果遗传变异总数为奇数 ，则简单的中位数估计为中比估计 ；如果为偶数 ，则将中位数内插在两个中值估计 之间。具体内容及解释可参见 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 回归允许了非零的回归方程截距：
使用上述公式可以接受所有工具变量的水平多效性是不平衡或方向性的。在这种情形下，MR-Egger 方法仍可以得到无偏因果效应估计，但必须满足水平多效性效应与工具变量对暴露变量的效应不相关。下文展示了多种 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
------------------------------------------------------------------------------
报告 统计量和异质性 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
------------------------------------------------------------------------------
使用 分布检验的 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，即无水平多效性，与混杂无关联，且与暴露变量 强相关，则因果效应估计是无偏的。加权中位数估计则是通过将每个工具变量的贡献用工具变量与结果变量的相关性的逆方差加权，从而使得较强的工具变量对估计值的贡献更大。
加权中值估计应用。 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
------------------------------------------------------------------------------
mrmodalNOME 假设估计应用
. 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
专题：面板数据 regife：面板交互固定效应模型-Interactive Fixed Effect 专题：IV-GMM Stata：无需工具变量的IV估计-kinkyreg- Stata：当工具变量小于内生变量时，该如何估计？-mmeiv IV估计中的外部有效性(External Validity) IV在哪里？奇思妙想的工具变量 twostepweakiv：弱工具变量有多弱？ 多个(弱)工具变量如何应对-IV-mivreg？ IV：工具变量不满足外生性怎么办？ IV-工具变量法：第一阶段系数符号确定时的小样本无偏估计 IV：可以用内生变量的滞后项做工具变量吗？ Stata: 工具变量法 (IV) 也不难呀！ IV专题: 内生性检验与过度识别检验 IV 经典：寻找 IV 的足迹——Card（1993） IV-估计：工具变量不外生时也可以用！ 专题：内生性-因果推断 Abadie新作：简明IV,DID,RDD教程和综述 工具变量-IV：排他性约束及经典文献解读 专题：Python-R-Matlab Python爬虫：从SEC-EDGAR爬取股东治理数据-Shareholder-Activism 专题：Markdown-LaTeX LaTeX小白入门：TeXLive安装及使用
🍓 课程推荐：2024 空间计量专题
主讲老师：范巧 (兰州大学)
课程时间：2024 年 10 月 2-4 日 (三天)
课程咨询：王老师 18903405450（微信）
New！ Stata 搜索神器：
lianxh和songblGIF 动图介绍
搜： 推文、数据分享、期刊论文、重现代码 ……
👉 安装：
. ssc install lianxh
. ssc install songbl
👉 使用：
. lianxh DID 倍分法
. songbl all
🍏 关于我们
连享会 ( www.lianxh.cn，推文列表) 由中山大学连玉君老师团队创办，定期分享实证分析经验。 直通车： 👉【百度一下： 连享会】即可直达连享会主页。亦可进一步添加 「知乎」,「b 站」,「面板数据」,「公开课」 等关键词细化搜索。
