Free Press

Table of Contents

Updated August 04, 2020. Der Code kann hier heruntergeladen werden.

Einführung

In manchen Situationen können wir den indirekten Effekt einer Variable auf ein Ergebnis oder Resultat betrachten. So können beispielsweise schlechte Lebensbedingungen im Elternhaus in der Kindheit die Lernergebnisse in der Schule verschlechtern, was sich wiederum negativ auf die spätere Lebensqualität auswirkt, z. B. auf das Lebenseinkommen. In einem anderen Fall könnten wir eine einzelne Variable betrachten, die zu mehreren Zeitpunkten erhoben wird, so dass es einen Effekt der Variable zum Zeitpunkt 1 auf den Zeitpunkt 2 und zum Zeitpunkt 2 auf den Zeitpunkt 3 gibt. Der Grundgedanke ist etwa folgender:

\

Mit anderen Worten: \(\mathcal{A}\) führt zu \(\mathcal{B}\), und dann führt \(\mathcal{B}\) zu \(\mathcal{C}\). Bei Mediationsmodellen setzen wir eine intervenierende Variable zwischen die normale Kovariate \(\rightarrow\) und den Ergebnispfad, den wir im Rahmen der Standardregression haben könnten, und diese Modelle ermöglichen es uns, solche Verhaltensweisen zu untersuchen. In der obigen Darstellung ist die intervenierende Variable oder der Mediator \(\mathcal{B}\). Es ist oft der Fall, dass wir immer noch einen direkten Effekt von \(\mathcal{A}\) auf \(\mathcal{C}\) haben, aber wie bei dem Modell im Allgemeinen wäre dies theoretisch motiviert.

Die Mediationsanalyse ist in den sozialwissenschaftlichen Disziplinen sehr beliebt, obwohl sie keineswegs auf diese beschränkt ist, und wird gewöhnlich unter dem Deckmantel der Strukturgleichungsmodellierung (SEM) durchgeführt, die ihrerseits eine spezifische Ausrichtung von grafischen Modellen im Allgemeinen ist1. Das grafische Modell eines Mediationsmodells könnte wie folgt aussehen.

In diesem Fall spiegeln a und b den indirekten Pfad der Wirkung von \(\mathrm{X}\) auf das Ergebnis über den Mediator wider, während c' der direkte Effekt von \(\mathrm{X}\) auf das Ergebnis ist, nachdem der indirekte Pfad entfernt wurde (c wäre der Effekt vor der Annahme des indirekten Effekts, und c – c' entspricht dem indirekten Effekt). Die Gesamtwirkung von \(\mathrm{X}\) ist die kombinierte indirekte und direkte Wirkung.

Ich sollte ein paar Dinge anmerken, die ich bei meiner Beratungstätigkeit in Dutzenden von Disziplinen sehe. Zunächst einmal scheint es so zu sein, dass nur sehr wenige Leute, die meinen, sie bräuchten ein Mediationsmodell, dies tatsächlich tun. Wenn Sie zum Beispiel Ihr Modell nicht in zeitlichen oder physikalischen Begriffen denken können, so dass \(\mathrm{X}\) notwendigerweise zum Vermittler führt, der dann notwendigerweise zum Ergebnis führt, brauchen Sie wahrscheinlich kein Vermittlungsmodell. Wenn Sie die Pfeile in beide Richtungen sehen können, brauchen Sie wahrscheinlich ebenfalls kein solches Modell. Auch wenn bei der Beschreibung Ihres Modells jeder denkt, dass es sich um eine Interaktion (auch bekannt als Moderation) handelt, brauchen Sie dies wahrscheinlich nicht. Und schließlich, wie man vermuten könnte, wenn es keine starke Korrelation zwischen Schlüsselvariablen (\(\mathrm{X}\)) und Mediator (Pfad a) gibt, und wenn es keine starke Korrelation zwischen Mediator und Ergebnis(sen) gibt (Pfad b), brauchen Sie dies wahrscheinlich nicht. Nichts wird Sie davon abhalten, eine Mediationsanalyse durchzuführen, aber ohne diese Voraussetzungen werden Sie mit ziemlicher Sicherheit ein schwaches und wahrscheinlich verwirrenderes Modell haben, als Sie es sonst hätten.

