中介分析的实现与多个 R 包横向测评

文摘科学 2024-10-27 12:06 上海

中介分析用于探讨自变量通过中介变量对因变量的间接影响，广泛应用于社会科学等领域。本文对比了五个R包：mediation、lavaan、psych、piecewiseSEM 和 brms。通过每个包的具体操作和结果展示，分析了其各自的优缺点：mediation 包简单易用，lavaan 支持潜变量，psych 操作便捷，piecewiseSEM 提供可视化和结构方程模型支持，brms 适合复杂模型。本文总结了不同R包在中介分析中的适用性，为研究者提供了实用的选择指南。

简介

在某些情况下，我们可能会考虑某个变量对结果或结果的间接影响。例如，童年时期家庭生活条件差可能会降低在校的学习成绩，进而对未来的生活质量（如终身收入）产生负面影响。另一个例子是，我们可能考虑在多个时间点收集的单个变量，这样就存在变量在时间1对时间2的影响，时间2对时间3的影响。基本思想是：


A → B → C

换句话说，A导致B，然后B导致C。通过中介模型，我们假设在我们可能在标准回归设置中具有的协变量→结果路径之间存在一个中介变量，这些模型允许我们研究这种行为。在上述情况下，中介变量是B。通常情况下，我们可能仍然有A对C的直接影响，但与一般模型一样，这应该是有理论依据的。

中介分析在社会科学学科中非常流行，尽管不限于此，通常是在结构方程模型（SEM）的指导下进行的，而SEM本身是更一般的图形模型的特定取向。中介模型的图形模型可能如下所示。

在这种情况下，a和b反映了X通过中介对结果的间接路径影响，而c'是X在间接路径被移除后对结果的直接影响（c是假设没有间接影响时的效果，c - c'等于间接效应）。X的总效应是间接和直接效应的组合。

需要注意的是，根据我在不同学科的咨询经验，很少有人真正需要中介模型。例如，如果你不能以时间或物理方式思考你的模型，使得X必然导致中介，然后必然导致结果，那么你可能不需要中介模型。如果你认为箭头可以朝任意方向，那么你可能不需要这样的模型。另外，如果在描述你的模型时，每个人都认为你在谈论交互作用（即调节作用），你可能也不需要中介模型。最后，如果关键变量（X）和中介变量之间没有强相关性（路径a），并且中介变量和结果之间也没有强相关性（路径b），你可能不需要中介模型。虽然没有什么能阻止你进行中介分析，但如果没有这些先决条件，你的模型几乎肯定会很弱，可能比你原本的模型更令人困惑。

简而言之，当变量之间存在强烈的因果关系时，中介分析效果最佳。即使在这种情况下，也应将这种模型与没有中介的简单模型进行比较。在任何情况下，在R中有一些非常简单的方法来研究这种模型，这是本文的目标，仅仅是演示如何开始。

数据

为了演示不同R包的中介模型，我们将使用mediation包附带的jobs数据集。以下是描述。

求职干预研究（JOBS II）：JOBS II是一个随机实地实验，研究了求职培训干预对失业工人的有效性。该项目旨在不仅提高失业者的再就业率，还提高求职者的心理健康状况。在JOBS II实地实验中，1801名失业工人接受了预筛选问卷，然后被随机分配到处理组和控制组。处理组的参与者参加了求职技能研讨会。在研讨会中，受访者学习了求职技能和应对求职过程中挫折的策略。控制组则收到一本描述求职技巧的小册子。在后续访谈中，两个关键的结果变量是基于霍普金斯症状清单的抑郁症状的连续测量，以及一个表示受访者是否已就业的二元变量。

以下是本次演示中使用的变量描述。还有其他可供您试用的变量。

econ_hard：处理前的经济困难程度，取值为1到5。
sex：性别的指示变量。1 = 女性
age：年龄（单位：岁）。
educ：教育程度的五个类别的因子。
job_seek：衡量求职自我效能感水平的连续量表，取值为1到5。中介变量。
depress2：处理后的抑郁症状测量。结果变量。
treat：参与者是否被随机选中参加JOBS II培训项目的指示变量。1 = 参与。


data(jobs, package = 'mediation')

模型

鉴于这些数据，中介模型和结果模型如下：

中介模型：

job_seek ~ treatment + econ_hard + sex + age

结果模型：

