孟德尔随机化（Mendelian Randomization, MR）是一种利用遗传变异作为工具变量来探索潜在因果关系的统计方法。这种方法可以绕开传统观察性研究中的一些常见偏误，如混杂因素、反向因果等问题。在R语言中，MendelianRandomization包提供了一系列功能，以支持单变量和多变量MR分析的实施。
主要功能和方法
单变量MR方法
-
IVW (Inverse Variance Weighted)：逆方差加权方法，这是MR分析中的标准方法，适用于估计一个单一遗传工具变量与表型之间的关系。
- 逆方差加权（IVW）方法之所以被这样命名，主要是因为它在估计过程中对每个工具变量（例如遗传变异）的贡献进行加权，而这种加权是基于各个工具变量影响暴露变量的估计方差的逆。简单来说，每个工具变量的影响被赋予一个权重，该权重是其方差的倒数。这样做的目的是为了增强估计的精确度，特别是当不同的工具变量对于暴露变量的影响估计精度不一致时。
- 如果你正在研究吸烟（暴露变量）和肺癌（结果变量）之间的关系，IVW方法会使用各种与吸烟相关的遗传变异（如烟草使用量相关的SNP），并根据这些变异的精确度加权，从而估计吸烟对肺癌风险的总体影响。
-
Median-based：中位数方法，当工具变量可能存在无效时，提供一个更稳健的因果估计。
- 在研究维生素D水平（暴露变量）与骨折风险（结果变量）的关系时，如果部分SNP与其他疾病相关联，中位数基方法可以通过专注于中间值来提供一个较为稳健的估计。
-
MR-Egger：可以提供对工具变量多效性的敏感性分析，同时也允许存在工具变量无效的情况。
- 在探究胆固醇水平（暴露变量）对心血管疾病（结果变量）影响的研究中，如果一些胆固醇相关的SNP也影响其他病状（如糖尿病），MR-Egger可以帮助识别并纠正这种多效性的影响。
-
Maximum Likelihood：极大似然估计法，考虑了SNP与暴露变量关联的不确定性，并可以处理样本重叠的问题。
- 研究体重（暴露变量）对心脏病（结果变量）的影响时，如果使用的SNP数据来源于部分相同的研究群体，极大似然方法能够适当考虑这种样本重叠，提供更准确的因果关系估计。
选择使用孟德尔随机化（MR）中的哪种单变量方法取决于数据的特点、工具变量的质量以及所面对的具体统计挑战。下面是每种方法的使用原则和推荐的使用顺序：
1. 逆方差加权（IVW）
使用原则：
- 当所有的工具变量都认为是有效的，即它们仅通过影响暴露变量来影响结果变量，没有通过任何其他途径影响结果。
- 没有证据显示工具变量存在多效性（即一个遗传变异影响多个表型）。
优先使用条件：
- 当你拥有多个信誉良好、与暴露变量强相关的SNPs时。
- 当遗传工具变量与暴露之间的关联非常明确且不存在显著的遗传异质性时。
2. MR-Egger
使用原则：
- 用于检测和校正多效性偏差，即某些SNP可能通过不相关的途径影响结果。
- 当存在工具变量无效性的疑虑时，可以用来估计这种偏差的大小和方向。
优先使用条件：
- 当研究中的SNP可能不完全满足排除限制（仅通过暴露影响结果）时。
- 当需要评估和校正多效性偏差时，作为分析的一部分进行敏感性测试。
3. 中位数基方法（Median-based）
使用原则：
- 当存在一部分工具变量可能无效（即它们可能同时影响其他表型）时，提供一种更为稳健的因果估计方法。
- 当数据中存在异常值或者工具变量效应大小的分布极为不均时。
优先使用条件：
- 在对IVW结果的稳健性有疑虑时，作为一种补充分析。
- 当部分SNP的有效性受到质疑，或者存在显著的遗传异质性时。
4. 极大似然（Maximum Likelihood）
使用原则：
- 考虑到SNP与暴露之间关联的不确定性，同时可以处理样本重叠的问题。
- 当需要精细调整样本重叠对因果估计的潜在影响时。
优先使用条件：
- 当数据来自可能存在样本重叠的不同研究时，或者当单个研究中用于暴露和结果的样本部分相同时。
- 当分析需要精确控制遗传工具的不确定性和复杂性时。
总体策略
- 开始：通常以IVW方法作为首选，因为它直接、简单且在工具变量全部有效时非常强大。
- 敏感性分析：使用MR-Egger检查多效性和工具变量无效性的影响；同时，运用中位数方法来评估结果的稳健性。
- 复杂情境处理：在面对样本重叠或高度复杂的数据结构时，采用极大似然方法。
多变量MR方法
- MV-IVW：多变量逆方差加权，同时考虑多个遗传工具变量。
- MV-Egger：多变量Egger回归，可以在多变量情境下评估多效性和工具变量无效的影响。
里面还有多种方法，注意带有mr_m为多变量的，不带有的为单变量
image.png
应用示例
使用MendelianRandomization包进行MR分析的R代码示例：
# 加载R包
library(MendelianRandomization)
# 输入数据准备
MRInputObject <- mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse)
# 使用逆方差加权方法
IVWObject <- mr_ivw(MRInputObject)
# 使用中位数方法进行稳健分析
MedianObject <- mr_median(MRInputObject)
# 使用MR-Egger方法分析潜在的多效性
EggerObject <- mr_egger(MRInputObject)
# 使用极大似然法考虑SNP-exposure关联的不确定性
MaxLikObject <- mr_maxlik(MRInputObject)
# 多变量MR分析
MRMVInputObject <- mr_mvinput(bx = cbind(ldlc, hdlc, trig), bxse = cbind(ldlcse, hdlcse, trigse), by = chdlodds, byse = chdloddsse)
MRMVObject <- mr_mvivw(MRMVInputObject)
数据可视化
MendelianRandomization包还提供了数据可视化功能，可以用于生成散点图、森林图等，帮助用户评估MR分析的结果：
# 绘制基于IVW方法的MR分析结果图
mr_plot(MRInputObject, line="ivw")
# 绘制多方法比较的图形
mr_plot(mr_allmethods(MRInputObject, method="all"))
示例数据读取
用了 MendelianRandomization 包中的 extract.pheno.csv 函数，该函数用于提取和整理特定的遗传和表型数据，以便进行孟德尔随机化分析。
# 定义指向包内预设数据文件的路径，这些文件包含PhenoScanner数据库的SNP信息
path.noproxy <- system.file("extdata", "vitD_snps_PhenoScanner.csv",
package = "MendelianRandomization")
path.proxies <- system.file("extdata", "vitD_snps_PhenoScanner_proxies.csv",
package = "MendelianRandomization")
# 提取不使用代理SNPs的数据，用于暴露变量“log(eGFR creatinine)”和结果变量“Tanner stage”
# 这些数据针对的是欧洲血统的研究
extract.pheno.csv(
exposure = "log(eGFR creatinine)", # 指定暴露变量
pmidE = 26831199, # 暴露数据的PubMed ID
ancestryE = "European", # 暴露数据的人种背景
outcome = "Tanner stage", # 指定结果变量
pmidO = 24770850, # 结果数据的PubMed ID
ancestryO = "European", # 结果数据的人种背景
file = path.noproxy # 指定使用不包含代理的数据文件
)
# 提取使用代理SNPs的数据，同样针对暴露变量“log(eGFR creatinine)”和结果变量“Tanner stage”
extract.pheno.csv(
exposure = "log(eGFR creatinine)",
pmidE = 26831199,
ancestryE = "European",
outcome = "Tanner stage",
pmidO = 24770850,
ancestryO = "European",
rsq.proxy = 0.6, # 设置代理SNPs的R²阈值为0.6，表示中等关联强度
file = path.proxies # 指定使用包含代理的数据文件
)
# 提取使用代理SNPs的数据，暴露变量“log(eGFR creatinine)”和不同的结果变量“Asthma”
extract.pheno.csv(
exposure = "log(eGFR creatinine)",
pmidE = 26831199,
ancestryE = "European",
outcome = "Asthma", # 指定不同的结果变量为“哮喘”
pmidO = 20860503, # 哮喘数据的PubMed ID
ancestryO = "European",
rsq.proxy = 0.6, # 同样设置代理SNPs的R²阈值为0.6
file = path.proxies # 使用包含代理的数据文件
)
从PhenoScanner提供的示例数据集中提取了用于孟德尔随机化的必要遗传信息。其中，rsq.proxy参数指定了接受的代理SNP与主SNP之间的最低相关性阈值，用于在主SNP数据不足时扩充可用的遗传工具变量。这使得能够在数据可用性受限的情况下，继续进行有效的MR分析。
单变量方法
IVW方法
逆方差加权（IVW）方法基于一个关键假设：所有的遗传变异都是有效的工具变量，比如它们只通过影响暴露变量（如LDL cholesterol）来影响结果变量（如心血管疾病风险）。在这种假设下，可以将每个SNP的估计效应加权平均，以得到因果效应的整体估计。权重通常是每个估计的逆方差，以提高精确度较高的估计的贡献度。
# 加载MendelianRandomization包
library(MendelianRandomization)
# 创建输入对象，包括每个SNP对暴露和结果的效应大小及其标准误
MRInputObject <- mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse)
print(MRInputObject)
# 使用IVW方法进行估算，不使用稳健回归，不对SNP中的outliers进行惩罚
IVWObject1 <- mr_ivw(MRInputObject, model= "default", robust = FALSE, penalized = FALSE,
correl = FALSE, weights = "simple", psi = 0, distribution = "normal", alpha = 0.05)
print(IVWObject1)
# 使用稳健回归，并对SNP中的outliers进行惩罚，以提高估计的稳健性
IVWObject2 <- mr_ivw(MRInputObject, model= "default", robust = TRUE, penalized = TRUE,
correl = FALSE, weights = "simple", psi = 0, distribution = "normal", alpha = 0.05)
print(IVWObject2)
median-based方法
在孟德尔随机化分析中，中位数基方法（Median-based）是用来提供对个别SNP异常值或无效工具变量的稳健因果估计。这种方法可以有效地减少异常数据点对总体估计结果的影响。下面的代码片段展示了如何使用不同的加权策略来执行中位数基方法，并注释解释了每个函数调用的目的和参数的意义。
# 加载所需的MendelianRandomization包
library(MendelianRandomization)
# 创建输入数据对象，包括SNP对LDL胆固醇和冠心病风险的影响及其标准误
MRInputObject <- mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse)
# 使用加权的中位数方法进行因果估计，加权以反映每个估计的不确定性
WeightedMedianObject1 <- mr_median(MRInputObject, weighting = "weighted", distribution = "normal", alpha = 0.05, iterations = 10000, seed = 314159265)
print(WeightedMedianObject1) # 打印结果，显示LDL升高与增加的CHD风险之间的显著关联
# 使用“penalized”加权法，这种方法在考虑权重时对潜在的异常值或强影响点进行惩罚
WeightedMedianObject2 <- mr_median(MRInputObject, weighting = "penalized", distribution = "normal", alpha = 0.05, iterations = 10000, seed = 314159265)
print(WeightedMedianObject2) # 打印结果
# 使用简单中位数方法，不进行任何加权
WeightedMedianObject3 <- mr_median(MRInputObject, weighting = "simple", distribution = "normal", alpha = 0.05, iterations = 10000, seed = 314159265)
print(WeightedMedianObject3) # 打印结果
# 增加迭代次数到100,000，以提高估计的精度和稳健性
WeightedMedianObject4 <- mr_median(MRInputObject, weighting = "weighted", distribution = "normal", alpha = 0.05, iterations = 100000, seed = 314159265)
print(WeightedMedianObject4) # 打印结果
参数解释
- weighting: 控制如何加权每个SNP的估计。可选值包括"weighted"（考虑标准误的逆），"penalized"（对潜在的异常值施加惩罚），和"simple"（不加权）。
- distribution: 指定使用的分布类型，这里使用正态分布（"normal"）。
- alpha: 显著性水平，常设置为0.05。
- iterations: 迭代次数，用于计算的迭代过程中提高估计的精度。
- seed: 随机数种子，确保结果的可重复性。
示例展示了如何在不同场景下应用中位数方法来评估LDL对冠心病风险的影响，同时考虑到了数据中可能存在的异常值和统计的稳健性。通过调整加权策略和迭代次数，可以根据具体研究需求优化因果估计的准确性和稳健性。
MR-Egger方法
MR-Egger 方法在孟德尔随机化分析中被用来检测和调整工具变量的潜在多效性或偏倚。该方法类似于常规的回归分析，但增加了对工具变量多效性偏倚的估计和调整。这可以帮助研究者评估工具变量是否仅通过暴露变量影响结果，还是也通过其他路径影响结果，这种偏倚被称为多效性偏倚或遗传混杂。
# 加载所需的MendelianRandomization包，这是执行MR分析的必要工具包
library(MendelianRandomization)
# 创建输入数据对象，包括SNP对LDL胆固醇和冠心病风险的影响及其标准误
MRInputObject <- mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse)
# 执行MR-Egger回归分析，未应用稳健回归或惩罚方法
EggerObject1 <- mr_egger(MRInputObject, robust = FALSE, penalized = FALSE, correl = FALSE, distribution = "normal", alpha = 0.05)
print(EggerObject1) # 打印结果，提供不考虑异常值调整的基本MR-Egger估计
# 执行MR-Egger回归分析，应用稳健回归和惩罚方法
EggerObject2 <- mr_egger(MRInputObject, robust = TRUE, penalized = TRUE, correl = FALSE, distribution = "normal", alpha = 0.05)
print(EggerObject2) # 打印结果，提供考虑异常值调整后的稳健MR-Egger估计
# 分析和讨论结果
# 使用稳健回归后，MR估计值的误差变小了，但显著性并未改变。
# 这表明LDL胆固醇水平的升高与增加的冠心病风险之间的因果关联是稳健的，即使考虑潜在的多效性偏倚。
参数解释
- robust: 是否使用稳健回归，帮助减少异常值和影响点对回归分析结果的影响。
- penalized: 是否应用惩罚项，用于减少模型的过拟合和增强模型预测的稳定性。
-
correl: 这个参数用于指示是否考虑SNP之间的相关性，通常在MR-Egger分析中设为
FALSE。 - distribution: 指定模型估计的分布，通常为正态分布（"normal"）。
- alpha: 显著性检验的阈值，常设置为0.05。
Maximum likelihood方法（极大似然估计法)
极大似然方法（Maximum Likelihood, ML）在孟德尔随机化（Mendelian Randomization, MR）分析中提供了一个强有力的工具，尤其是在处理样本重叠和SNP与暴露之间关联不确定性的情况下。这个方法可以通过调整psi参数来适应不同程度的样本重叠，从而提供更为精确的因果推断。
极大似然方法的优势
- SNP-Exposure关联的不确定性：极大似然估计法允许研究者考虑到每个SNP与暴露之间关联强度的不确定性，这在简单的逆方差加权方法（IVW）中通常是被忽略的。
-
样本重叠情况：当使用来自相同样本的遗传和表型数据时，可能会出现样本重叠，这会影响到估计的准确性。极大似然方法通过
psi参数调整，可以有效管理这一点。
参数psi的作用
-
psi = 0：表示样本之间完全不重叠，即传统的独立双样本MR研究场景。 -
psi > 0：表示样本之间有重叠，psi的大小代表暴露和结果变量之间的相关性，即观察性研究中暴露和结局的相关系数。
这里展示了如何使用不同的psi值进行极大似然估计，以探索样本重叠对因果估计的影响：
library(MendelianRandomization)
# 极大似然估计，样本不重叠
MaxLikObject1 <- mr_maxlik(MRInputObject, model = "default", correl = FALSE, psi = 0, distribution = "normal", alpha = 0.05)
print(MaxLikObject1)
# 极大似然估计，样本轻微重叠，psi = 0.3
MaxLikObject2 <- mr_maxlik(MRInputObject, model = "default", correl = FALSE, psi = 0.3, distribution= "normal", alpha = 0.05)
print(MaxLikObject2)
# 极大似然估计，样本中等重叠，psi = 0.6
MaxLikObject3 <- mr_maxlik(MRInputObject, model = "default", correl = FALSE, psi = 0.6, distribution= "normal", alpha = 0.05)
print(MaxLikObject3)
# 极大似然估计，样本高度重叠，psi = 0.9
MaxLikObject4 <- mr_maxlik(MRInputObject, model = "default", correl = FALSE, psi = 0.9, distribution= "normal", alpha = 0.05)
print(MaxLikObject4)
通过这种方式，研究者可以评估样本重叠在MR分析中的实际影响，确保因果估计的可靠性和精确性。
多样本
MRMVInputObject <- mr_mvinput(bx = cbind(ldlc, hdlc, trig),
bxse = cbind(ldlcse, hdlcse, trigse),
by = chdlodds, byse = chdloddsse) #MVMR的input格式会和单变量的有所不同
MRMVInputObject
MRMVObject1 <- mr_mvivw(MRMVInputObject,
model = "default",
correl = FALSE,
distribution = "normal",
alpha = 0.05)
MRMVObject2 <- mr_mvegger(MRMVInputObject,
orientate = 1,
correl = FALSE,
distribution = "normal",
alpha = 0.05)
绘图
结果取决于传递给 mr_plot 的对象类型。当对象是 MRInput对象，该函数使用 plot 命令（如果 interactive 设置为 FALSE）或 plotly语法（如果 interactive 设置为 TRUE）来绘制关联估计值。当。。。的时候object 是一个 MRMVInput 对象，功能类似，除了我们绘制估计的关联结果在 y 轴上，关联的拟合值与来自x轴上的逆方差加权方法。如果 interactive 设置为 FALSE，则为静态图被生产。通过将标签设置为 TRUE，基因变体的名称出现在点上方。这会产生一个视觉上不太吸引人的图表，但更容易识别个人遗传变异。如果 interactive 设置为 TRUE，则绘图是交互式的，用户可以悬停在各个点上查看相关遗传变异的名称及其关联估计。当对象是 MRAll 对象时，该函数生成一个 ggplot 来比较因果估计通过不同的方法提出。
mr_plot(mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse),
line="egger", orientate = TRUE)
image.png
mr_plot(mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse),
line="ivw", interactive=FALSE) # produces a static graph
image.png
mr_plot(mr_allmethods(mr_input(bx = ldlc, bxse = ldlcse,
by = chdlodds, byse = chdloddsse), method="all", iterations = 50))
image.png