Kurz gesagt, Mediation funktioniert am besten, wenn es stark implizierte kausale Verbindungen zwischen den Variablen gibt. Selbst dann sollte ein solches Modell mit einem einfacheren Modell ohne Mediation2 verglichen werden. Auf jeden Fall gibt es ein paar sehr einfache Möglichkeiten, solche Modelle in R zu untersuchen, und das ist das Ziel hier, nur um zu demonstrieren, wie man anfangen kann.

Daten

Für die Demonstration von Mediationsmodellen mit den verschiedenen Paketen werden wir den Jobs-Datensatz verwenden, der mit dem Mediationspaket geliefert wird. Hier ist die Beschreibung.

Job Search Intervention Study (JOBS II). JOBS II ist ein randomisiertes Feldexperiment, das die Wirksamkeit einer Berufsbildungsmaßnahme für Arbeitslose untersucht. Das Programm soll nicht nur die Wiederbeschäftigung von Arbeitslosen erhöhen, sondern auch die psychische Gesundheit der Arbeitssuchenden verbessern. Im Rahmen des JOBS-II-Feldversuchs erhielten 1 801 Arbeitslose einen Fragebogen zur Voruntersuchung und wurden dann nach dem Zufallsprinzip einer Behandlungs- und einer Kontrollgruppe zugewiesen. Die Teilnehmer der Behandlungsgruppe nahmen an Workshops zur beruflichen Qualifizierung teil. In den Workshops lernten die Befragten Fähigkeiten zur Arbeitssuche und Bewältigungsstrategien für den Umgang mit Rückschlägen bei der Arbeitssuche. Die Teilnehmer der Kontrollgruppe erhielten eine Broschüre mit Tipps für die Arbeitssuche. Bei den Folgebefragungen waren die beiden wichtigsten Ergebnisvariablen ein kontinuierliches Maß für depressive Symptome auf der Grundlage der Hopkins Symptom Checklist und eine binäre Variable, die angab, ob der Befragte eine Beschäftigung gefunden hatte.

Hier ist eine Beschreibung der Variablen in dieser Demonstration. Es gibt noch weitere Variablen, mit denen Sie vielleicht herumspielen möchten.

• econ_hard: Grad der wirtschaftlichen Notlage vor der Behandlung mit Werten von 1 bis 5.
• sex: Indikatorvariable für das Geschlecht. 1 = weiblich
• age: Alter in Jahren.
• educ: Faktor mit fünf Kategorien für das Bildungsniveau.
• job_seek: Eine kontinuierliche Skala zur Messung des Grades der Selbstwirksamkeit bei der Arbeitssuche mit Werten von 1 bis 5. Die Mediatorvariable.
• depress2: Maß für depressive Symptome nach der Behandlung. Die Ergebnisvariable.
• treat: Indikatorvariable dafür, ob der Teilnehmer zufällig für das JOBS II Trainingsprogramm ausgewählt wurde. 1 = Zuweisung zur Teilnahme.

data(jobs, package = 'mediation')

Modell

Angesichts dieser Daten lauten die Modelle für den Mediator und das Ergebnis wie folgt:

\

Wir erwarten also, dass das Job-Skills-Training einen negativen Effekt auf die Depression hat (d.h. eine Steigerung des Wohlbefindens), aber zumindest ein Teil davon wäre auf einen positiven Effekt auf die Arbeitssuche zurückzuführen.

Als grafisches Modell könnten wir es kurz und bündig wie folgt darstellen.

Pakete

