我附上我自己做的R中的code,以及如何读取结果中的效应量和P值,若有错误请大家指正:
#mediation: R Package for Causal Mediation Analysis#
library(mediation)
library(sandwich)
tab<-function(x1x){
print(class(x1x))
print(table(x1x))
print(paste0("缺失值:",sum(is.na(x1x))))
}
#示例1#
data(jobs)
a=jobs
a=data.frame(a)
tab(a$control)#中介变量:2分类
tab(a$nonwhite)#自变量:2分类
tab(a$work1)#结局变量:2分类
#
a$nonwhite_1=as.numeric(a$nonwhite)-1
tab(a$nonwhite_1)
med.fit <- glm(control ~ nonwhite, data = a,family = binomial("probit"))
out.fit <- glm(work1 ~ control +nonwhite, data = a,family = binomial("probit"))
med23.out1 <- mediate(med.fit1, out.fit1, treat = "nonwhite", mediator = "control1",control.value = "white0",treat.value = "non.white1")
summary(med23.out1)
plot(med23.out1)
#读取结果中ACME_average的效应量和P值#
med23.out1[["d.avg"]]#ACME_aver_value
med23.out1[["d.avg.ci"]]#ACME_aver_ci
med23.out1[["d.avg.p"]]#ACME_aver_p
#读取结果中ADE_average的效应量和P值#
med23.out1[["z.avg"]]#ACME_aver_value
med23.out1[["z.avg.ci"]]#ACME_aver_ci
med23.out1[["z.avg.p"]]#ACME_aver_p
|