中介分析教程

如果有两个性状（A 与 B）同时与另一性状 C 关联，通过中介分析可以检验性状 C 是否通过性状 A 影响性状 B 或相反，本文将带大家进行一个简单的中介分析。

理论知识

原理

中介分析在流行病学研究中的重要性在于揭示暴露对结果的影响途径。中介分析通常用于评估暴露的影响被解释的程度，或者在没有假设中介变量时（也称为中间变量）解释的程度时，通过这种方式，可以定义暴露对结果的总影响，由一组给定的中介变量的暴露的影响（间接效应）以及这些相同的中介变量无法解释的暴露的影响（直接效应）

通常，中介变量可以解释的部分为间接效应，由这些中介变量无法解释的部分是直接效应，人们期望总效应可以分解为直接效应和间接效应，也就是说，假设暴露因素有 15% 的风险，其中直接效应的风险为 10%，则间接效应的风险为 5%，换句话说，三分之一的总效应可以被中介变量解释，其余三分之二的总效应被其他的途径解释

[scode type="yellow"]
一般情况下，只有在自变量对因变量具有影响（具有总效应）时才考虑检验中介效应，但是自变量对因变量没有影响不代表不存在中间变量，只是此时的影响叫做间接效应，中介效应都是间接效应，但间接效应不一定是中介效应。在实际生活中，会出现自变量与中间变量、中间变量与因变量之间的效应方向相反，从而使总效应为 0 的情况，这种情况叫做效应遮盖，也就是广义中介效应。
[/scode]

关于中介分析的综述已经有很多，这里列举两篇，一篇中文一篇英文，供参考：

1. 文献1：Mediation analysis in epidemiology: methods, interpretation and bias

2. 文献2：温忠麟,.张雷,侯杰泰,刘红云.中介效应检验程序及其应用

中介效应的估计

中介效应的计算特别简单，简述如下：

1. 首先估计总效应。总效应当为自变量对因变量的总体效应，一般为：

\[Y = \hat{c}X + e_1 \]

2. 接着估计自变量和中介变量的关联：

\[M = \hat{a}X + e_2 \]

3. 最后估计直接效应，也就是使用中介变量矫正后的自变量对因变量的效应，也称直接效应：