Wir werden uns die folgenden Pakete ansehen, um zu zeigen, wie man Mediationsanalysen in R durchführen kann:

• Mediation
• lavaan
• psych
• brms

Während diese im Mittelpunkt stehen werden, werde ich auch einige andere Alternativen erwähnen, einschließlich Python und Stata.

Mediation

Wir werden mit dem Mediation-Paket beginnen, da es im Grunde nicht mehr Programmierkenntnisse erfordert, als man bereits durch die Ausführung von Standard-Regressionsmodellen in R besitzt. Das Paket liefert den durchschnittlichen kausalen Mediationseffekt, der in der Hilfedatei und in Imais Artikel3 wie folgt definiert ist:

Der durchschnittliche kausale Mediationseffekt (ACME) stellt den erwarteten Unterschied im potenziellen Ergebnis dar, wenn der Mediator den Wert annimmt, der sich unter der Behandlungsbedingung im Gegensatz zur Kontrollbedingung einstellen würde, während der Behandlungsstatus selbst konstant gehalten wird.

Beachten Sie, wie sich diese Definition auf erwartete oder vorhergesagte Werte unter der Bedingung des Behandlungswertes konzentriert. Das Konzept der kontrafaktischen Werte, d. h. wie würde die Beobachtung unter den umgekehrten Bedingungen aussehen, hat in der Modellierung eine lange Tradition. Man kann sich das so vorstellen: Wenn man in der Behandlungsgruppe ist, hat man einen bestimmten Wert für den Mediator, und angesichts dessen hat man dann einen bestimmten erwarteten Wert für das Ergebnis. Wir könnten jedoch die gleiche Beobachtung auch in der Kontrollgruppe machen und den Effekt auf das Ergebnis durch den Mediator genauso bewerten. Wir können die möglichen Ergebnisse bewerten, während wir die Behandlung konstant halten. Die Betrachtung von Ergebnisveränderungen in Abhängigkeit vom Wert des Mediators macht keine Annahmen über den Modelltyp. Auf diese Weise kann das Mediationspaket verschiedene Modelle für den Mediator und das Ergebnis einbeziehen. Zum Beispiel könnte der Mediator binär sein, was ein logistisches Regressionsmodell erfordert, während das Ergebnismodell ein Überlebensmodell sein könnte.

In unserem Beispiel bleiben wir bei standardmäßigen (normalen) linearen Modellen. Beachten Sie auch, dass unsere Behandlung zwar eine binäre Variable ist, dies aber auf den kontinuierlichen Fall verallgemeinert werden kann, bei dem wir das Ergebnis einer Bewegung um eine Einheit auf der „Behandlung“ betrachten. Damit das Mediationspaket funktioniert, führen wir einfach unsere jeweiligen Modelle für den Mediator und das Ergebnis durch und verwenden dann die Funktion mediate, um das Endergebnis zu erhalten.

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)

Die obigen Ergebnisse zeigen, dass sich der ACME statistisch nicht von Null oder keiner Mediation unterscheidet. Der durchschnittliche direkte Effekt ist negativ, aber ebenfalls statistisch nicht bemerkenswert, ebenso wenig wie der Gesamteffekt (indirekter + direkter Effekt). Ebenfalls angegeben ist der so genannte „vermittelte Anteil“, der das Verhältnis des indirekten Effekts zum Gesamteffekt darstellt. Dies ist jedoch keine Proportion, sondern kann sogar negativ sein und ist daher meist eine bedeutungslose Zahl.

Pros

• Standard R-Modelle und -Syntax
• Mehrere Modelltypen für Mediator und Ergebnis
• Liefert mehrere Ergebnisse gleichzeitig
• Gute Dokumentation und zugehörige Artikel sind frei verfügbar
• Kann ‚moderierte‘ Mediation durchführen

Einschränkungen

• Verwendung von MASS4
• Einfache Modelle mit zufälligen Effekten
• Funktionalität kann bei einigen Modellkomplexitäten eingeschränkt sein
• Keine Fähigkeiten für latente Variablen