depression ~ treatment + econ_hard + sex + age + job_seek

因此，我们期望求职技能培训对抑郁症状有负面影响（即提高福祉），但至少部分原因是由于对求职的积极影响。

作为一个图形模型，我们可以简洁地将其描述如下：

R包介绍

我们将介绍以下R包，以演示如何在R中进行中介分析：

mediation
lavaan
psych
brms

尽管这些将是重点，但我还会提到一些其他的替代方案，包括Python和Stata。

mediation 包

我们将从mediation包开始，因为它基本上不需要比运行R中标准回归模型更多的编程能力。该包提供了平均因果中介效应（ACME），其定义如下（来自帮助文件）：

平均因果中介效应（ACME）表示当中介变量取其在处理条件下实现的值而不是在控制条件下实现的值时，潜在结果的预期差异，同时处理状态本身保持不变。

请注意，这一定义侧重于根据处理值的条件预期或预测值。这种反事实的概念，或者说在相反设置下观测值会是什么样子，在建模中已经有很长的历史。可以这样理解，如果一个人在处理组中，他们会有一个特定的中介变量值，并且鉴于此，他们会有一个特定的结果预期值。然而，我们也可以假设同一观测值也在控制组中，并评估通过中介对结果的影响。同样，我们可以在保持处理不变的情况下评估潜在结果。考虑给定中介变量值的结果变化不假设模型类型。这就是mediation包能够将不同的模型用于中介变量和结果的原因。例如，中介变量可以是二元的，需要逻辑回归模型，而结果模型可能是生存模型。

在我们的示例中，我们将坚持使用标准的（正态）线性模型。还要注意，尽管我们的处理是二元变量，但这可以推广到连续情况，在这种情况下，我们考虑对“处理”的单位变动的结果。

要使用mediation包，只需运行中介变量和结果的各自模型，然后使用mediate函数获得最终结果。


library(mediation)
model_mediator <- lm(job_seek ~ treat + econ_hard + sex + age, data = jobs)model_outcome  <- lm(depress2 ~ treat + econ_hard + sex + age + job_seek, data = jobs)
# 通过准贝叶斯近似进行估计mediation_result <- mediate(  model_mediator,   model_outcome,   sims = 500,  treat = "treat",  mediator = "job_seek")

结果基于对系数的多元正态分布的模拟。算法总结如下（来自Imai等人，2010）：

拟合观察到的结果和中介变量的模型。
从模型参数的采样分布中模拟模型参数。
重复以下三个步骤：

模拟中介变量的潜在值，
给定中介变量的模拟值，模拟潜在结果，
计算因果中介效应。

计算总结统计量，例如点估计和置信区间。

通过这种方法，我们可以根据模拟的抽样值获得平均差异和相应的分位数。


summary(mediation_result)plot(mediation_result)

输出：


Estimate    95% CI Lower    95% CI Upper    p-value
ACME        -0.016          -0.038          0.009     0.220
ADE         -0.045          -0.127          0.047     0.292
Total Effect-0.061          -0.149          0.027     0.188
Prop. Med.  0.226           -3.222          1.596     0.344

上述结果表明，ACME在统计上不显著异于零，或没有中介作用。平均直接效应（ADE）为负，但同样在统计上不显著。总效应（间接效应+直接效应）也不显著。还提供了所谓的“中介效应占比”，即间接效应与总效应的比率。然而，这不是一个比例，甚至可能为负，因此大多是一个无意义的数字。

优点

使用标准的R模型和语法
对中介变量和结果的多种模型类型
同时提供多种结果
良好的文档和相关的文章可免费获取
可以进行“调节的”中介

局限

依赖于MASS包
对简单的随机效应模型支持有限
在一些模型复杂性方面功能可能有限
不支持潜变量

piecewiseSEM包

piecewiseSEM 包的功能与 mediation 包类似，但它支持更多的模型类型，并能生成附加输出（如标准化结果）和更多选项（如多组分析、相关残差），还可以对模型进行可视化。


library(piecewiseSEM)
model_mediator <- lm(job_seek ~ treat + econ_hard + sex + age, data = jobs)model_outcome  <- lm(depress2 ~ treat + econ_hard + sex + age + job_seek, data = jobs)
mediation_result <-  psem(model_mediator, model_outcome, data = jobs)
summary(mediation_result)

运行结果如下：