\[Y = \hat{c'} X + \hat{b}M + e_3 \]

4. 现在，总效应为 \(c\)，直接效应为 \(c'\)，则中介效应为 \(c - c'\)，中介效应的占比为 \(\frac {c - c'} {c}\)

中介效应的检验

相比中介效应的估计，中介效应的检验要稍微复杂一些，Sobel 检验比较常用

下文检验方法中的参数均来自下图：

中介分析的实战

了解了中介分析的原理以及检验的统计学方法后，以一个实例来说明中介分析的流程。

设计模拟数据

使用 R 自带的鸢尾花的数据进行模拟研究，这里我们定义萼片长度 Sepal.Length 为自变量，命名为 X，定义其对蜜蜂的吸引程度 为中介变量，命名为 M，我们再定义一个因变量，我们假设其含义为最终被蜜蜂采蜜的可能性 ，命名为 Y，为了制造一些噪音和中介效应，我们定义这三个变量有以下的关系：

\[M = 0.35X + e_1\\ Y = 0.35M + e_2 \]

如此一来，我们的预期中介效应就为 12.25%，接着我们来进行检验

创建模拟数据

df=iris
set.seed(12334)

colnames(df)[1] = "X"

df$e1 = runif(nrow(df), min = min(df$X), max = max(df$X))
df$M = df$X * 0.35 + df$e1 * 0.65

df$e2 = runif(nrow(df), min = min(df$M), max = max(df$M))
df$Y = df$M * 0.35 + df$e2 * 0.65

检验总效应

检验的方程为：

\[Y = \hat{c}X + e_1 \]

检验的代码为：

fit_total = lm(Y ~ X, df)
summary(fit_total)

输出如下：

Call:
lm(formula = Y ~ X, data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.15930 -0.45815 -0.01242  0.44662  1.20905 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  5.29106    0.32791  16.136   <2e-16 ***
X            0.12984    0.05557   2.337   0.0208 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.5616 on 148 degrees of freedom
Multiple R-squared:  0.03558,	Adjusted R-squared:  0.02907 
F-statistic:  5.46 on 1 and 148 DF,  p-value: 0.02079

可以看到总效应为 0.12984，接近我们的预期 12.25%，p 值为 0.0208，第一步检验是显著的，接着进行下面的步骤

检验因变量对中介变量的效应

检验的方程为：

\[M = \hat{a}X + e_2 \]

检验的代码为：

fit_mediator = lm(M ~ X, df)
summary(fit_mediator)

输出如下：

Call:
lm(formula = M ~ X, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.2494 -0.5082  0.0123  0.5483  1.0799 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4.32300    0.37901  11.406  < 2e-16 ***
X            0.30429    0.06422   4.738 5.02e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.6492 on 148 degrees of freedom
Multiple R-squared:  0.1317,	Adjusted R-squared:  0.1258 
F-statistic: 22.45 on 1 and 148 DF,  p-value: 5.019e-06

可以看到中介变量对因变量的效应为 0.30429，接近我们的预期 35%，p 值为 5.019e-06，第二步检验也是显著的，接着进行下面的步骤

检验直接效应

检验的方程为：

\[Y = \hat{c'}X + \hat{b}M + e_3 \]

检验的代码为：

fit_direct=lm(Y ~ X + M, df)
summary(fit_direct)

输出如下：

Call:
lm(formula = Y ~ X + M, data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.90734 -0.46316  0.00764  0.39751  0.87026 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.68314    0.40721   9.045 8.17e-16 ***
X            0.01667    0.05402   0.309    0.758    
M            0.37194    0.06443   5.773 4.46e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.5088 on 147 degrees of freedom
Multiple R-squared:  0.2138,	Adjusted R-squared:  0.2031 
F-statistic: 19.99 on 2 and 147 DF,  p-value: 2.092e-08

可以看到使用中介变量进行矫正后，直接效应为 0.016，很微弱，p 值为 0.758，不显著，而中介变量对因变量的效应值为 0.37194，接近我们的预期 35%，而 p 值为 4.46e-08，这个结果说明使用中介变量进行校正后，自变量对于因变量的影响没有了显著性，这种情况称为完全中介

sobel检验

如果系数a和系数b不都显著的话，则由以下公式计算中介效应的显著性，这个方法由Sobel提出的，所以称为Sobel检验：

\[T = \frac {\hat{a}\hat{b}} {\sqrt{\hat{a}^2s_b^2 + \hat{b}^2s_a^2}} \]

其中，\(\hat{a}\) 是自变量对中介变量的效应，\(\hat{b}\) 是使用中介变量进行校正后的中介变量对因变量的效应，\(s_a\) 和 \(s_b\) 分别是 \(\hat{a}\) 和 \(\hat{b}\) 的标准误，在上面的结果中指的是 Std. Error 这一列，在 R 中操作如下：

a = coefficients(fit_mediator)["X"]
s_a = summary(fit_mediator)$coefficients["X", "Std. Error"]
b = coefficients(fit_mediator)["M"]
s_b = summary(fit_direct)$coefficients["M", "Std. Error"]
SE = sqrt(a^2 * s_b^2 + b^2 * s_a^2)
# 自由度为n-k-1，k是解释变量的个数，这里为2，因此df = 147
t_statistic = a * b / SE
p = 2 * pt(
	abs(t_statistic),
	df = df.residual(fit_mediator),
	lower.tail = FALSE
)

总结

在上面的过程中，自变量对于因变量的总体效应为 0.12984，自变量对于中介变量的影响系数、中介变量对因变量的影响系数都接近我们模拟的 0.35，而在矫正了中介变量的影响后，自变量对于因变量的影响（直接效应）变得很小（0.01667）且不显著，因此，我们可以说自变量对于因变量的影响是完全通过中介变量进行 的，称为完全中介效应

由于总效应是 0.12984，直接效应为 0.01667，所以此时的中介效应为 0.12984 - 0.01667 = 0.11317，中介效应占比 87.16%

使用 R 包简化操作

上面的操作略显繁琐，那么有没有 R 包封装了这一过程呢？当然有，那便是 mediation ：

install.packages("mediation")
library(mediation)

results = mediate(
	model.m  = fit_mediator,
	model.y  = fit_direct,
	treat    = "X",
	mediator = "M",
	boot     = TRUE
)
summary(results)

结果如下：

Causal Mediation Analysis 

Nonparametric Bootstrap Confidence Intervals with the Percentile Method

               Estimate 95% CI Lower 95% CI Upper p-value    
ACME             0.1132       0.0575         0.17  <2e-16 ***
ADE              0.0167      -0.0904         0.12   0.752    
Total Effect     0.1298       0.0152         0.24   0.018 *  
Prop. Mediated   0.8716       0.3796         3.68   0.018 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Sample Size Used: 150 

Simulations: 1000

这里结果中的 ACME 代表平均中介效应，等于 \(\hat{a}*\hat{b}\)，ADE 代表直接效应，Total Effect 是总效应，Prop. Mediated 是中介效应占比

一般将自变量对中介变量的效应 \(\beta_1\) 和中介变量对因变量的效应 \(\beta_2\) 的乘积当作平均中介效应

采用中介分析的研究

以下研究使用了中介分析

数据类型	URL	中介方法	说明
eQTL-meQTL	nature communications	Sobel	对每一个共定位的数据对重新进行了关联分析并检验中介效应
eQTL-pQTL	nature	孟德尔随机化	遗传变异作为工具变量，血浆蛋白作为暴露（即pQTL），疾病作为 outcome
pQTL-eQTL	nature communications	孟德尔随机化	使用 GWAS 的独立位点（0.1）作为工具变量，CHD 作为 outcome，蛋白质作为暴露（pQTL），分别用了多位点 MR 和 wald 检验计算多独立位点和单独立位点的情况的显著性，使用工具为 MRbase