lavaan

In dem speziellen Fall, in dem sowohl Mediations- als auch Ergebnismodelle lineare Standardmodelle mit einer Normalverteilung für die Zielvariable sind, entspricht der indirekte Effekt dem Produkt aus den Pfaden a und b im vorherigen Diagramm. Der direkte Effekt ist der c'-Pfad. Vergleicht man den alleinigen direkten Effekt, den wir c nennen könnten, mit dem geschätzten direkten Effekt im Mediationsmodell c', so ergibt sich c - c' = a*b. Was bereits erwähnt wurde, wird nun klarer: Wenn entweder a oder b nahezu Null sind, dann kann auch der indirekte Effekt nur nahezu Null sein, so dass es ratsam ist, solche Beziehungen vorher zu untersuchen.

Dieser Produkt-der-Wege-Ansatz (oder Differenz der Koeffizienten) ist derjenige, den wir mit dem lavaan-Paket verfolgen werden, und tatsächlich ist dies zum Zeitpunkt der Erstellung dieses Artikels unsere einzige Möglichkeit, dies zu tun. lavaan ist speziell auf die Modellierung von Strukturgleichungen ausgerichtet, wie z. B. Faktorenanalyse, Wachstumsmodelle und Mediationsmodelle, wie wir sie hier durchführen, und wird für solche Modelle sehr empfohlen. Es ist zwar auf den Standardfall eines linearen Modells zur Bewertung der Mediation beschränkt, aber es ist das einzige unserer Tools, das latente Variablen ohne weiteres einbeziehen kann5. Wir könnten zum Beispiel unser Depressionsergebnis als latente Variable betrachten, die den einzelnen Fragebogenelementen zugrunde liegt. Darüber hinaus könnten wir auch mehrere Mediatoren und mehrere Outcomes einbeziehen.

Um die Dinge so zu halten, wie wir sie besprochen haben, werde ich die a, b und c' Pfade in Lavaan so bezeichnen, wie sie zuvor dargestellt wurden. Ansonsten ist lavaan sehr einfach zu bedienen und verwendet im Falle von beobachteten Variablen die Standard-R-Formelnotation für die Modelle. Darüber hinaus definieren wir die interessierenden Effekte, die wir berechnen wollen, mit dem Operator :=. Wir spezifizieren das Modell in seiner Gesamtheit als einfache Zeichenkette und verwenden dann die Funktion sem, um die Analyse durchzuführen.

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

Wir sehen die gleiche Ausgabe wie zuvor und können unseren Parameter indirect mit dem ACME vergleichen, den wir zuvor hatten, der Effekt direct wird mit dem ADE verglichen, und der total wird mit dem vorherigen Gesamteffekt verglichen. Die Werte sind im Wesentlichen gleich.

Beachten Sie auch, dass die Ausgabe den \(R^2\) Wert für beide Modelle anzeigt. Im Fall von job_seek können wir sehen, dass der Grund, warum wir nicht viel in der Art der Mediation finden, darin liegt, dass die beteiligten Kovariaten keine Variation des Mediators zu Beginn erklären. Eine Voruntersuchung hätte uns in diesem Fall die Mühe erspart.

Pros

• Kann mit mehreren Mediatoren umgehen
• Kann mit mehreren ‚Behandlungen‘ umgehen
• Kann mit mehreren Ergebnissen umgehen
• Kann latente Variablen verwenden
• Einige mehrstufige Unterstützung
• Kann moderierte Mediation und mediated Moderation (allerdings nicht für latente Variablen)

Einschränkungen

• Erfordert zusätzliche Kodierung zur Schätzung des indirekten Effekts
• Einfache zufällige Effekte
• Während die Modelle binäre oder ordinale Variablen für den Mediator/die Ergebnisse einbeziehen können, gibt es keine einfache Möglichkeit, den indirekten Effekt in der Art des Mediationspakets in diesen Fällen zu berechnen.

