x为多分类变量，m、y为连续变量，计算相对中介效应
image.png
image.png
image.png
image.png
rm(list = ls())
#install.packages("mediation")
library(mediation)
data(jobs)
head(jobs)
# 自变量为类别变量，中介变量和因变量为连续变量的中介分析，同样都是拟合线性模型。
# 如果自变量为二分类别变量，可利用回归分析按照逐步法进行中介分析
# 如果自变量为k个类别(k≥3)建议使用相对中介和整体中介分析方法。
# 相对中介分析过程：
# 根据研究目的选择自变量的某个水平为参照水平，因为自变量的其它k-1个水平都要与参照水平进行对照，从而得到相对于参照水平的中介效应、直接效应和总效应；因此需要对自变量进行虚拟编码。
# 相对总效应的的大小等于k-1个相对直接效应和相对中介效应的大小之和。
# 需要说明的是，类别自变量的水平数越多，需要检验的相对中介效应、直接效应和总效应就越多，这可能会导致检验的第一类错误率增大，Bootstrap置信区间的置信度(用1–α/(k-1) 代替通常的 1–α)来控制检验的第一类错误率。
# 相对中介分析检验的是多类别自变量的某一个水平相对于参照水平而言是否存在显著的中介效应，整体中介分析检验的是k-1 个相对中介效应是否全部为 0。
# 整体中介分析相比相对中介分析存在至少三个优势:
# 第一，整体中介分析的结果不受参照水平的影响;
# 第二，整体中介分析的统计功效高于相对中介分析;
# 第三，整体中介分析得到的某些模糊结论(存在若干显著的相对中介效应)要比相对中介分析的某些结论(自变量的某一个水平相对于参照水平的中介效应不显著) 更有意义。
整体中介分析包括整体中介效应、整体直接效应和整体总效应三部分：
# 多类别自变量的中介分析流程
# 1）整体中介分析。如果整体中介效应不显著，就表示k-1 个相对中介效应全部为 0，则分析结束；否则进入步骤2
# 2）相对中介分析，弄清具体哪一个或哪一些相对中介效应显著；如果相对中介效应不显著，则分析结束；否则进入步骤3
# 3）报告相应的相对直接效应检验的显著性结果，即使直接效应不显著，也避免使用完全中介的概念。
# 仍然depress1作为自变量，econ_hard作为中介变量，depress2作为因变量；treat作为控制变量
# depress1转换成因子变量
jobs$depress1<-ifelse(jobs$depress1<1.5,
"d1",
ifelse(jobs$depress1>=1.5 & jobs$depress1<2,
"d2",
ifelse(jobs$depress1>=2 & jobs$depress1<2.5,
"d3","d4")))
jobs$depress1<-factor(jobs$depress1,levels = c("d1","d2","d3","d4"))
table(jobs$depress1)
第一个模型，X与Y的总效应
model.0 <- lm(depress2 ~ depress1+treat, data=jobs)
summary(model.0)
第二个模型，X与M的效应
model.1 <- lm(econ_hard~ depress1+treat, data=jobs)
summary(model.1)
第三个模型，M与Y的效应，调整X和Z
model.2 <- lm( depress2~ econ_hard+depress1+treat, data=jobs)
summary(model.2)
第四个模型，M与Y的模型，不调整X
model.3 <- lm( depress2~ econ_hard+treat, data=jobs)
summary(model.3)
编写一个函数，计算整体总效应、整体直接效应和整体中介效应是否有意义
Calculate_effect <- function(model0, model1, model2, model3, n,k,data){ #n为样本量,k为自变量分类数,data为数据集
#计算总体总效应检验
F_value <- (n-k)*summary(model0)$r.squared/(k*(1-summary(model0)$r.squared))
df1 <- k-1
df2 <- n-k
p_value1 <- 1 - pf(F_value, df1, df2)
print(paste0("总效应检验：P=", p_value1))
#计算总体的直接效应检验
delta_R2 <- summary(model2)$r.squared - summary(model3)$r.squared
F_value <- (n-k-1)*delta_R2/((k-1)*(1-summary(model2)$r.squared))
df1 <- k-1
df2 <- n-k-1
p_value2 <- 1 - pf(F_value, df1, df2)
print(paste0("总直接效应检验：P=", p_value1))
#计算中介效应检验
a <- summary(model.2)
b <- as.data.frame(a$coefficients)
M.effect <- summary(model1)$adj.r.squared*b$Estimate[2]
#借助bootstrap方法计算中介效应
library(boot)
# 定义计算 R^2*b统计量的函数
r_squared <- function(data, indices) {
# 使用 bootstrap 抽样的索引
d <- data[indices, ] # 选择样本
model.1 <- update(model1, data=d)# 拟合模型
model.2 <- update(model2, data=d)
return(summary(model.1)$adj.r.squared * coef(model.2)[2]) # 返回统计量的值
}
# 使用自助法计算
set.seed(123) # 设置随机种子以便复现结果
boot_results <- boot(data, statistic = r_squared, R = 500) # 500 次抽样
# 计算置信区间
out <- boot.ci(boot_results, type = "perc") # 百分位数置信区间
print(paste0("中介效应：", M.effect," 95%CI:",out$percent[4]," - ",out$percent[5]))
}
Calculate_effect(model.0,model.1,model.2,model.3,899,4,jobs)
使用’mediation’包中的mediate()进行相对中介效应分析
# 以b1作为基线水平，分别计算b2、b3、b4相对于b1的相对中介效应分析，为了防止检验的第一类错误增大，增加Bootstrap置信区间的置信度(用0.98代替通常的0.95)来控制检验的第一类错误率。
#b2的相对中介效应分析
contcont2 <- mediate(model.1,
model.2,
sims=500,
treat="depress1",
mediator="econ_hard",
boot=TRUE,
control.value="d1", ##b1作为基线水平
treat.value="d2",
conf.level = 0.98) #改！
summary(contcont2)
#b3的相对中介效应分析
contcont3 <- mediate(model.1,
model.2,
sims=500,
treat="depress1",
mediator="econ_hard",
boot=TRUE,
control.value="d1",
treat.value="d3",
conf.level = 0.98)
summary(contcont3)
#b4的相对中介效应分析
contcont4 <- mediate(model.1,
model.2,
sims=500,
treat="depress1",
mediator="econ_hard",
boot=TRUE,
control.value="d1",
treat.value="d4",
conf.level = 0.98)
summary(contcont4)
# 三个相对中介效应输出结果表明：3个相对中介效应都为部分中介效应。