Structural Equation Model of mediation_result 

Call:
  job_seek ~ treat + econ_hard + sex + age
  depress2 ~ treat + econ_hard + sex + age + job_seek

    AIC      BIC
 26.000   88.417

---  
Tests of directed separation:

 No independence claims present. Tests of directed separation not possible.

Global goodness-of-fit:

  Fisher's C = 0 with P-value = 1 and on 0 degrees of freedom

---  
Coefficients:
  
  Response Predictor Estimate Std.Error  DF Crit.Value P.Value Std.Estimate    
  job_seek     treat   0.0656    0.0515 894     1.2748  0.2027       0.0425    
  job_seek econ_hard   0.0532    0.0246 894     2.1612  0.0309       0.0720   *
  job_seek       sex  -0.0076    0.0487 894    -0.1567  0.8755      -0.0052    
  job_seek       age   0.0046    0.0023 894     1.9779  0.0482       0.0658   *
  depress2     treat  -0.0403    0.0435 893    -0.9255  0.3550      -0.0291    
  depress2 econ_hard   0.1485    0.0208 893     7.1323  0.0000       0.2248 ***
  depress2       sex   0.1068    0.0411 893     2.5957  0.0096       0.0818  **
  depress2       age   0.0006    0.0020 893     0.3306  0.7410       0.0104    
  depress2  job_seek  -0.2400    0.0282 893    -8.4960  0.0000      -0.2682 ***

  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05

可以使用 piecewiseSEM 的绘图功能快速生成模型的可视化图：


plot(mediation_result)

该图由 Diagrammer 和 graphviz 支持，可以用于进一步的定制化处理。需要注意的是，目前尚无自动计算间接效应的功能，需要手动进行计算。此外，原有的 semEff 包以前支持 psem 对象的间接效应计算，但目前不可用。

优点

支持标准 R 模型和语法
支持多种类型的模型（如线性和广义线性模型）
可以生成结构方程模型样式的结果（如模型拟合、标准化系数、AIC）
可以快速生成模型的可视化图

局限性

不支持自动计算间接效应
不支持潜变量

lavaan 包

在中介模型和结果模型都是具有正态分布目标变量的标准线性模型的特定情况下，间接效应等于前面图中a和b路径的乘积。直接效应是c'路径。独立直接效应（我们可以称之为c）与在中介模型中估计的直接效应c'的比较是c - c' = a*b。之前提到的内容现在可能更清楚了，如果a或b接近零，那么间接效应只能接近零，因此谨慎地事先研究这种关系是明智的。

这种路径乘积（或系数差异）的方法是我们在lavaan包中将采用的方法，事实上，截止本文撰写时，这是我们唯一的方法。lavaan专门用于结构方程模型，例如因子分析、增长模型和我们在此进行的中介模型，强烈建议用于此类模型。虽然它在评估中介方面仅限于标准线性模型案例，但它是我们工具中唯一能够轻松结合潜变量的工具。例如，我们可以将抑郁结果作为潜在变量，潜在于各个问卷项目之下。此外，我们还可以结合多个中介和多个结果。

为了保持我们一直讨论的内容，我将在lavaan中按照之前描述的方式标记a、b和c'路径。否则，lavaan非常容易使用，对于观测变量，使用标准的R公式表示模型。除此之外，我们使用:=操作符定义我们要计算的感兴趣的效应。我们将整个模型作为一个简单的字符字符串指定，然后使用sem函数进行分析。


library(lavaan)
sem_model = '  job_seek ~ a*treat + econ_hard + sex + age  depress2 ~ c*treat + econ_hard + sex + age + b*job_seek   # 直接效应  direct := c   # 间接效应  indirect := a*b   # 总效应  total := c + (a*b)'
model_sem = sem(sem_model, data=jobs, se='boot', bootstrap=500)summary(model_sem, rsq=T)  # 与mediation包中的ACME比较

输出：


lavaan 0.6-6 ended normally after 25 iterations

  Estimator                                         ML
  Optimization method                           NLMINB
  Number of free parameters                         11
                                                      
  Number of observations                           899
                                                      
Model Test User Model:
                                                      
  Test statistic                                 0.000
  Degrees of freedom                                 0

Parameter Estimates:

  Standard errors                            Bootstrap
  Number of requested bootstrap draws              500
  Number of successful bootstrap draws             500