piecewiseSEM

Das piecewiseSEM-Paket arbeitet sehr ähnlich wie das Mediation-Paket. Das Schöne daran im Vergleich zum Mediation-Paket ist, dass piecewiseSEM mit zusätzlichen Modelltypen umgehen kann, sowie zusätzliche Ausgaben (z.B. standardisierte Ergebnisse), zusätzliche Optionen (z.B. Multigroup, korrelierte Residuen) und eine Visualisierung des Modells bietet.

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

Wir können die Plot-Funktionen nutzen, um eine schnelle Visualisierung des Modells zu erstellen.

plot(mediation_result)

Unglücklicherweise gibt es derzeit keine automatische Möglichkeit, die indirekten Effekte zu berechnen, so dass man die Ergebnisse von Hand booten müsste.

Pros

• Standard R-Modelle und -Syntax
• Mehrere Modelltypen sowohl für Mediator als auch Ergebnis
• Einige SEM-ähnliche Ergebnisse (z.B.. Fit, standardisierte Koeffizienten, AIC)
• Schnelle Darstellung der Ergebnisse
• Kann mit mehreren Mediatoren, ‚Behandlungen‘, und Outcomes

Einschränkungen

• Berechnet nicht automatisch indirekte Effekte
• Keine latenten Variablen

psych

Das Psych-Paket macht sich die Tatsache zunutze, dass man im Falle des linearen Standardmodells die Ergebnisse über die entsprechenden Regressionsmodelle allein auf der Grundlage der Kovarianzmatrizen erhalten kann. Es ist lavaan sehr ähnlich, verwendet aber einen Ansatz der gewöhnlichen kleinsten Quadrate im Gegensatz zur Maximum-Likelihood-Methode. Das Schöne daran ist eine Syntax, die es ermöglicht, sich nur auf den interessierenden Effekt zu konzentrieren oder alles einzubeziehen, was sehr nützlich ist, wenn man auch an den indirekten Effekten für wirtschaftliche Notlagen, Alter und Geschlecht interessiert ist.

Für diese Demo verwenden wir die bereinigte Version mit - anstelle von + für die Nicht-Behandlungseffekte. Dies bedeutet lediglich, dass sie in die Modelle einbezogen werden, aber keine Ergebnisse zu ihnen angezeigt werden. Der Mediator ist mit () gekennzeichnet. Ein weiterer Bonus ist eine schnelle Darstellung der Ergebnisse, die den Unterschied zwischen den unbereinigten und den bereinigten direkten Effekten sowie das entsprechende Bootstrap-Intervall zeigt.

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

Gleiche Ergebnisse, andere Verpackung, aber möglicherweise der bisher einfachste Weg, da nur ein Funktionsaufruf erforderlich ist. Das Psych-Paket kann auch mehrere Vermittler und Ergebnisse als Bonus behandeln.

Pros

• Einfachste Syntax, im Grunde ein Ein-Zeilen-Modell
• Schnelle Darstellung der Ergebnisse
• Kann mit mehreren Mediatoren, ‚Behandlungen‘, und Outcomes
• Kann ‚moderierte‘ Mediation durchführen

Einschränkungen

• Beschränkt auf das lineare Standardmodell (lm)
• Verwendung von MASS

brms

Für unsere nächste Demo kommen wir zu dem meiner Meinung nach mächtigsten Paket, brms. Der Name steht für Bayesian Regression Modeling with Stan, und Stan ist eine leistungsstarke probabilistische Programmiersprache für Bayesianische Analysen. Ich werde nicht näher auf die Bayes’sche Analyse eingehen, aber Sie können sich gerne mein Dokument dazu ansehen.

