1. 中介分析
中介分析/介导分析 (Mediation analysis)：假设自变量X和因变量Y中间存在一个中介变量(Mediator)M，如果X对M有显著影响，同时M对Y有显著影响，那么就可以称在X和Y之间存在由M引起的中介效应。这个分析过程就叫中介分析。
如果X对Y也有意义，那么就是部分中介，否则就是完全中介。
2. mediation包
mediation包主要用来做中介分析，其核心是mediate函数。
mediate(model.m, model.y, sims = 1000, boot = FALSE,
boot.ci.type = "perc", treat = "treat.name", mediator = "med.name",
covariates = NULL, outcome = NULL, control = NULL,
conf.level = 0.95, control.value = 0, treat.value = 1,
long = TRUE, dropobs = FALSE, robustSE = FALSE, cluster = NULL,
group.out = NULL, use_speed = FALSE, ...)
参数说明（?mediate）
model.m
a fitted model object for mediator. Can be of class 'lm', 'polr', 'bayespolr', 'glm', 'bayesglm', 'gam', 'rq', 'survreg', or 'merMod'.
model.y
a fitted model object for outcome. Can be of class 'lm', 'polr', 'bayespolr', 'glm', 'bayesglm', 'gam', 'vglm', 'rq', 'survreg', or 'merMod'.
sims
number of Monte Carlo draws for nonparametric bootstrap or quasi-Bayesian approximation.
boot
a logical value. if 'FALSE' a quasi-Bayesian approximation is used for confidence intervals; if 'TRUE' nonparametric bootstrap will be used. Default is 'FALSE'.
boot.ci.type
a character string indicating the type of bootstrap confidence intervals. If "bca" and boot = TRUE, bias-corrected and accelerated (BCa) confidence intervals will be estimated. If "perc" and boot = TRUE, percentile confidence intervals will be estimated. Default is "perc".
treat
a character string indicating the name of the treatment variable used in the models. The treatment can be either binary (integer or a two-valued factor) or continuous (numeric).
mediator
a character string indicating the name of the mediator variable used in the models.
covariates
a list or data frame containing values for a subset of the pre-treatment covariates in 'model.m' and 'model.y'. If provided, the function will return the estimates conditional on those covariate values.
outcome
a character string indicating the name of the outcome variable in ‘model.y’. Only necessary if 'model.y' is of class 'survreg'; otherwise ignored.
control
a character string indicating the name of the control group indicator. Only relevant if 'model.y' is of class 'gam'. If provided, 'd0', 'z0' and 'n0' are allowed to differ from 'd1', 'z1' and 'n1', respectively.
conf.level
level of the returned two-sided confidence intervals. Default is to return the 2.5 and 97.5 percentiles of the simulated quantities.
control.value
value of the treatment variable used as the control condition. Default is 0.
treat.value
value of the treatment variable used as the treatment condition. Default is 1.
long
a logical value. If 'TRUE', the output will contain the entire sets of simulation draws of the the average causal mediation effects, direct effects, proportions mediated, and total effect. Default is 'TRUE'.
dropobs
a logical value indicating the behavior when the model frames of 'model.m' and 'model.y' (and the 'cluster' variable if included) are composed of different observations. If 'TRUE', models will be re-fitted using common data rows. If 'FALSE', error is returned. Default is 'FALSE'.
robustSE
a logical value. If 'TRUE', heteroskedasticity-consistent standard errors will be used in quasi-Bayesian simulations. Ignored if 'boot' is 'TRUE' or neither 'model.m' nor 'model.y' has a method for vcovHC in the sandwich package. Default is 'FALSE'.
cluster
a variable indicating clusters for standard errors. Note that this should be a vector of cluster indicators itself, not a character string for the name of the variable.
group.out
a character string indicating the name of the lmer/glmer group on which the mediate output is based. Can be used even when a merMod function is applied to only one of the mediator or the outcome. If merMod functions are applied to both the mediator and the outcome, default is the group name used in the outcome model; if the mediator group and the outcome group are different and the user is interested in the mediate output based on the mediator group, then set group.out to the group name used in the mediator merMod model. If a merMod function is applied to only one of the mediator or the outcome, group.out is automatically set to the group name used in the merMod model.
use_speed
a logical value indicating whether, if nonparametric bootstrap is used, lm and glm models should be re-fit using functions from the speedglm package. Ignored if 'boot' is 'FALSE' or if neither 'model.m' nor 'model.y' is of class 'lm' or 'glm'. Default is 'FALSE'.
...
other arguments passed to vcovHC in the sandwich package: typically the 'type' argument, which is ignored if 'robustSE' is 'FALSE'. Arguments to the boot in the boot package may also be passed, e.g. 'parallel' and 'ncpus'.
假设自变量X和因变量Y之间有中介变量M，mediate 函数需要两个统计模型 (X可以是一系列变量)：
- M对X的模型(model.m, outcome model, out.fit)
- Y对X的模型(model.y, mediator model, med.fit)
下图是可以放在mediation函数里的模型，*表示使用mediate函数的时候需要bootstrap非参，即设置boot = TRUE
Types of statistical models that can be used with the mediate function.
3. 演示
# **For illustration purposes a small number of simulations are used**
library(mediation)
data(jobs)
View(jobs)
Example 1: Linear Outcome and Mediator Models
#job_seek是mediator，b是x对m的模型，c是x+m对y的模型
b <- lm(job_seek ~ treat + econ_hard + sex + age, data=jobs)
c <- lm(depress2 ~ treat + job_seek + econ_hard + sex + age, data=jobs)
# Estimation via quasi-Bayesian approximation
contcont <- mediate(b, c, sims=50, treat="treat", mediator="job_seek")
summary(contcont)
# Causal Mediation Analysis
# Quasi-Bayesian Confidence Intervals
# Estimate 95% CI Lower 95% CI Upper p-value
# ACME -0.0181 -0.0355 0.01 0.20
# ADE -0.0324 -0.1107 0.04 0.52
# Total Effect -0.0505 -0.1277 0.04 0.24
# Prop. Mediated 0.3043 -1.0818 5.45 0.28
# Sample Size Used: 899
# Simulations: 50
plot(contcont)
结果说明
ACME: stands for average causal mediation effects.间接因果效应，表示X通过M对Y的效应大小
通过med.sum$d0和med.sum$d0.p可以获得ACME的效应和p值
ADE: stands for average direct effects.直接效应，表示X直接对Y的作用大小
通过med.sum$z0和med.sum$z0.p可以获得ADE的效应和p值
Total Effect: stands for the total effect (direct + indirect) of the IV on the DV. X对Y的直接和间接作用总和
Prop. Mediated: describes the proportion of the effect of the IV on the DV that goes through the mediator. X通过M对Y的作用的比例
# Estimation via nonparametric bootstrap (设置boot=TRUE)
contcont.boot <- mediate(b, c, boot=TRUE, sims=50, treat="treat", mediator="job_seek")
summary(contcont.boot)
# Causal Mediation Analysis
# Nonparametric Bootstrap Confidence Intervals with the Percentile Method
# Estimate 95% CI Lower 95% CI Upper p-value
# ACME -0.0157 -0.0466 0.00 0.28
# ADE -0.0403 -0.1239 0.03 0.24
# Total Effect -0.0560 -0.1373 0.02 0.16
# Prop. Mediated 0.2811 -0.1447 1.38 0.36
# Sample Size Used: 899
# Simulations: 50
# Allowing treatment-mediator interaction
d <- lm(depress2 ~ treat + job_seek + treat:job_seek + econ_hard + sex + age, data=jobs)
contcont.int <- mediate(b, d, sims=50, treat="treat", mediator="job_seek")
summary(contcont.int)
# Causal Mediation Analysis
# Quasi-Bayesian Confidence Intervals
# Estimate 95% CI Lower 95% CI Upper p-value
# ACME (control) -0.0154 -0.0420 0.01 0.36
# ACME (treated) -0.0111 -0.0321 0.01 0.36
# ADE (control) -0.0387 -0.1127 0.06 0.28
# ADE (treated) -0.0344 -0.1086 0.06 0.28
# Total Effect -0.0498 -0.1189 0.03 0.16
# Prop. Mediated (control) 0.1914 -1.9203 4.12 0.52
# Prop. Mediated (treated) 0.1432 -1.2864 3.24 0.52
# ACME (average) -0.0133 -0.0391 0.01 0.36
# ADE (average) -0.0365 -0.1115 0.06 0.28
# Prop. Mediated (average) 0.1673 -1.4223 3.65 0.52
# Sample Size Used: 899
# Simulations: 50
# Allowing ``moderated mediation'' with respect to age
b.int <- lm(job_seek ~ treat*age + econ_hard + sex, data=jobs)
d.int <- lm(depress2 ~ treat*job_seek*age + econ_hard + sex, data=jobs)
contcont.age20 <- mediate(b.int, d.int, sims=50, treat="treat", mediator="job_seek",
covariates = list(age = 20))
contcont.age70 <- mediate(b.int, d.int, sims=50, treat="treat", mediator="job_seek",
covariates = list(age = 70))
summary(contcont.age20)
# Causal Mediation Analysis
# Quasi-Bayesian Confidence Intervals
# (Inference Conditional on the Covariate Values Specified in `covariates')
# Estimate 95% CI Lower 95% CI Upper p-value
# ACME (control) -0.0456 -0.1281 0.01 0.28
# ACME (treated) -0.0321 -0.0801 0.02 0.28
# ADE (control) -0.0463 -0.1792 0.13 0.68
# ADE (treated) -0.0328 -0.1561 0.12 0.68
# Total Effect -0.0784 -0.2333 0.07 0.48
# Prop. Mediated (control) 0.2520 -3.0706 6.01 0.60
# Prop. Mediated (treated) 0.1591 -2.8914 3.62 0.60
# ACME (average) -0.0389 -0.0971 0.01 0.28
# ADE (average) -0.0395 -0.1667 0.12 0.68
# Prop. Mediated (average) 0.2055 -2.9062 5.30 0.60
# Sample Size Used: 899
# Simulations: 50
summary(contcont.age70)
# Causal Mediation Analysis
# Quasi-Bayesian Confidence Intervals
# (Inference Conditional on the Covariate Values Specified in `covariates')
# Estimate 95% CI Lower 95% CI Upper p-value
# ACME (control) 0.0151 -0.0930 0.13 0.92
# ACME (treated) 0.0152 -0.0563 0.08 0.80
# ADE (control) -0.0712 -0.2875 0.15 0.64
# ADE (treated) -0.0711 -0.3103 0.13 0.72
# Total Effect -0.0560 -0.2937 0.16 0.64
# Prop. Mediated (control) 0.0310 -4.8853 2.83 0.92
# Prop. Mediated (treated) -0.0674 -3.4487 1.61 0.88
# ACME (average) 0.0152 -0.0725 0.10 0.72
# ADE (average) -0.0711 -0.2954 0.14 0.64
# Prop. Mediated (average) -0.0182 -4.0125 1.65 0.88
# Sample Size Used: 899
# Simulations: 50
# Continuous treatment
jobs$treat_cont <- jobs$treat + rnorm(nrow(jobs)) # (hypothetical) continuous treatment
b.contT <- lm(job_seek ~ treat_cont + econ_hard + sex + age, data=jobs)
c.contT <- lm(depress2 ~ treat_cont + job_seek + econ_hard + sex + age, data=jobs)
contcont.cont <- mediate(b.contT, c.contT, sims=50,
treat="treat_cont", mediator="job_seek",
treat.value = 4, control.value = -2)
summary(contcont.cont)
# Causal Mediation Analysis
# Quasi-Bayesian Confidence Intervals
# Estimate 95% CI Lower 95% CI Upper p-value
# ACME -0.00183 -0.06296 0.07 0.88
# ADE 0.08331 -0.12027 0.25 0.40
# Total Effect 0.08148 -0.11249 0.28 0.36
# Prop. Mediated -0.01808 -2.11776 1.33 0.84
# Sample Size Used: 899
# Simulations: 50
# Categorical treatment
# }
# NOT RUN {
b <- lm(job_seek ~ educ + sex, data=jobs)
c <- lm(depress2 ~ educ + job_seek + sex, data=jobs)
# compare two categories of educ --- gradwk and somcol
model.cat <- mediate(b, c, treat="educ", mediator="job_seek", sims=50,
control.value = "gradwk", treat.value = "somcol")
summary(model.cat)
# Causal Mediation Analysis
# Quasi-Bayesian Confidence Intervals
# Estimate 95% CI Lower 95% CI Upper p-value
# ACME 0.03009 0.00094 0.06 0.04 *
# ADE -0.08757 -0.21384 0.01 0.16
# Total Effect -0.05748 -0.18474 0.06 0.40
# Prop. Mediated -0.25516 -7.25997 4.41 0.44
# ---
# Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# Sample Size Used: 899
# Simulations: 50
