关于MR-Egger方法的注意事项（1）
关于MR-Egger方法的注意事项（1）
知识回顾
在往期的内容中，我和大家简单介绍过MR研究中IVW和MR-Egger这两种方法的区别，具体参见孟德尔随机化之IVW和MR-Egger方法简介。
主要问题
前一阵子，有小伙伴提出来用“TwoSampleMR”包里的MR-Egger方法算出来的结果和直接在回归模型里添加截距项的结果不同。我看了一下，这里主要是因为暴露数据里的beta值存在负数，要想彻底理解这个问题，我们有必要看一下计算的源代码。
“TwoSampleMR”包里的MR-Egger计算代码如下（该代码可以在R语言中加载好TwoSampleMR包后直接输入mr_egger_regression并回车即可获取）：
function (b_exp, b_out, se_exp, se_out,parameters)
{
stopifnot(length(b_exp) == length(b_out))
stopifnot(length(se_exp) == length(se_out))
stopifnot(length(b_exp) == length(se_out))
#上面三行stopifnot是排错步骤，针对的是数据不完整的情况
nulllist <- list(b = NA, se = NA, pval = NA, nsnp = NA, b_i = NA,
se_i = NA, pval_i = NA, Q = NA, Q_df = NA, Q_pval = NA,
mod = NA, smod = NA, dat = NA)
if(sum(!is.na(b_exp) & !is.na(b_out) & !is.na(se_exp) &
!is.na(se_out)) < 3) { # MR-Egger的计算需要至少3个有效的IV，否则报错
return(nulllist)
}
sign0 <- function(x) {
x[x == 0] <- 1 # 区别于sign()函数的不同之处，把beta为0的变成1，基础函数中sign(0) = 0
return(sign(x))
}
to_flip <- sign0(b_exp) == -1 # 确定exposure中beta<0的要翻转
b_out = b_out * sign0(b_exp) # 对于要翻转的SNP，其outcome中的beta值要取相反数
b_exp = abs(b_exp) # 暴露的beta值全变为正值
dat<- data.frame(b_out = b_out, b_exp = b_exp, se_exp = se_exp,
se_out = se_out, flipped = to_flip)
mod<- lm(b_out ~ b_exp, weights = 1/se_out^2) # 回归分析，逆方差加权，保留截距项
smod <- summary(mod)
if(nrow(coefficients(smod)) > 1) {
b <- coefficients(smod)[2, 1] # 获取beta值
se <- coefficients(smod)[2, 2]/min(1, smod$sigma) # 获取se
pval <- 2 * pt(abs(b/se), length(b_exp) - 2, lower.tail = FALSE)
b_i <- coefficients(smod)[1, 1] # 获取截距项
se_i <- coefficients(smod)[1, 2]/min(1, smod$sigma) # 获取截距项的se
pval_i <- 2 * pt(abs(b_i/se_i), length(b_exp) - 2, lower.tail =FALSE) # 获取截距项的p值（采用的是T检验）
Q <- smod$sigma^2 * (length(b_exp) - 2) # 计算MR结果的异质性
Q_df <- length(b_exp) - 2 # 计算自由度
Q_pval <- pchisq(Q, Q_df, lower.tail = FALSE) # 使用卡方检验获取异质性p值
}
else {
warning("Collinearities in MR Egger, try LD pruning the exposurevariables.")
return(nulllist)
}
return(list(b = b, se = se, pval = pval, nsnp = length(b_exp),
b_i = b_i, se_i = se_i, pval_i = pval_i, Q = Q, Q_df = Q_df,
Q_pval = Q_pval, mod = smod, dat = dat))
}上面这个计算方法中，它借用了R里的sign()基础函数来重新定义了sign0()这个函数，其目的就是把beta.exposure为0的符号变为1（不过米老鼠觉得没有必要）。
接下来，我们看看这里最关键的部分“to_flip”，这一块就是借用之前新定义的sign0()函数来把beta.exposure为负值的调整为正值，相应的beta.outcome也会取一个相反数，这样就保证了每个SNP的effect allele对exposure的效应是一致的（也即都增加暴露的风险），同时不改变利用单个SNP进行Wald ratio估计的效果。
接下来我将以示例文件具体展示一下比较的过程（示例文件可以点击阅读原文下载demo1.csv即可）：
# 利用MR-Egger方法计算（有flip的过程）
library(TwoSampleMR)
data = read.csv('demo1.csv',header=T)
mr_egger_regression(b_exp =data$beta.exposure, b_out = data$beta.outcome, se_exp = data$se.exposure,se_out = data$se.outcome)其主要输出结果如下（注意mr_egger_regression是TwoSampleMR包里的函数）：
# 直接回归（包含截距项，但没有flip的过程）
summary(lm(beta.outcome ~ 1 +beta.exposure, weights = 1/se.outcome^2, data = data))其输出结果如下：
通过上述比较我们不难发现，两个结果是有差异的，其中第一个是MR-Egger方法计算出来的结果，第二个就是普通的线性回归。从数学上来看，MR-Egger在进行回归分析前把位于二、三象限的点通过中心（原点）对称的方式分别转化到四、一象限了（二象限对称到四象限，三象限对称到一象限），经过这么转换后SNP对exposure的效应方向也就一致了，也便于解释接下来的多效性问题。
因此，当你的暴露文件中beta值有负数的时候不要直接拟合，可以先翻转符号后再拟合，这样才是MR-Egger计算方法。希望大家注意！
当然，你也可以直接使用TwoSampleMR提供的MR-Egger方法，这是没有问题的。
PS: 之前有小伙伴提出无法从我的Gitee里下载MRinstruments包，米老鼠查看了一下，是我之前把这个MRInstruments库设成私密文件了，在这里向大家道个歉。现在我已经将其设置为公开了，大家可以自由下载（Gitee的下载速度比Github快很多！）。