Im Allgemeinen machen wir es wie bisher, indem wir das Mediatormodell und das Ergebnismodell angeben. brms macht nichts Besonderes für die Mediationsanalyse, aber seine Hypothesenfunktion kann es uns ermöglichen, den Product-of-Paths-Ansatz zu testen. Darüber hinaus liefert das sjstats-Paket die Ergebnisse im Wesentlichen auf die gleiche Weise wie das Mediations-Paket, und im Übrigen ist das Mediations-Paket im Grunde genommen ohnehin ein Versuch einer Bayes’schen Lösung mit frequentistischen Methoden. Hätten wir unterschiedliche Verteilungen für das Ergebnis und den Mediator, wäre es relativ einfach, diese durchschnittlichen Vorhersagewerte und ihre Unterschiede zu erhalten, da Bayes’sche Ansätze immer über posteriore Vorhersageverteilungen nachdenken. Auf jeden Fall ist hier der Code.

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()

In der Ausgabe ist alles mit jobseek_* ein Ergebnis für das Mediatormodell, während depress2_* für das Ergebnis steht. An diesem Punkt haben wir die gleiche alte Geschichte, aber mit dem Bayes’schen Ansatz können wir uns mehr interessante Dinge ansehen. So können wir zum Beispiel sehen, dass wir die Schiefe des Depressionsergebnisses nicht wirklich gut abbilden. Unsere vorhergesagten Werte im Vergleich zu den beobachteten stimmen nicht ganz überein. Beim Mediator sind wir etwas besser, aber bei einigen unserer modellbasierten Vorhersagen vielleicht immer noch ein wenig zu hoch.

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

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

Pros

• Einfache Syntax
• Extrem leistungsstark- Modelle sind größtenteils auf die eigene Vorstellungskraft beschränkt
• Tut im Grunde das, was das Mediation-Paket annähert
• Alle Vorzüge der Bayes’schen Inferenz: Diagnosen, posteriore Vorhersageprüfungen, Modellvergleiche, etc.

Einschränkungen

• Schwächer zu schätzen
• ‚Manuelle‘ Berechnungen, die erforderlich sind, um über das lineare Standardmodell hinauszugehen, aber dies ist bereits ein gängiger Ansatz aus der Bayes’schen Perspektive
• Ein gewisses Maß an Vertrautheit mit dem Bayes’schen Ansatz ist erforderlich

Mehr Komplexität

Einige der genannten Pakete können komplexere Modelle verarbeiten oder bieten zusätzliche Ansätze zur Untersuchung indirekter Effekte.

Interaktionen

Einige Modelle beinhalten Interaktionen entweder für das Mediationsmodell oder für das Ergebnis, und leider wird dies oft als vermittelte Moderation oder moderierte Mediation bezeichnet. Ich persönlich sehe keinen Vorteil darin, einem ansonsten einfachen Konzept (wenn auch nicht so einfachen Modell) mehrdeutige Namen zu geben, aber das Schiff ist schon lange abgefahren. Ich werde nicht ins Detail gehen, aber die Idee ist, dass Sie irgendwo im Modell einen Interaktionsterm haben könnten, und die Interaktion könnte die Behandlungsvariable, den Mediator oder beides betreffen.

