Free Press

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)# Estimation via quasi-Bayesian approximation?mediatemediation_result <- mediate( model_mediator, model_outcome, sims = 500, treat = "treat", mediator = "job_seek")detach(package:mediation)detach(package:MASS)

summary(mediation_result)plot(mediation_result)

-0.047

0.047

library(lavaan)sem_model = ' job_seek ~ a*treat + econ_hard + sex + age depress2 ~ c*treat + econ_hard + sex + age + b*job_seek # direct effect direct := c # indirect effect indirect := a*b # total effect total := c + (a*b)'model_sem = sem(sem_model, data=jobs, se='boot', bootstrap=500)summary(model_sem, rsq=T) # compare with ACME in mediation

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 0Parameter Estimates: Standard errors Bootstrap Number of requested bootstrap draws 500 Number of successful bootstrap draws 500Regressions: 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.000Variances: 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.000R-Square: Estimate job_seek 0.011 depress2 0.120Defined 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

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---Individual R-squared: Response method R.squared job_seek none 0.01 depress2 none 0.12

plot(mediation_result)

残念ながら、現在のところ間接効果を自動的に計算する方法はありませんので、手作業で結果をブートスト ラップしなければならないでしょう。

Pros

• Standard R models and syntax
• Mediator and outcomeの両方で複数のタイプのモデル
• SEM-style results (e.g…………….)……….。 fit, standardized coefficients, AIC）
• Quick plot of results
• Can handle multiple mediator, ‘treatments’,

Limitations

• Do not automatically calculate indirect effects
• No latent variable capabilities

psych

The psych package takes advantage of fact that one can get results via appropriate regression models based on covariance matrices alone.標準線形モデルケースにおいて、分散行列にのみ基づいて、結果が得られるという事実。 最尤法ではなく、普通の最小二乗法を用いていますが、lavaanと非常によく似ています。 また、経済的困難、年齢、性別の間接効果に興味がある場合は、この構文が便利です。

