中介分析教程

如果有兩個性狀（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