目录
- Brief summary
- 中介效应的检验
- 分析方法
- Mediation包
- BruceR包
- 结果解读
- Part1
- Part2
- 结果整理（for linux）
Brief summary
网上关于中介效应的资料挺多的，这里就不对原理进行过多解释。
简单来说，我们分析自变量 X 对因变量 Y 产生的影响，如果变量 X 通过影响变量 M 来影响变量 Y ，那么这个变量 M 就是中介变量。
例如租客 (X) 通过中介公司 (M) 找到合适的房子 (Y)，中介公司就在其中扮演了中介变量的角色，中介变量发挥的作用就称为中介效应了。
对于该分析我个人有一个认知，就是X和Y必须有关联。这个关联可以是数据上的（correlation or association，calculated from statistical method）或者是学科上的（因为我的研究方向是疾病多组学，可以说是生物学上的关联）。
如果X和Y毫无关系，这个研究是否还有必要呢，这可能是值得思考的
中介效应的检验
参考来源中介效应检验程序及其应用，温忠麟，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包
如果变量不多，或者研究的变量关系很明确，可以用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