このデモでは、非処置効果に+の代わりに-を使用したクリーンアップ版を使用します。 これは、モデルには含まれているが、それに関する結果は表示されないことを意味します。 Mediatorは()で識別されます。

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 ageTotal effect(c) of treat* on depress2* = -0.06 S.E. = 0.05 t = -1.24 df= 895 with p = 0.21Direct effect (c') of treat* on depress2* removing job_seek* = -0.04 S.E. = 0.15 t = 14.91 df= 893 with p = 4.6e-45Indirect 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.01R = 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 ProbIntercept 2.21 0.15 14.91 893 4.60e-45treat -0.04 0.04 -0.93 893 3.55e-01job_seek -0.24 0.03 -8.50 893 8.14e-17R = 1.07 R2 = 1.15 F = -3510.4 on 2 and 893 DF p-value: 1 Total effect estimates (c) depress2* se t df Probtreat -0.06 0.05 -1.24 895 0.215 'a' effect estimates job_seek se t df ProbIntercept 3.67 0.13 29.33 894 5.65e-133treat 0.07 0.05 1.27 894 2.03e-01 'b' effect estimates depress2* se t df Probjob_seek -0.24 0.03 -8.5 894 7.83e-17 'ab' effect estimates (through mediators) depress2* boot sd lower uppertreat -0.02 -0.02 0.01 -0.04 0.01

結果は同じでパッケージは異なりますが、関数の呼び出しが1回で済むので、おそらく最も簡単な方法です。 psychパッケージは、ボーナスとして複数のmediatorと結果を扱うことができます。

Pros

• Easiest syntax, basically a one line model
• Quick plot of results
• Can handle multiple mediators, ‘treatments’.All Rights Reserved, そして結果
• 「緩和された」調停が可能

Limitations

• Limited to standard linear model (lm)
• Use of MASS

brms

次のデモでは、最もパワフルだと私が感じた brms というパッケージについてお話したいと思います。 この名前は Bayesian Regression Modeling with Stan の略で、Stan はベイズ解析のための強力な確率的プログラミング言語です。 brms は媒介分析のために特別なことはしませんが、仮説関数によりパスの積のアプローチをテストすることができます。 さらに、sjstatsパッケージは基本的にmediationパッケージと同じように結果を出してくれますし、mediationパッケージは基本的に頻度論的手法によるベイズ的解決を試みています。 もし結果とmediatorに異なる分布があったとしても、ベイズ的アプローチは常に事後予測分布について考えるので、これらの平均予測値やその差を得るのは比較的簡単です。

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 = 4000Population-Level Effects: Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESSjobseek_Intercept 3.67 0.12 3.43 3.91 1.00 6699 3749depress2_Intercept 2.21 0.15 1.92 2.50 1.00 6174 3091jobseek_treat 0.07 0.05 -0.03 0.17 1.00 6322 2709jobseek_econ_hard 0.05 0.02 0.00 0.10 1.00 6266 2656jobseek_sex -0.01 0.05 -0.10 0.09 1.00 5741 2655jobseek_age 0.00 0.00 0.00 0.01 1.00 6539 2846depress2_treat -0.04 0.04 -0.12 0.04 1.00 5458 3102depress2_job_seek -0.24 0.03 -0.30 -0.18 1.00 5950 2938depress2_econ_hard 0.15 0.02 0.11 0.19 1.00 7543 3102depress2_sex 0.11 0.04 0.03 0.19 1.00 5599 2699depress2_age 0.00 0.00 -0.00 0.00 1.00 4555 2887Family Specific Parameters: Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESSsigma_jobseek 0.73 0.02 0.69 0.76 1.00 6639 3276sigma_depress2 0.61 0.01 0.59 0.64 1.00 6145 2987Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESSand Tail_ESS are effective sample size measures, and Rhat is the potentialscale reduction factor on split chains (at convergence, Rhat = 1).

# using brms we can calculate the indirect effect as follows# hypothesis(med_result, 'jobseek_treat*depress2_job_seek = 0')# sjstats provides similar printing as the mediation package# print(sjstats::mediation(med_result), digits=4)sjstats::mediation(med_result) %>% kable_df()

出力で、jobseek_*があるものはmediatedモデルの結果、depress2_*は結果の場合です。 この時点では同じような話ですが、ベイズ的アプローチではもっと楽しいことがあります。 例えば、うつ病の結果の歪度をうまくとらえていないことがわかります。 我々の予測値と観測値は全く一致していません。 媒介因子については少し良くなっていますが、モデルベースの予測ではまだ少し高いかもしれません。

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

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

Pros

• Straightforward syntax
• Extremely powerful- （非常に強力な）。 モデルはほとんど想像力に制限される
• Mediation パッケージが近似するものを基本的に行う
• ベイズ推論のすべての特典がある。 診断、事後予測チェック、モデル比較、など。

Limitations

• Slower to estimate
• ‘By-hand’ calculations needed for going beyond the standard linear model.これは、標準的な線形モデルを超えるために必要な計算です。 しかし、これはベイズの観点からはすでに一般的なアプローチです
• Some comfort with Bayesian approach required

More complexity

Some of the packages mentioned can handle more complex models or provide additional approaches to investigate indirect effects.NetSupport.

相互作用

モデルによっては、媒介モデルや結果に対して相互作用を伴うものがあり、残念ながらこれは媒介モデレーションやモデレートメディテーションと呼ばれることが多いようです。 私自身は、そうでなければ（まだそれほど単純なモデルではないにせよ）単純な概念であるかもしれないものに曖昧な名前をつけることにメリットを感じませんが、その船はずいぶん前に出航しています。 詳細は省きますが、モデルのどこかに相互作用項があり、その相互作用は治療変数、媒介変数、またはその両方を含むかもしれない、ということです。

lm のような標準的なモデリング ツールとその拡張を使用しているため、上記のすべてのパッケージで相互作用を組み込むことは簡単ですが、パスの積のタイプのアプローチは保持されません (a*b != c')。 また、治療、仲介者、結果の間の非線形の関係を調査したい場合もあります。 また、反復測定など相関のある観測データがある場合もあります。 lavaanはバイナリや順序変数のモデルの限られたセットを実際に行うことができますが、適切な間接推定値を得るには、非常に退屈な手作業が必要になります。 あまり数が多くない場合はこれらを削除できることもありますが、それ以外の場合は何か手を打ちたいものです。 lavaan、psych、および brms パッケージは、この状況に対処するための 1 つ以上の方法 (たとえば、多重代入) を提供します。

Alternatives

私たちは、モデルをノード、それらを結ぶ円弧/辺/パスのネットワークとして描写しています。 私たちの議論は、矢印がフィードバック ループのない一方向にしか進まない Directed Acyclic Graphs (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'))# For simpler output we'll use treatment and sex as numeric (explained later)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) )# extract path coefficients if desired# parameters = bn.fit(model, jobs_trim)# parameters$job_seek# parameters$econ_hard# parameters$depress2

library(bnlearn)model = gs(jobs_trim, whitelist = whitelist, blacklist = blacklist)plot(model)

プロットを見ると、少し変化していることがわかりますね。 例えば、年齢は求職の自己効力感にのみ関係し、性別はうつ病にのみ影響します。

もし、パスを以前の例と同じものに限定しても、同じ結果が得られます。

library(bnlearn)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をデフォルトの最尤推定量と一緒に使うのと本質的に同等です。

因子としてtreatmentとsexを使った場合、bnlearnは因子の値によって異なる条件付きモデルを生成することになります。 つまり、treatment == 'treatment'のときとtreatment == controlのときで別々のモデルを持つことになります。 我々の場合、これは、例えば、lm( job_seek ~ treat * (econ_hard + sex + age))のように、すべてが治療と相互作用することを許可することと同じで、同様に、うつ病モデルもそうです。 これは、あらゆるバイナリ変数（例えば、性別を含む）に拡張される可能性があります。 Mediator がバイナリ変数である場合、これはおそらく我々がやりたいことでしょう。

Python

CSCAR director Kerby Shedden は媒介モデルに関する Python ワークショップを行っているので、ここで statsmodels の実装を紹介します。 これはImaiのアプローチに従っているので、mediationパッケージのPython版と見なすことができます。 出力は、基本的に、因子変数として処理を使ったものと同じで、各処理カテゴリについて個別の結果を得ます。 これは我々のデモには不要なので、「平均」結果を以前のmediationパッケージの結果と比較すればよいのです。

import statsmodels.api as smfrom statsmodels.stats.mediation import Mediationimport numpy as npimport pandas as pdoutcome_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-valueACME (control) -0.016 -0.048 0.014 0.332ACME (treated) -0.016 -0.048 0.014 0.332ADE (control) -0.043 -0.130 0.044 0.308ADE (treated) -0.043 -0.130 0.044 0.308Total effect -0.059 -0.144 0.029 0.208Prop. mediated (control) 0.241 -1.710 2.254 0.364Prop. mediated (treated) 0.241 -1.710 2.254 0.364ACME (average) -0.016 -0.048 0.014 0.332ADE (average) -0.043 -0.130 0.044 0.308Prop. mediated (average) 0.241 -1.710 2.254 0.364

Stata

最後に、Stataのsemコマンドを使ってオプションを提供しています。 Stataはこの例で間接効果を簡単に得ることができますが、共変量ごとに行うので、出力は控えめに言っても少し冗長になります6。 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)

