Brief summary

网上关于中介效应的资料挺多的，这里就不对原理进行过多解释。

简单来说，我们分析自变量 X 对因变量 Y 产生的影响，如果变量 X 通过影响变量 M 来影响变量 Y ，那么这个变量 M 就是中介变量。
例如租客 (X) 通过中介公司 (M) 找到合适的房子 (Y)，中介公司就在其中扮演了中介变量的角色，中介变量发挥的作用就称为中介效应了。

对于该分析我个人有一个认知，就是X和Y必须有关联。这个关联可以是数据上的（correlation or association，calculated from statistical method）或者是学科上的（因为我的研究方向是疾病多组学，可以说是生物学上的关联）。
如果X和Y毫无关系，这个研究是否还有必要呢，这可能是值得思考的

这里分享几个可以帮助理解的链接：
https://zhuanlan.zhihu.com/p/376007591
https://www.jianshu.com/p/bdf084f42e80
https://www.docin.com/p-1622412323.html

---

中介效应的检验

参考来源中介效应检验程序及其应用，温忠麟，2004
这篇文章很好的介绍了中介效应的概念，站里好像就有，感兴趣的可以看一下~

分析方法

Mediation包

在网络上能看到最多的应该就是用mediation包进行的中介分析，简单好用

df<- read.csv('~/PD/PD分期/olink_sig_original.csv', check.names = F)
df[is.na(df)]<- 0 #补充NA数据，这一步看自我需求
#因子化我的分类变量
df$Gender<- as.factor(df$Gender)
df$Group<- as.factor(df$Group)

建立线性回归，我的X是CX3CL1, M是Meta52, Y是Group（二分类变量）

#建立线性回归，我的X是CX3CL1, M是Meta52, Y是Group（二分类变量）
a<- lm(Meta52 ~ CX3CL1+Age+Gender,df) #lm(M~X,df)
b<- glm(Group ~Meta52+CX3CL1+Age+Gender, df,family = binomial) #glm(Y~X+M)
#install.packages("mediation")
library(mediation)
set.seed(123) #保证结果可以复现
result = mediate(a,b,treat="CX3CL1",mediator = "Meta52",boot = T)#默认1000次抽样
#这里需要注意变量的名称，如果名称比较特殊，有中文或者符号很容易出错，建议提前处理变量名
summary(result)

重点关注的结果：
ACME(average）间接效应
Prop.Mediated (average) 中介效应占比
可以从间接效应的p-value看出来，且置信区间没有0，这个中介模型是成立的。中介变量占该模型94%的效应 （一般不会那么高，我的数据可能存在过拟合的可能性了）

---

BruceR包

关于这个包的中介模型用法详情可以参考下面这个知乎文章（也介绍了调节效应、链式中介等），真的救我于水火之中
https://zhuanlan.zhihu.com/p/376007591

如果变量不多，或者研究的变量关系很明确，可以用mediation包跑。但我面临的数据是十几个X变量，100多个M中介变量，粗略下来大概1000+个中介模型，真的是跑断腿也搞不定。虽然我有尝试写mediate（）的循环，但怎么也解决不了treat和mediator要固定输入变量名的问题。

还好上帝开了一扇新的窗，让我发现了BruceR这个神包。这个包内嵌很多个常用的包（ggplot2、lmerTest、dplyr等等等）最厉害的是用来做中介效应/调节效应的PROCESS函数，自动化为你解决该分析，不用自己建立回归方程，只需要准备好一份数据集

#继续用上面的数据
#使用方法
PROCESS(data, y = "Group", x = "CX3CL1",
        meds = "Meta52", #中介变量，可以是多个，用c()
        covs = c("Age","Gender"), #协变量
        ci="boot", nsim=1000, seed =1)
        
#写循环
X.names<- colnames(df1[5:18]) #所有X变量的名字
M.names<- colnames(df1[19:111]) #所有M变量的名字
result1<- list()#建一个空的list放结果
k = 1
for (i in X.names){
  for (j in M.names) {
    k = 1
    result1[[i]] <- list()
    result1[[i]][[k]] <- PROCESS(data, y = "Group", x = i,
                                 meds = j,
                                 covs = c("Age","Gender"),
                                 ci="boot", nsim=1000, seed =1)
    k = k +1
  }
}

结果解读

Part1

************ PART 1. Regression Model Summary ************

Model Summary

图一列出了建立回归所用的变量和formula信息
图二，不带括号的是斜率b，括号内的数字是标准差SE，星号代表显著

Part2

PART 2. Mediation/Moderation Effect Estimate

这是最重要最需要看的部分了
通过indirect的P-value<0.05 且CI不包括0，可以知道间接效应成立
中介效应=indirect effect / total effect = （-0.211)/(-0.222) = 95%
是一个完全中介效应

结果整理（for linux）

因为我的数据计算量很大，运程是在linux服务器的R跑的，所以对结果的提取是perl + shell scirpt。
附上在shell运行的完整代码, 使用的时候请自行更改

将下面代码储存为一个.R 文件，假设命名为MedAnalysis.R

library(bruceR)
library(texreg)
data<- read.csv('your file.csv')
data$Group<- as.factor(data$Group)
data$Gender<- as.factor(data$Gender)

micro.names<- colnames(df1[12:18])
meta<- colnames(data[20:128])

result1<- list()
k = 1
for (i in micro.names){
  for (j in meta) {
    k = 1
    result1[[i]] <- list()
    result1[[i]][[k]] <- PROCESS(data, y = "Group", x = i,
                                 meds = j,
                                 covs = c("Age","Gender"),
                                 ci="boot", nsim=5000, seed =1)
    k = k +1
  }
}

同时我们准备好跑完后用于提取关键信息（如简介效应pvalue<0.05的）的perl代码，同样存成pl脚本文件。假设命名为 sel.pl

#!/usr/bin/perl -w
use strict;

open IN,$ARGV[0] or die $!;
while(<IN>)
{
	chomp;
	my $line1=<IN>;
	chomp $line1;
	my $line2=<IN>;
	chomp $line2;
	if($line1=~/\*/)
	{
		print "$_\n$line1\n$line2\n";
	}
}
close IN;

准备好R脚本和perl脚本，可以开始后台运行了

1. nohup Rscript MedAnalysis.R > result.out &  #放到后台跑，等结果
2. less result.out |grep "irect" result.out >result_sel.out#在如上面图一图二一大串结果中，只提取part2 有用的pvalue信息
3. perl sel.pl result_sel.out >result_sel_res.out #只保留indirect显著的组合

R跑完产生的.out文件， 每组都是重复的这些，很多不需要的信息

提取后的result_sel.out 文件

提取后的result_sel_res.out文件，只有indirect <0.05的组合

最后就根据结果，进一步的筛选重要的中介模型了

————————————————————————
感觉整理的不是很清晰完美，有机会再修改修改
DONE