Es genügt zu sagen, dass, da wir Standardmodellierungswerkzeuge wie lm und Erweiterungen davon verwenden, die Einbeziehung von Interaktionen für alle der oben genannten Pakete trivial ist, aber der Produkt-der-Wege-Ansatz gilt nicht (a*b != c').

Verallgemeinerte lineare Modelle

In einigen Fällen kann unser Mediator oder Ergebnis binär, zählend oder etwas sein, wo die Annahme einer Normalverteilung nicht die beste Idee ist. Oder wir wollen nichtlineare Beziehungen zwischen Behandlung/Mediator/Ergebnis untersuchen. Oder wir haben Daten mit korrelierten Beobachtungen wie wiederholte Messungen oder ähnliches. Das Mediation-Paket ist besonders stolz darauf, aber brms kann alles, was es kann, und noch mehr, auch wenn Sie möglicherweise etwas mehr Arbeit leisten müssen, um das Ergebnis tatsächlich zu berechnen. lavaan kann tatsächlich eine begrenzte Anzahl von Modellen für binäre und ordinale Variablen durchführen, aber um die geeignete indirekte Schätzung zu erhalten, wäre ein sehr mühsamer manueller Ansatz erforderlich.

Fehlende Daten

Wenn man mit solchen Daten zu tun hat, vor allem in den Sozialwissenschaften, fehlen oft Daten zu einer der Kovariaten. Manchmal können wir diese weglassen, wenn es nicht zu viele sind, aber in anderen Fällen werden wir etwas dagegen tun wollen. Die Pakete lavaan, psych und brms bieten eine oder mehrere Möglichkeiten, mit dieser Situation umzugehen (z.B. multiple Imputation).

Alternativen

Wir haben die Modelle als Netzwerke von Knoten dargestellt, mit Bögen/Kanten/Pfaden, die sie verbinden. Unsere Diskussion dreht sich um so genannte gerichtete azyklische Graphen (DAG), bei denen die Pfeile nur in eine Richtung gehen können und es keine Rückkopplungsschleifen gibt. Das Ergebnis einer beliebigen Ergebnisvariablen ist eine Funktion der ihr vorangehenden Pfeile und bedingt unabhängig von anderen. Einige theoretische Modelle können dies lockern, und andere können überhaupt keine Pfeile haben, d.h. ungerichtet sein, so dass wir nur an den Verbindungen interessiert sind (z.B. bei einigen sozialen Netzwerken).

bnlearn

Das bnlearn-Paket ermöglicht die Untersuchung von gerichteten, teilweise gerichteten und ungerichteten Graphen. In Bezug auf DAGs können wir es verwenden, um die Vermittlungsmodelle, die wir diskutiert haben, im Wesentlichen zu duplizieren. Das Schöne daran ist, dass dieses Paket Pfade effizient auf ihre Einbeziehung testet, anstatt sie anzunehmen, aber wir können immer noch theoretische Beschränkungen auferlegen, wenn es nötig ist. Wir können dann nicht nur auf prinzipielle Weise mit Bayes’schen Netzwerken und Pearls kausaler Graphentheorie als Grundlage nach den interessierenden Pfaden suchen, sondern wir haben auch Werkzeuge, um eine Überanpassung durch Kreuzvalidierung zu vermeiden.

Für das erste Modell stellen wir sicher, dass Pfade zwischen Behandlung – Mediator, Behandlung – Ergebnis und Mediator – Ergebnis existieren (die Whitelist). Unsinnige Pfade wie Pfeile zur Behandlung (die zufällig zugewiesen wurde), zum Geschlecht, zur wirtschaftlichen Notlage und zum Alter (die schwarze Liste) werden nicht zugelassen. Ansonsten werden wir sehen, was die Daten nahelegen.

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)

Wir sehen in der Grafik, dass sich die Dinge ein wenig verändert haben. Zum Beispiel bezieht sich das Alter jetzt nur noch auf die Selbstwirksamkeit bei der Arbeitssuche, und das Geschlecht wirkt sich nur noch auf die Depression aus.

Wenn wir die Pfade so einschränken, dass sie nur das sind, was sie in unseren vorherigen Beispielen sind, erhalten wir die gleichen Ergebnisse.

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 

Das Wichtigste ist, dass die geschätzten Parameter dem entsprechen, was wir mit den vorherigen Paketen erhalten haben. Es ist im Wesentlichen äquivalent zur Verwendung von lavaan mit dem Standard-Maximum-Likelihood-Schätzer.