Regressions:
                   Estimate  Std.Err  z-value  P(>|z|)
  job_seek ~                                          
    treat      (a)    0.066    0.049    1.332    0.183
    econ_hard         0.053    0.024    2.242    0.025
    sex              -0.008    0.047   -0.163    0.870
    age               0.005    0.002    1.934    0.053
  depress2 ~                                          
    treat      (c)   -0.040    0.044   -0.905    0.365
    econ_hard         0.149    0.022    6.908    0.000
    sex               0.107    0.038    2.831    0.005
    age               0.001    0.002    0.332    0.740
    job_seek   (b)   -0.240    0.030   -8.079    0.000

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)
   .job_seek          0.524    0.030   17.610    0.000
   .depress2          0.373    0.022   17.178    0.000

R-Square:
                   Estimate
    job_seek          0.011
    depress2          0.120

Defined Parameters:
                   Estimate  Std.Err  z-value  P(>|z|)
    direct           -0.040    0.045   -0.904    0.366
    indirect         -0.016    0.012   -1.324    0.185
    total            -0.056    0.046   -1.224    0.221

我们看到与之前相同的输出，可以将我们的间接参数与之前的ACME进行比较，直接效应与ADE比较，总效应与之前的总效应比较。数值基本相同。

还要注意的是，输出显示了两个模型的R²值。就job_seek而言，我们可以看到，我们没有发现太多的中介作用的原因是涉及的协变量一开始就没有解释中介变量的任何变化。在这种情况下，初步调查本可以节省我们的麻烦。

优点

可以处理多个中介
可以处理多个“处理”变量
可以处理多个结果
可以使用潜变量
支持一些多层次模型
可以进行调节的中介和被中介的调节（尽管不适用于潜变量）

局限

需要额外的编码来估计间接效应
仅支持单个随机效应
虽然模型可以结合中介/结果的二元或有序变量，但在这些设置中没有直接的方法来计算类似于mediation包中的间接效应

psych 包

psych包利用了在标准线性模型情况下，可以仅根据协方差矩阵通过适当的回归模型获得结果的事实。它与lavaan非常相似，尽管使用的是普通最小二乘法，而不是最大似然法。这里的好处是语法允许你仅关注感兴趣的效应，或者包括所有内容，如果你对经济困难、年龄和性别的间接效应感兴趣，这很不错。

对于这个演示，我们将使用使用-而不是+的简化版本。这仅意味着它们被包括在模型中，但没有显示关于它们的结果。中介变量用()标识。另一个好处是可以快速绘制结果，显示未调整和调整后的直接效应之间的差异，以及适当的自举区间。



library(psych)
mediation_psych = mediate(  depress2 ~ treat + (job_seek) - econ_hard - sex - age,   data = jobs,  n.iter = 500)
mediation_psych

输出：


Mediation/Moderation Analysis 
Call: mediate(y = depress2 ~ treat + (job_seek) - econ_hard - sex - 
    age, data = jobs, n.iter = 500)

The DV (Y) was  depress2* . The IV (X) was  treat* . The mediating variable(s) =  job_seek* . Variable(s)  partialled out were econ_hard sex age