Summary

間接効果を持つモデルは、データ分析に採用するために慎重な理論的考察を必要とします。 しかし、モデルがデータの状況に適していれば、Rの様々なパッケージから非常に簡単に結果を得ることができます。さらに、間接効果のある分析を行うために構造方程式モデリングパッケージを使用する必要はなく、実際には、標準のR構文を使用して遠くまで行くことができます。

Package comparison summarized

The following table may help one decide which package to use their needs given their theoretical considerations.とあるように、観測変数、すなわち潜在変数がない場合、SEMツールは必要ないし、推奨さえしない。

1. SEMについては調停分析などもっと詳しいドキュメントがあるんだけどね。︎

2. なぜかあまり実務では見かけませんし、保証されていないのにこのようなモデルに従順なデータにするために何をしたのだろうと思ってしまいます。︎

3. 今井さんの記事はホームページで公開されています。︎

4. MASSはこの時点で10年以上他のものに取って代わられ、ロードするとTidyverseや他のパッケージを混乱させる傾向がほとんどです。 MASS は素晴らしいパッケージですが (昔は素晴らしかった)、パッケージで使用する場合、1つか2つの関数を使用するためだけに環境にロードしない (または他のパッケージをロードしない) ことが良いと思います。 私はほとんどmvrnorm（多変量正規分布）とglm.nbで使われているのを見るだけですが、その機能を持つ他のパッケージがあれば、さらなる利点がありますし、Rコミュニティで最もよく使われているうちの一つであるdplyr関数をマスクしないことです。︎

5. brms はそれに取り組んでいます。︎

6. コード内のオプションは、あり得るものを抑制/最小化するためにあります。︎