Wenn wir Behandlung und Geschlecht als Faktoren verwenden, wird bnlearn bedingte Modelle erzeugen, die sich je nach dem gewählten Faktorwert unterscheiden. Mit anderen Worten, man hätte ein separates Modell für wenn treatment == 'treatment' und eines für wenn treatment == control. In unserem Fall wäre dies gleichbedeutend damit, alles mit der Behandlung interagieren zu lassen, z. B. lm( job_seek ~ treat * (econ_hard + sex + age)), und dies gilt auch für das Depressionsmodell. Dies würde für jede binäre Variable gelten (z. B. auch für das Geschlecht). Wenn der Mediator eine binäre Variable ist, ist dies wahrscheinlich das, was wir tun wollen.

Python

CSCAR-Direktor Kerby Shedden hat einen Python-Workshop über Mediationsmodelle gegeben, daher zeige ich hier die Implementierung von statsmodels. Sie folgt dem Imai-Ansatz und kann daher als die Python-Version des Mediation-Pakets angesehen werden. Die Ausgabe ist im Wesentlichen die gleiche wie bei der Verwendung der Behandlung als Faktorvariable, bei der Sie separate Ergebnisse für jede Behandlungskategorie erhalten. Dies ist für unsere Demo nicht notwendig, so dass Sie einfach die „durchschnittlichen“ Ergebnisse mit den vorherigen Ergebnissen des Mediationspakets vergleichen können.

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

Schließlich biete ich eine Option in Stata mit dem Befehl sem. Stata macht es einfach, die indirekten Effekte in diesem Beispiel zu erhalten, aber es tut dies für jede Kovariate, so dass die Ausgabe gelinde gesagt ein wenig langatmig ist6. Wer mit Stata arbeitet, braucht kein separates SEM-Paket, um diese Art von Ergebnissen zu erhalten.

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)

Zusammenfassung

Modelle mit indirekten Effekten erfordern sorgfältige theoretische Überlegungen, um sie für die Datenanalyse zu verwenden. Wenn das Modell jedoch für Ihre Datensituation geeignet ist, ist es recht einfach, mit einer Vielzahl von Paketen in R Ergebnisse zu erhalten. Außerdem muss man kein Strukturgleichungsmodellierungspaket verwenden, um eine Analyse mit indirekten Effekten durchzuführen, und man kann sogar mit der Standard-R-Syntax weit kommen. Für rein beobachtete, d. h. nicht latente Variablen ist kein SEM-Tool erforderlich oder sogar empfehlenswert.

\

Paketvergleich zusammengefasst

Die folgende Tabelle kann bei der Entscheidung helfen, welches Paket angesichts der theoretischen Überlegungen für die eigenen Bedürfnisse zu verwenden ist.

1. Ich habe ein viel ausführlicheres Dokument über SEM, einschließlich Mediationsanalyse.︎

2. Aus irgendeinem Grund sieht man das in der Praxis nicht oft, und man fragt sich, was getan wurde, um die Daten für ein solches Modell zugänglich zu machen, wenn es nicht gerechtfertigt war.︎

3. Imai stellt seine Artikel auf seiner Website zur Verfügung.︎

4. MASS ist seit über einem Jahrzehnt durch andere Pakete ersetzt worden, und es neigt dazu, Tidyverse und andere Pakete durcheinander zu bringen, wenn es geladen wird. Es ist ein gutes Paket (und war früher großartig), aber wenn Sie es in einem Paket verwenden wollen, wäre es gut, es (oder andere Pakete) nicht in der Umgebung zu laden, nur um eine oder zwei Funktionen zu verwenden. Ich sehe es meistens nur für mvrnorm (multivariate Normalverteilung) und glm.nb verwendet, aber es gibt andere Pakete mit dieser Funktionalität, die zusätzliche Vorteile bieten würden, und nicht die dplyr-Funktionen maskieren, die zu den am häufigsten verwendeten in der R-Gemeinschaft gehören.︎

5. brms arbeitet daran.︎

6. Die Optionen im Code sind dazu da, zu unterdrücken/minimieren, was sein kann.︎