Total effect(c) of  treat*  on  depress2*  =  -0.06   S.E. =  0.05  t  =  -1.24  df=  895   with p =  0.21
Direct effect (c') of  treat*  on  depress2*  removing  job_seek*  =  -0.04   S.E. =  0.15  t  =  14.91  df=  893   with p =  4.6e-45
Indirect effect (ab) of  treat*  on  depress2*  through  job_seek*   =  -0.02 
Mean bootstrapped indirect effect =  -0.02  with standard error =  0.01  Lower CI =  -0.04    Upper CI =  0.01
R = 1.07 R2 = 1.15   F = -3510.4 on 2 and 893 DF   p-value:  1 

 To see the longer output, specify short = FALSE in the print statement or ask for the summary

summary(mediation_psych)

输出：


Call: mediate(y = depress2 ~ treat + (job_seek) - econ_hard - sex - 
    age, data = jobs, n.iter = 500)

Direct effect estimates (traditional regression)    (c') 
          depress2*   se     t  df     Prob
Intercept      2.21 0.15 14.91 893 4.60e-45
treat         -0.04 0.04 -0.93 893 3.55e-01
job_seek      -0.24 0.03 -8.50 893 8.14e-17

R = 1.07 R2 = 1.15   F = -3510.4 on 2 and 893 DF   p-value:  1 

 Total effect estimates (c) 
      depress2*   se     t  df  Prob
treat     -0.06 0.05 -1.24 895 0.215

 'a'  effect estimates 
          job_seek   se     t  df      Prob
Intercept     3.67 0.13 29.33 894 5.65e-133
treat         0.07 0.05  1.27 894  2.03e-01

 'b'  effect estimates 
         depress2*   se    t  df     Prob
job_seek     -0.24 0.03 -8.5 894 7.83e-17

 'ab'  effect estimates (through mediators)
      depress2*  boot   sd lower upper
treat     -0.02 -0.02 0.01 -0.04  0.01

相同的结果，不同的包装，但可能是迄今为止最简单的方法，因为只需要一个函数调用。psych包还处理多个中介和结果，这是一个额外的好处。

优点

最简单的语法，基本上是一个行的模型
快速绘制结果
可以处理多个中介、“处理”和结果
可以进行“调节的”中介

局限

仅限于标准线性模型（lm）
依赖于MASS包

brms 包

在下一个演示中，我们将介绍我认为最强大的包brms。它的名称代表使用Stan进行贝叶斯回归建模，Stan是一个用于贝叶斯分析的强大概率编程语言。我不会详细讨论贝叶斯分析，但可以查看我的相关文档。

我们通常像以前一样，指定中介模型和结果模型。brms在中介分析方面没有做任何特别的事情，但它的hypothesis函数可以让我们测试路径乘积的方法。此外，sjstats包基本上会以与mediation包相同的方式为我们提供结果，实际上，mediation包基本上是使用频率方法的贝叶斯解决方案。如果我们确实对结果和中介有不同的分布，我们将相对容易地获得这些平均预测值及其差异，因为贝叶斯方法始终考虑后验预测分布。无论如何，以下是代码。


library(brms)
model_mediator <- bf(job_seek ~ treat + econ_hard + sex + age)model_outcome  <- bf(depress2 ~ treat + job_seek + econ_hard + sex + age)
med_result = brm(  model_mediator + model_outcome + set_rescor(FALSE),   data = jobs)save(med_result, file = 'data/mediation_brms.RData')
load('data/mediation_brms.RData')summary(med_result)

输出：


Family: MV(gaussian, gaussian) 
 Links: mu = identity; sigma = identity
        mu = identity; sigma = identity 
Formula: job_seek ~ treat + econ_hard + sex + age 
         depress2 ~ treat + job_seek + econ_hard + sex + age 
   Data: jobs (Number of observations: 899) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Population-Level Effects: 
                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
jobseek_Intercept      3.67      0.12     3.43     3.91 1.00     6699     3749
depress2_Intercept     2.21      0.15     1.92     2.50 1.00     6174     3091
jobseek_treat          0.07      0.05    -0.03     0.17 1.00     6322     2709
jobseek_econ_hard      0.05      0.02     0.00     0.10 1.00     6266     2656
jobseek_sex           -0.01      0.05    -0.10     0.09 1.00     5741     2655
jobseek_age            0.00      0.00     0.00     0.01 1.00     6539     2846
depress2_treat        -0.04      0.04    -0.12     0.04 1.00     5458     3102
depress2_job_seek     -0.24      0.03    -0.30    -0.18 1.00     5950     2938
depress2_econ_hard     0.15      0.02     0.11     0.19 1.00     7543     3102
depress2_sex           0.11      0.04     0.03     0.19 1.00     5599     2699
depress2_age           0.00      0.00    -0.00     0.00 1.00     4555     2887

Family Specific Parameters: 
               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma_jobseek      0.73      0.02     0.69     0.76 1.00     6639     3276
sigma_depress2     0.61      0.01     0.59     0.64 1.00     6145     2987

Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

使用brms，我们可以按以下方式计算间接效应：



# 使用hypothesis函数hypothesis(med_result, 'jobseek_treat*depress2_job_seek = 0')
# sjstats包提供了与mediation包类似的输出# print(sjstats::mediation(med_result), digits=4)sjstats::mediation(med_result) %>% kable_df()

输出：

effect

value

hdi.low

hdi.high


direct	-0.039	-0.112	0.031
indirect	-0.015	-0.036	0.005
mediator	-0.240	-0.286	-0.193
total	-0.055	-0.133	0.017
proportion mediated	0.277	-0.813	1.366

在输出中，任何带有jobseek_*的都是中介模型的结果，而depress2_*是结果模型的结果。此时，我们有相同的结果，但使用贝叶斯方法，我们可以看到更多有趣的内容。例如，我们可以看到我们实际上并没有很好地捕捉到抑郁结果的偏度。我们的预测值与观察值并不完全匹配。对于中介变量，我们要好一些，但基于我们模型的预测可能仍然有点高。


pp_check(med_result, resp = 'depress2') + ggtitle('Depression Outcome')pp_check(med_result, resp = 'jobseek') + ggtitle('Mediator')

淡色的线条（yrep）是后验预测抽样值。

优点

语法简洁
功能极其强大——模型主要受想象力限制
基本上完成了mediation包近似的工作
贝叶斯推断的所有优点：诊断、后验预测检查、模型比较等

局限

估计速度较慢
超过标准线性模型时需要“手动”计算，但从贝叶斯的角度来看，这已经是常见的方法
需要对贝叶斯方法有一定的熟悉程度
更复杂

替代方法

我们一直将模型描述为节点的网络，节点之间由弧/边/路径连接。我们的讨论围绕着所谓的有向无环图（DAG），其中箭头只能朝一个方向，没有反馈循环。任何结果变量的结果是之前箭头的函数，并且在条件上与其他变量独立。一些理论模型可能会放宽这一点，其他模型可能根本没有箭头，即无向的，这样我们只对连接感兴趣（例如，在一些社交网络中）。

bnlearn

bnlearn包允许研究有向、部分有向和无向图。在DAG方面，我们可以用它来基本上复制我们一直在讨论的中介模型。不过，值得注意的是，该包将有效地测试路径是否包含，而不是假设它们，但我们仍然可以根据需要施加理论约束。不仅如此，我们还可以以贝叶斯网络和Pearl的因果图理论为基础，以原则性的方式搜索感兴趣的路径，我们还将有工具通过交叉验证进一步避免过拟合。

对于初始模型，我们将确保在处理-中介、处理-结果和中介-结果之间存在路径（白名单）。我们将禁止非理性的路径，例如指向处理的箭头（处理是随机分配的）、性别、经济困难和年龄（黑名单）。否则，我们将看看数据表明什么。


whitelist = data.frame(  from = c('treat', 'treat', 'job_seek'),  to   = c('job_seek', 'depress2', 'depress2'))
blacklist = expand.grid(  from = colnames(mediation_result$model.y$model),  to   = c('treat', 'sex', 'age', 'econ_hard'))


# 为了简化输出，我们将处理和性别作为数值型（后面解释）library(dplyr)
jobs_trim = jobs %>%   select(depress2, treat, econ_hard, sex, age, job_seek) %>%   mutate(    treat = as.numeric(jobs$treat),    sex = as.numeric(jobs$sex)    )
library(bnlearn)
model = gs(jobs_trim, whitelist = whitelist, blacklist = blacklist)
plot(model)

我们在图中看到事情发生了一些变化。例如，年龄现在仅与求职自我效能感相关，性别仅对抑郁有影响。

如果我们将路径限制为仅与之前的示例相同，我们将得到相同的结果。


whitelist = data.frame(  from = c('treat', 'age', 'sex', 'econ_hard', 'treat', 'job_seek', 'age', 'sex', 'econ_hard'),  to   = c('job_seek', 'job_seek','job_seek','job_seek', 'depress2', 'depress2', 'depress2', 'depress2', 'depress2'))
blacklist = expand.grid(  from = colnames(mediation_result$model.y$model),  to   = c('treat', 'sex', 'age', 'econ_hard')) 
model = gs(jobs_trim, whitelist = whitelist, blacklist = blacklist)plot(model)
parameters = bn.fit(model, jobs_trim)
parameters$depress2$coefficients
 
(Intercept)         treat     econ_hard           sex           age      job_seek 
2.2076414333 -0.0402647000  0.1485433818  0.1068048699  0.0006488642 -0.2399549527 

parameters$job_seek$coefficients

(Intercept)        treat    econ_hard          sex          age 
3.670584908  0.065615003  0.053162413 -0.007637336  0.004586492

主要需要注意的是，估计的参数等于我们从之前的包中得到的结果。它本质上等同于使用lavaan的默认最大似然估计。

如果我们将处理和性别作为因子，bnlearn将根据因子取值的不同生成不同的条件模型。换句话说，对于treatment == 'treatment'和treatment == control，将有一个单独的模型。在我们的例子中，这将等同于允许一切与处理交互，例如lm(job_seek ~ treat * (econ_hard + sex + age))，抑郁模型也是如此。这可能会扩展到任何二元变量（例如，包括性别）。如果中介变量是二元变量，这可能是我们想要做的。

使用 Python

我在这里展示下statsmodels的实现。它遵循Imai的方法，因此可以被视为mediation包的Python版本。输出与将处理作为因子变量时得到的基本相同，在这种情况下，您会获得每个处理类别的单独结果。对于我们的演示，这是不必要的，因此您可以将“平均”结果与之前的mediation包结果进行比较。



import statsmodels.api as smfrom statsmodels.stats.mediation import Mediationimport numpy as npimport pandas as pd
outcome_model = sm.OLS.from_formula("depress2 ~ treat + econ_hard + sex + age + job_seek",                                    data = jobs)
mediator_model = sm.OLS.from_formula("job_seek ~ treat + econ_hard + sex + age",                                     data = jobs)
med = Mediation(outcome_model, mediator_model, "treat", "job_seek")
med_result = med.fit(n_rep = 500)
print(np.round(med_result.summary(), decimals = 3))

输出：


                          Estimate  Lower CI bound  Upper CI bound  P-value
ACME (control)              -0.016          -0.048           0.014    0.332
ACME (treated)              -0.016          -0.048           0.014    0.332
ADE (control)               -0.043          -0.130           0.044    0.308
ADE (treated)               -0.043          -0.130           0.044    0.308
Total effect                -0.059          -0.144           0.029    0.208
Prop. mediated (control)     0.241          -1.710           2.254    0.364
Prop. mediated (treated)     0.241          -1.710           2.254    0.364
ACME (average)              -0.016          -0.048           0.014    0.332
ADE (average)               -0.043          -0.130           0.044    0.308
Prop. mediated (average)     0.241          -1.710           2.254    0.364

使用 Stata

最后，我提供了在Stata中使用其sem命令的选项。Stata使得在此示例中获取间接效应变得容易，但它对每个协变量都这样做，因此输出有点冗长。对于使用Stata的人来说，他们不需要单独的SEM包来获得这些结果。

use "data\jobs.dta"
sem (job_seek <- treat econ_hard sex age) (depress2 <- treat econ_hard sex age job_seek), cformat(%9.3f) pformat(%5.2f)
estat teffects, compact cformat(%9.3f) pformat(%5.2f)

输出：

（由于篇幅限制，此处不展示Stata的完整输出）

总结

具有间接效应的模型需要谨慎的理论考虑才能用于数据分析。然而，如果模型适合您的数据情况，从R的各种包中获得结果是相当容易的。此外，进行具有间接效应的分析并不需要使用结构方程模型包，事实上，使用标准的R语法可以达到很大的效果。对于严格的观测变量（即没有潜变量），不需要使用SEM工具，甚至不推荐。

R包功能比较总结

以下表格可能有助于根据您的理论考虑来决定使用哪个包来满足您的需求。

功能	mediation包	lavaan包	piecewiseSEM包	psych包	brms包
自动计算间接效应	•		•	•	*
多个处理变量	•	•	•	•	•
多个中介变量	•	•	•	•	•
多个结果变量		•	•	•	•
超越标准线性模型	•	•	•		•
随机效应	•	•	•		•
处理缺失值		•		*	•
支持潜变量		•		*

注释：

* 表示大致支持，带有一些注意事项。
• 表示支持。

http://mp.weixin.qq.com/s?__biz=MzkwNTM5MzE4MA==&mid=2247488335&idx=1&sn=f096ed72a2a33c97eefcb35af2f13e60

真实世界数据

介绍真实世界数据，真实世界研究和真实世界证据

中介分析的实现与多个 R 包横向测评

目录

简介

数据

模型

R包介绍

mediation 包

piecewiseSEM包

优点

局限性

lavaan 包

psych 包

brms 包

更多复杂的场景

交互作用

广义线性模型

缺失数据

替代方法

bnlearn

使用 Python

使用 Stata

总结

R包功能比较总结