温忠麟、叶宝娟，2014，《中介效应分析:方法和模型发展》，《心理科学进展》第5期。 连享会：https://zhuanlan.zhihu.com/p/99435552
中介效应地基本方程：
\[ \begin{gathered} Y=c X+e_{1} & (1) \\ M=a X+e_{2} & (2) \\ Y=c^{\prime} X+b M+e_{3} & (3) \end{gathered} \]
其中\(X\)是自变量，\(M\)是中介变量，\(Y\)是因变量；\(c\)是总效应，\(a*b\)是间接效应，\(c^{\prime}是直接效应\)。
数据基本描述：这是一组有关大型百货公司销售人员的数据，我们用来讨论经理的激励与员工工作表现之间的关系，基本假设是：经理的激励 (perceived support from managers) 可能通过影响员工的工作满意度 (job satisfaction) 而影响员工的工作表现 (job performance)。
dat <- haven::read_dta("http://www.stata-press.com/data/r15/gsem_multmed.dta")
dim(dat)## [1] 1500 4head(dat)## # A tibble: 6 × 4
## branch support satis perform
## <dbl> <dbl> <dbl> <dbl>
## 1 1 0.100 0.200 4.44
## 2 2 0.600 -0.5 3.64
## 3 3 0.200 -0.100 5.40
## 4 4 0.300 -0.100 4.34
## 5 5 0.900 -0.300 5.37
## 6 6 -0.200 0.700 7.07数据分布可以看到，我们的样本量为 1500， 其中所需使用的变量为：
# 分析x和y之间的关系
library(estimatr)
library(tidyverse)
lm_robust(perform ~ support, data = dat) %>% summary()##
## Call:
## lm_robust(formula = perform ~ support, data = dat)
##
## Standard error type: HC2
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
## (Intercept) 4.9984 0.02049 243.96 0.000e+00 4.9582 5.0385 1498
## support 0.8218 0.04063 20.23 1.273e-80 0.7421 0.9014 1498
##
## Multiple R-squared: 0.2157 , Adjusted R-squared: 0.2152
## F-statistic: 409.1 on 1 and 1498 DF, p-value: < 2.2e-16结果显示员工的工作表现与经理的激励显著相关，回归系数 c=0.822，可以进行下一步检验。
# 分析 x 和 m 之间的关系
lm_robust(satis ~ support, data = dat) %>% summary()##
## Call:
## lm_robust(formula = satis ~ support, data = dat)
##
## Standard error type: HC2
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
## (Intercept) 0.01926 0.01544 1.247 2.124e-01 -0.01103 0.04955 1498
## support 0.22889 0.02953 7.752 1.653e-14 0.17098 0.28681 1498
##
## Multiple R-squared: 0.03618 , Adjusted R-squared: 0.03553
## F-statistic: 60.1 on 1 and 1498 DF, p-value: 1.653e-14回归结果显示， 经理的激励显著增加员工的工作满意度，系数 a=0.229。
lm_robust(perform ~ satis + support, data = dat) %>% summary()##
## Call:
## lm_robust(formula = perform ~ satis + support, data = dat)
##
## Standard error type: HC2
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
## (Intercept) 4.9811 0.01509 330.12 0.000e+00 4.9515 5.0107 1497
## satis 0.8984 0.02358 38.09 1.436e-222 0.8522 0.9447 1497
## support 0.6161 0.03058 20.15 4.459e-80 0.5561 0.6761 1497
##
## Multiple R-squared: 0.5756 , Adjusted R-squared: 0.575
## F-statistic: 1070 on 2 and 1497 DF, p-value: < 2.2e-16在加入工作满意度后，员工的表现和经理激励之间的显著关系没有发生变化，但是系数由第一步的 \(c=0.822\) 减小到 \(c’=0.616\)，员工的工作满意度和员工的表现之间显著相关，系数\(b=0.898\)，说明工作满意度在经理激励和员工表现之间起到了部分中介的作用。
完整结果解读:
另附bootstrap计算中介效应的方法：
get_abc <- function(data){
fit_c <- coef(lm(perform ~ support, data = data))[[2]] # 总效应
fit_a <- coef(lm(satis ~ support, data = data))[[2]] # x对m的效应
fit_cprime_b <- lm(perform ~ satis + support, data = data)
fit_b <- coef(fit_cprime_b)[[2]] # m对y的效应
fit_cprime <- coef(fit_cprime_b)[[3]] # 直接效应
return(list(c = fit_c, a = fit_a, b = fit_b, cprime = fit_cprime))
}
# 1000次
set.seed(123)
res <- map_dfr(1:500,
~get_abc(data = (dat %>% sample_n(size = nrow(dat), replace = T))))
# 输出结果
res %>%
mutate(a_times_b = a * b) %>%
select(c, cprime, a_times_b) %>%
map_dfr(~list(Mean = mean(.),
upper = quantile(., .025),
lower = quantile(., .975)),
.id = "effects")## # A tibble: 3 × 4
## effects Mean upper lower
## <chr> <dbl> <dbl> <dbl>
## 1 c 0.819 0.741 0.890
## 2 cprime 0.613 0.557 0.678
## 3 a_times_b 0.205 0.155 0.257其中，第一行是总效应，第二行是直接效应，第三行是间接效应。
当然，我们也可以直接将这些变量作为显变量，纳入结构方程，利用sem package直接计算，参见：https://rpubs.com/Plumber/862703
对于模型：
\[ Y=\alpha_{F}+\beta_{F} X+\gamma_{F} Z+\delta_{F} C+\epsilon \tag{1} \]
其中，\(Y\)是因变量，\(X\)是自变量，\(Z\)是中介变量，\(C\)是控制变量。因此，\(\beta_F\)就是直接效应。
总效应可以在reduced model中进行估计：
\[ Y=\alpha_{R}+\beta_{R} X+\delta_{R} C+\varepsilon \tag{2} \]
那么相减就可以算出间接效应：
\[ \beta_{I}=\beta_{R}-\beta_{F} \tag{3} \]
在一般线性回归中，都是用的上面这种计算两个嵌套模型(1)和(2)，然后对系数进行相减的方法。
但是在广义线性模型中，对嵌套模型的系数进行比较和分析并没有那么简单。因为即使\(Z\)和\(X\)完全不相关，引入\(Z\)也必然改变\(X\)的系数。
广义线性模型，例如logit模型，可以视为对某一潜变量进行建模\(Y^*\)建模的模型：
\[ \begin{aligned} Y^{*} &=\alpha_{F}+\beta_{F} X+\gamma_{F} Z+\delta_{F} C+\epsilon & (1) \\ Y^{*} &=\alpha_{R}+\beta_{R} X+\delta_{R} C+\varepsilon & (2) \end{aligned} \]
一般来说，我们得到的估计量是:
\[ b_{F}=\frac{\beta_{F}}{\sigma_{F}} \quad \text { and } \quad b_{R}=\frac{\beta_{R}}{\sigma_{R}} \]
其中\(\sigma_F\)和\(\sigma_R\)被称为尺度参数，如果直接相减，得到的是\(b_{R}-b_{F}=\frac{\beta_{R}}{\sigma_{R}}-\frac{\beta_{F}}{\sigma_{F}}\)
KHB分解遵循以下步骤:
第一步：估计直接效应：
\[ Y^{*}=\alpha_{F}+\beta_{F} X+\gamma_{F} Z+\delta_{F} C+\epsilon \]
第二步：将\(Z\)中与\(X\)线性相关的部分剥离：
\[ R=Z-(a+b X) \]
第三步：用\(R\)替代\(Z\)，来估计总效应：
\[ Y^{*}=\widetilde{\alpha}_{R}+\widetilde{\beta}_{R} X+\widetilde{\gamma}_{R} R+\widetilde{\delta}_{R} C+\epsilon \]
第二步和第三步的作用，是保证两个模型的尺度参数相同，则有：
\[ \widetilde{b}_{R}-b_{F}=\frac{\widetilde{\beta}_{R}}{\widetilde{\sigma}_{R}}-\frac{\beta_{F}}{\sigma_{F}}=\frac{\beta_{R}-\beta_{F}}{\sigma_{F}} \] 也就可以计算直接效应和间接效应的占比了：
# data
dat <- haven::read_dta("http://fmwww.bc.edu/repec/bocode/d/dlsy_khb.dta")
dim(dat)## [1] 1896 8head(dat)## # A tibble: 6 × 8
## edu upsec univ fgroup fses abil intact boy
## <dbl+lbl> <dbl+lb> <dbl+l> <dbl+l> <dbl> <dbl> <dbl+l> <dbl+l>
## 1 1 [Compulsory schooling] 0 [no] 0 [no] 1 [Low… -2.01 -1.21 1 [yes] 1 [yes]
## 2 1 [Compulsory schooling] 0 [no] 0 [no] 1 [Low… -2.09 0.972 1 [yes] 0 [no]
## 3 1 [Compulsory schooling] 0 [no] 0 [no] 1 [Low… -1.87 0.302 1 [yes] 1 [yes]
## 4 1 [Compulsory schooling] 0 [no] 0 [no] 1 [Low… -2.09 0.604 1 [yes] 0 [no]
## 5 3 [University] 1 [yes] 1 [yes] 1 [Low… -2.09 1.91 1 [yes] 1 [yes]
## 6 1 [Compulsory schooling] 0 [no] 0 [no] 1 [Low… -2.09 0.509 1 [yes] 0 [no]其中：
univ: university completion, 作为二分因变量fses: parental social status，作为自变量abil: academic ability， 作为中介变量# full model(直接效应)
full_fit <- glm(univ ~ fses + abil + boy + intact, data = dat, family = "binomial") %>%
broom::tidy() %>%
filter(term == 'fses')
# reduced model（间接效应）
R <- lm(abil ~ fses, data = dat)
reduced_fit <- glm(univ ~ fses + R$residuals + intact + boy, data = dat, family = binomial("logit")) %>%
broom::tidy() %>%
filter(term == 'fses')
# 第一行是总效应，第二行是直接效应
bind_rows(reduced_fit, full_fit)## # A tibble: 2 × 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 fses 0.548 0.0780 7.02 2.16e-12
## 2 fses 0.382 0.0778 4.91 9.29e- 7# 减一下就是间接效应
reduced_fit[2] - full_fit[2]## estimate
## 1 0.1661177在R中没有进行KHB分解的现成软件包。从可靠性考虑，还是推荐STATA中的khb包，同样是上面的例子
. use http://fmwww.bc.edu/repec/bocode/d/dlsy_khb.dta
. khb logit univ fses || abil, concomitant(intact boy)
Decomposition using the KHB-Method
Model-Type: logit Number of obs = 1896
Variables of Interest: fses Pseudo R2 = 0.19
Z-variable(s): abil
Concomitant: intact boy
------------------------------------------------------------------------------
univ | Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
fses |
Reduced | .5459815 .0779806 7.00 0.000 .3931424 .6988206
Full | .3817324 .0778061 4.91 0.000 .2292353 .5342295
Diff | .1642491 .0293249 5.60 0.000 .1067734 .2217247
------------------------------------------------------------------------------
Reduced就是总效应，Full就是直接效应，Diff就是间接效应。