Free Press

Table of Contents

Updated August 04, 2020. Code kan hier worden gedownload.

Inleiding

In sommige situaties kunnen we het indirecte effect van een variabele op een uitkomst of resultaat in overweging nemen. Slechte levensomstandigheden thuis in de kindertijd kunnen bijvoorbeeld leiden tot slechtere leerresultaten op school, die vervolgens een negatief effect hebben op de latere kwaliteit van het leven, bijvoorbeeld het levensinkomen. In een ander geval kan het gaan om één variabele die op meerdere tijdstippen wordt verzameld, zodat er een effect is van de variabele op tijdstip 1 op tijdstip 2, en tijdstip 2 op tijdstip 3. Het basisidee is ongeveer als volgt:

Met andere woorden, \(\mathcal{A}}) leidt tot \(\mathcal{B}}), en vervolgens leidt \(\mathcal{B}}) tot \(\mathcal{C}}). Met bemiddelingsmodellen stellen we een interveniërende variabele tussen het normale covariaat en het uitkomstpad dat we in de standaard regressie setting zouden kunnen hebben, en deze modellen stellen ons in staat dergelijk gedrag te onderzoeken. In het bovenstaande is de interveniërende variabele, of mediator, de covariabele. Het is vaak nog mogelijk dat er een direct effect is van \(\mathcal{A}} op \(\mathcal{C}}), maar zoals met het model in het algemeen, zou dit theoretisch gemotiveerd zijn.

Mediatie-analyse is zeer populair in de sociale wetenschappen, hoewel geenszins beperkt tot die disciplines, en wordt gewoonlijk uitgevoerd onder het mom van structurele vergelijking modellering (SEM), die zelf een specifieke oriëntatie is van grafische modellen meer in het algemeen1. Het grafische model van een bemiddelingsmodel zou er als volgt uit kunnen zien.

In dit geval geven a en b het indirecte pad weer van het effect van \(\mathrm{X}\) op de uitkomst via de mediator, terwijl c' het directe effect is van \(\mathrm{X}\) op de uitkomst nadat het indirecte pad is verwijderd (c zou het effect zijn voordat het indirecte effect werd vastgesteld, en c – c' is gelijk aan het indirecte effect). Het totale effect van c is het gecombineerde indirecte en directe effect.

Ik moet een paar dingen opmerken op basis van wat ik zie in consulting in tientallen disciplines. Om te beginnen lijkt het erop dat maar heel weinig mensen die denken dat ze een bemiddelingsmodel nodig hebben, dat ook werkelijk hebben. Als je bijvoorbeeld niet in staat bent je model in temporele of fysische termen te zien, zodat het noodzakelijkerwijs leidt tot de bemiddelaar, die vervolgens noodzakelijkerwijs leidt tot het resultaat, dan heb je waarschijnlijk geen bemiddelingsmodel nodig. Als je de pijlen in beide richtingen kunt zien gaan, heb je zo’n model waarschijnlijk ook niet nodig. Ook als iedereen bij de beschrijving van uw model denkt dat u het over een interactie hebt (ook wel moderatie genoemd), hebt u dit wellicht niet nodig. En tenslotte, zoals men zou kunnen vermoeden, als er geen sterke correlatie is tussen de sleutelvariabelen (pad a) en de mediator (pad b), dan heb je dit waarschijnlijk niet nodig. Hoewel niets u ervan weerhoudt om een bemiddelingsanalyse uit te voeren, zult u zonder dergelijke voorwaarden vrijwel zeker een zwak en waarschijnlijk verwarrender model hebben dan u anders zou hebben.

Kortom, bemiddeling werkt het best als er sterk geïmpliceerde causale verbanden tussen de variabelen zijn. Zelfs dan moet een dergelijk model worden vergeleken met een eenvoudiger model van geen bemiddeling2. In elk geval zijn er een paar zeer eenvoudige manieren om dergelijke modellen in R te onderzoeken, en dat is het doel hier, gewoon om te demonstreren hoe je aan de slag kunt.

Data

Voor de demonstratie van bemiddelingsmodellen met de verschillende pakketten, zullen we de jobs-dataset gebruiken die bij het bemiddelingspakket wordt geleverd. Hier volgt de beschrijving.

Job Search Intervention Study (JOBS II). JOBS II is een gerandomiseerd veldexperiment dat de doeltreffendheid onderzoekt van een jobtraininginterventie voor werkloze werknemers. Het programma is ontworpen om niet alleen de herintreding van werklozen te verhogen, maar ook om de geestelijke gezondheid van de werkzoekenden te verbeteren. In het praktijkexperiment JOBS II ontvingen 1.801 werklozen een vragenlijst vóór de screening en werden zij vervolgens willekeurig ingedeeld in een behandelings- en controlegroep. Degenen in de behandelingsgroep namen deel aan workshops in beroepsvaardigheden. In de workshops leerden de respondenten vaardigheden voor het zoeken naar werk en strategieën om met tegenslagen in het sollicitatieproces om te gaan. De deelnemers in de controlegroep ontvingen een boekje met tips voor het zoeken naar werk. In de follow-upinterviews waren de twee belangrijkste uitkomstvariabelen een continue meting van depressieve symptomen op basis van de Hopkins Symptom Checklist, en een binaire variabele, die aangaf of de respondent werk had gevonden.

Hier volgt een beschrijving van de variabelen in deze demonstratie. Er zijn nog andere beschikbaar waarmee u misschien ook wilt spelen.

• econ_hard: Niveau van economische hardheid vóór de behandeling met waarden van 1 tot 5.
• sex: Indicatorvariabele voor geslacht. 1 = vrouwelijk
• leeftijd: Leeftijd in jaren.
• educ: Factor met vijf categorieën voor opleidingsniveau.
• job_seek: Een continue schaal die het niveau van zelfeffectiviteit bij het zoeken naar werk meet met waarden van 1 tot 5. De mediërende variabele.
• depress2: Meting van depressieve symptomen na de behandeling. De uitkomstvariabele.
• treat: Indicatorvariabele voor de vraag of de deelnemer willekeurig was geselecteerd voor het JOBS II-opleidingsprogramma. 1 = toewijzing aan deelname.

data(jobs, package = 'mediation')

Model

Gezien deze gegevens zijn de modellen voor de mediator en het resultaat als volgt:

Dus verwachten we dat de training in beroepsvaardigheden een negatief effect heeft op depressie (d.w.z. een toename in welzijn), maar ten minste een deel daarvan zou te danken zijn aan een positief effect op het zoeken naar werk.

In een grafisch model zouden we dit als volgt beknopt kunnen weergeven.

Pakketten

We zullen de volgende pakketten bekijken om te laten zien hoe men bemiddelingsanalyses in R kan uitvoeren:

• mediation
• lavaan
• psych
• brms

Hoewel deze de focus zullen zijn, zal ik ook enkele andere alternatieven vermelden, waaronder Python en Stata.

mediation

We beginnen met het mediation pakket, omdat het in principe niet meer programmeerkennis vereist om uit te voeren dan men al bezit van het uitvoeren van standaard regressie modellen in R. Het pakket levert het gemiddelde causale bemiddelingseffect, dat als volgt wordt gedefinieerd uit het helpbestand en de artikelen van Imai3:

Het gemiddelde causale bemiddelingseffect (ACME) staat voor het verwachte verschil in de potentiële uitkomst wanneer de bemiddelaar de waarde aanneemt die zich onder de behandelingsconditie zou realiseren ten opzichte van de controleconditie, terwijl de behandelingsstatus zelf constant wordt gehouden.

Merk op dat deze definitie is toegespitst op verwachte of voorspelde waarden die voorwaardelijk zijn voor de behandelingswaarde. Deze notie van contrafeiten, of hoe zou de waarneming eruit zien onder de tegenovergestelde setting, heeft een lange geschiedenis in de modellering op dit punt. Zie het zo: als iemand in de behandelingsgroep zit, zou hij een specifieke waarde voor de mediator hebben, en, gegeven dat, zou hij dan een specifieke verwachte waarde voor het resultaat hebben. Wij kunnen echter dezelfde waarneming doen als iemand uit de controlegroep, en het effect op het resultaat door de bemiddelaar op dezelfde manier beoordelen. Wij kunnen de potentiële uitkomsten beoordelen terwijl wij de behandeling constant houden. Het denken over uitkomstveranderingen gegeven de waarde van de mediator maakt geen veronderstelling over het modeltype. Dit is hoe het mediatiepakket in staat is om verschillende modellen voor de mediator versus de uitkomst te incorporeren. De mediator zou bijvoorbeeld binair kunnen zijn, waarvoor een logistisch regressiemodel nodig is, terwijl het uitkomstmodel een overlevingsmodel zou kunnen zijn.

In ons voorbeeld zullen we het houden bij standaard (normale) lineaire modellen. Merk ook op dat, hoewel onze behandeling een binaire variabele is, dit generaliseert naar het continue geval, waarin we het resultaat van een beweging van één eenheid op de “behandeling” beschouwen. Om het bemiddelingspakket te laten werken, voeren we gewoon onze respectieve modellen voor de bemiddelaar en het resultaat uit, en gebruiken dan de functie bemiddelen om het eindresultaat te krijgen.

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)

De bovenstaande resultaten tonen aan dat de ACME statistisch niet verschilt van nul, of geen bemiddeling. Het gemiddelde directe effect is negatief, maar ook dit is statistisch niet opmerkelijk, evenmin als het totale effect (indirect + direct effect). Ook wordt de soi disante “proportion mediated” gegeven, dat wil zeggen de verhouding van het indirecte effect tot het totale effect. Dit is echter geen verhouding, en kan zelfs negatief zijn, en is dus meestal een nietszeggend getal.

Pros

• Standaard R-modellen en syntax
• Meerdere typen modellen voor zowel mediator als uitkomst
• Levert meerdere resultaten tegelijk
• Goede documentatie en bijbehorende artikelen zijn vrij beschikbaar
• Kan ‘gemodereerde’ mediatie doen

Beperkingen

• Gebruik van MASS4
• Eenvoudige random effects modellen
• Functionaliteit wellicht beperkt bij sommige modelcomplexiteiten
• Geen latente variabele mogelijkheden

lavaan

In het specifieke geval dat zowel bemiddelings- als uitkomstmodellen standaard lineaire modellen zijn met een normale verdeling voor de doelvariabele, is het indirecte effect gelijk aan het product van de paden a en b in het vorige diagram. Het directe effect is het pad c'. Een vergelijking van het op zichzelf staande directe effect, dat we c zouden kunnen noemen, met dit geschatte directe effect in het bemiddelingsmodel c', is zodanig dat c - c' = a*b. Wat eerder werd gezegd, is nu misschien duidelijker: als a of b bijna nul zijn, kan het indirecte effect alleen maar bijna nul zijn, zodat het verstandig is dergelijke relaties van tevoren te onderzoeken.

Deze product-of-paths (of verschil in coëfficiënten) benadering is die welke we zullen volgen met het lavaan pakket, en in feite, vanaf dit schrijven, is dat onze enige manier om het te doen. lavaan is specifiek gericht op structurele vergelijking modellering, zoals factoranalyse, groeimodellen, en bemiddelingsmodellen zoals we hier uitvoeren, en wordt sterk aanbevolen voor dergelijke modellen. Hoewel het beperkt is tot het standaard lineaire model om bemiddeling te beoordelen, is het de enige van onze tools die gemakkelijk latente variabelen kan opnemen5. Wij zouden bijvoorbeeld onze depressie-uitkomst als latente variabele kunnen opnemen die ten grondslag ligt aan de individuele vragenlijstitems. Bovendien zouden we ook meerdere mediatoren en meerdere uitkomsten kunnen opnemen.

Om de zaken te houden zoals we ze hebben besproken, zal ik de paden a, b en c' in lavaan labelen volgens de manier waarop ze eerder zijn afgebeeld. Voor de rest is lavaan zeer eenvoudig te gebruiken, en in het geval van waargenomen variabelen, gebruikt het standaard R formule notatie voor de modellen. Daarbuiten definiëren we de effecten van belang die we willen berekenen met de := operator. We specificeren het model in zijn geheel als een eenvoudige tekenreeks, en gebruiken dan de sem-functie om de analyse uit te voeren.

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

We zien dezelfde uitvoer als voorheen en kunnen onze indirect parameter vergelijken met de ACME die we eerder hadden, het direct effect wordt vergeleken met de ADE, en de total vergelijkt met het vorige totale effect. De waarden zijn in wezen gelijk.

Merk ook op dat de uitvoer de R^2 waarde voor beide modellen laat zien. In het geval van job_seek kunnen we zien dat de reden waarom we niet veel bemiddeling vinden, is dat de betrokken covariaten om te beginnen geen variatie in de bemiddelaar verklaren. Vooronderzoek zou ons in dit geval de moeite hebben bespaard.

Pros

• Kan meerdere mediatoren aan
• Kan meerdere ‘behandelingen’ aan
• Kan meerdere uitkomsten aan
• Kan latente variabelen gebruiken
• Bevat enige ondersteuning voor multilevel
• Kan gemodereerde mediatie doen en gemedieerde moderatie (maar niet voor latente variabelen)

Limitations

• Vraagt extra codering om het indirecte effect te schatten
• Single random effects
• Weliswaar kunnen de modellen binaire of ordinale variabelen voor de mediator/uitkomsten bevatten, is er geen eenvoudige manier om het indirecte effect op de wijze van het mediatiepakket in die settings te berekenen.

piecewiseSEM

Het piecewiseSEM pakket werkt zeer vergelijkbaar met het mediation pakket. Het aardige hiervan ten opzichte van het mediation pakket is dat piecewiseSEM extra typen modellen aankan, alsmede extra output kan leveren (b.v. gestandaardiseerde resultaten), extra opties (b.v. multigroep, gecorreleerde residuen), en visualisatie van het model.

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

We kunnen de plotmogelijkheden gebruiken om het model snel te visualiseren.

plot(mediation_result)

Er is momenteel helaas geen automatische manier om de indirecte effecten te berekenen, zodat men de resultaten met de hand zou moeten bootstrappen.

Pros

• Standaard R-modellen en syntax
• Meerdere soorten modellen voor zowel mediator als uitkomst
• Een aantal resultaten in SEM-stijl (bijv. fit, gestandaardiseerde coëfficiënten, AIC)
• Snelle plot van resultaten
• Kan meerdere mediatoren, ‘behandelingen’ en uitkomsten

Beperkingen

• Berekent niet automatisch indirecte effecten
• Geen latente variabele mogelijkheden

psych

Het psych pakket maakt gebruik van het feit dat in het geval van een standaard lineair model, men de resultaten kan verkrijgen via de juiste regressiemodellen op basis van de covariantiematrices alleen. Het is zeer vergelijkbaar met lavaan, hoewel het een gewone kleinste kwadraten aanpak gebruikt in tegenstelling tot maximum likelihood. Het mooie hier is een syntaxis die je toelaat je alleen te concentreren op het effect dat je interesseert, of alles mee te nemen, wat leuk is als je ook geïnteresseerd bent in de indirecte effecten voor economische tegenspoed, leeftijd, en geslacht.

Voor deze demo gebruiken we de opgeschoonde versie met -, in plaats van +, voor de niet-behandelingseffecten. Dit betekent alleen dat ze in de modellen zijn opgenomen, maar dat de resultaten niet worden getoond. De mediator wordt aangeduid met (). Een andere bonus is een snelle plot van de resultaten, die het verschil laat zien tussen de niet-gecorrigeerde en de gecorrigeerde directe effecten, en het bijbehorende bootstrapped interval.

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

Zelfde resultaten, andere verpakking, maar mogelijk de gemakkelijkste route tot nu toe, omdat er slechts één functie-aanroep voor nodig was. Het psych pakket behandelt ook meerdere bemiddelaars en uitkomsten als een bonus.

Pros

• Gemakkelijkste syntaxis, in feite een model van één regel
• Snel plot van resultaten
• Kan meerdere mediatoren, ‘behandelingen’, en uitkomsten
• Kan ‘gemodereerde’ mediatie doen

Limitaties

• Beperkt tot standaard lineair model (lm)
• Gebruik van MASS

brms

Voor onze volgende demo komen we bij wat ik het krachtigste pakket vind, brms. De naam staat voor Bayesian Regression Modeling with Stan, en Stan is een krachtige probabilistische programmeertaal voor Bayesiaanse analyse. Ik zal niet in detail treden over Bayesiaanse analyse, maar voel je vrij om mijn document te bekijken dat dat wel doet.

In het algemeen doen we zoals we eerder hebben gedaan, het specificeren van het mediator model en het uitkomst model. brms doet niets speciaals voor mediatie analyse, maar zijn hypothese functie kan ons in staat stellen om de product-of-paths benadering te testen. Bovendien zal het sjstats pakket in wezen de resultaten leveren op dezelfde manier als het mediatie pakket dat voor ons doet, en wat dat betreft is het mediatie pakket in feite toch al een poging tot een Bayesiaanse oplossing met behulp van frequentistische methoden. Als we verschillende verdelingen voor de uitkomst en de mediator zouden hebben, zouden we relatief gemakkelijk deze gemiddelde voorspellingswaarden en hun verschillen kunnen verkrijgen, aangezien Bayesiaanse benaderingen altijd denken aan posterior voorspellingsverdelingen. Hoe dan ook, hier is de 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 de uitvoer is alles met jobseek_* een resultaat voor het mediator model, terwijl depress2_* voor de uitkomst is. Op dit punt hebben we hetzelfde oude verhaal, maar met de Bayesiaanse benadering hebben we leukere dingen om naar te kijken. We kunnen bijvoorbeeld zien dat we de scheefheid van de depressie-uitkomst niet goed vastleggen. Onze voorspelde waarden versus de waargenomen waarden komen niet helemaal overeen. We zijn iets beter voor de mediator, maar misschien nog steeds een beetje te hoog met sommige van onze modelmatige voorspellingen.

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

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

Pros

• Rechte syntaxis
• Extreem krachtig- Modellen zijn meestal beperkt tot iemands verbeelding
• Doet in feite wat het mediation pakket benadert
• Alle voordelen van Bayesiaanse inferentie: diagnostiek, posterior predictieve controles, modelvergelijking, enz.

Limitaties

• Lager te schatten
• ‘Met de hand’ berekeningen nodig om verder te gaan dan het standaard lineaire model, maar dit is al een gebruikelijke aanpak vanuit het Bayesiaanse perspectief
• Enig comfort met de Bayesiaanse aanpak vereist

Meer complexiteit

Enkele van de genoemde pakketten kunnen complexere modellen aan of bieden extra benaderingen om indirecte effecten te onderzoeken.

Interacties

Sommige modellen omvatten interacties, hetzij voor het bemiddelingsmodel, hetzij voor de uitkomst, en helaas wordt dit vaak aangeduid als gemedieerde moderatie of gemodereerde bemiddeling. Persoonlijk zie ik er het voordeel niet van in om dubbelzinnige namen te geven aan wat anders een rechtlijnig concept zou kunnen zijn (zij het nog steeds een niet zo rechtlijnig model), maar dat schip is al lang geleden vertrokken. Ik ga niet in op de details, maar het idee is dat je ergens in het model een interactieterm kunt hebben, en de interactie kan betrekking hebben op de behandelingsvariabele, de mediator, of beide.

Het volstaat te zeggen dat, aangezien we standaard modelleergereedschappen als lm en uitbreidingen daarvan gebruiken, het opnemen van interacties voor alle bovenstaande pakketten triviaal is, maar dat de product-van-paden-benadering niet opgaat (a*b != c').

Generaliseerde lineaire modellen

In sommige gevallen kan onze mediator of uitkomst binair zijn, tellen, of iets waarbij het aannemen van een normale verdeling misschien niet het beste idee is. Of we willen niet-lineaire verbanden tussen de behandeling/bemiddelaar/uitkomst onderzoeken. Of we kunnen gegevens hebben die gecorreleerde waarnemingen bevatten, zoals herhaalde metingen of iets dergelijks. Het mediation pakket is hier trots op, maar brms kan alles wat het kan en meer, hoewel je misschien iets meer werk moet doen om het resultaat daadwerkelijk te berekenen. lavaan kan eigenlijk een beperkte set modellen doen voor binaire en ordinale variabelen, maar het verkrijgen van de juiste indirecte schatting zou een zeer vervelende handmatige aanpak vereisen.

Missende gegevens

Vaak wanneer we te maken hebben met dergelijke gegevens, vooral in de sociale wetenschappen, ontbreken vaak gegevens over een van de covariaten. Soms kunnen we die weglaten als het er niet te veel zijn, maar in andere gevallen zullen we er toch iets aan willen doen. De pakketten lavaan, psych, en brms bieden een of meer manieren om met deze situatie om te gaan (b.v. meervoudige imputatie).

Alternatieven

We hebben de modellen afgebeeld als netwerken van knooppunten, met bogen/edges/paden die ze verbinden. Onze discussie draait om zogeheten Directed Acyclic Graphs (DAG), waarbij de pijlen slechts in één richting kunnen gaan en er geen terugkoppelingslussen zijn. Het resultaat van elke uitkomstvariabele is een functie van de pijlen die eraan voorafgaan, en voorwaardelijk onafhankelijk van andere. Sommige theoretische modellen kunnen dit versoepelen, en andere kunnen helemaal geen pijlen hebben, d.w.z. ongericht zijn, zodat we alleen geïnteresseerd zijn in de verbindingen (b.v. met sommige sociale netwerken).

bnlearn

Het bnlearn pakket laat onderzoek toe van gerichte, gedeeltelijk gerichte, en ongerichte grafieken. In termen van DAGs, kunnen we het gebruiken om in wezen de bemiddelingsmodellen te dupliceren die we hebben besproken. Het mooie is echter dat dit pakket op efficiënte wijze paden toetst op inclusie in plaats van ze aan te nemen, maar we kunnen nog steeds theoretische beperkingen opleggen indien nodig. Niet alleen kunnen we dan op een principiële manier zoeken naar de paden die van belang zijn met Bayesiaanse netwerken en Pearl’s causale grafentheorie als basis, we zullen ook hulpmiddelen hebben om overfitting verder te vermijden via cross-validatie.

Voor het initiële model zullen we ervoor zorgen dat er paden bestaan tussen behandeling – mediator, behandeling – uitkomst, en mediator – uitkomst (de whitelist). Onzinnige paden, zoals pijlen naar de behandeling (die willekeurig werd toegewezen), geslacht, economische tegenspoed en leeftijd, verbieden we (de zwarte lijst). Voor het overige zullen we zien wat de gegevens suggereren.

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)

We zien in de plot dat de zaken enigszins zijn veranderd. Zo heeft leeftijd nu alleen nog maar betrekking op de zelfeffectiviteit bij het zoeken naar werk, en heeft geslacht alleen nog maar een effect op depressie.

Als we de paden zouden beperken tot wat ze in onze vorige voorbeelden zijn, zouden we dezelfde resultaten krijgen.

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 

Het belangrijkste om op te merken is dat de geschatte parameters gelijk zijn aan wat we met de vorige pakketten kregen. Het is in wezen gelijk aan het gebruik van lavaan met de standaard maximum likelihood estimator.

Als we behandeling en geslacht als factoren gebruiken, zal bnlearn voorwaardelijke modellen produceren die verschillend zijn, afhankelijk van de factorwaarde die wordt genomen. Met andere woorden, er zou een afzonderlijk model zijn voor treatment == 'treatment' en één voor treatment == control. In ons geval zou dit identiek zijn aan het toestaan van een wisselwerking tussen alles en behandeling, bv. lm( job_seek ~ treat * (econ_hard + sex + age)), en evenzo voor het depressiemodel. Dit zou kunnen worden uitgebreid tot elke binaire variabele (bv. ook geslacht). Als de mediator een binaire variabele is, is dit waarschijnlijk wat we zouden willen doen.

Python

CSCAR-directeur Kerby Shedden heeft een Python-workshop gegeven over bemiddelingsmodellen, dus ik laat hier de statsmodels-implementatie zien. Het volgt de Imai benadering en kan dus gezien worden als de Python versie van het mediation pakket. De output is in wezen hetzelfde als wat je zou krijgen met behandeling als een factorvariabele, waar je aparte resultaten krijgt voor elke behandelingscategorie. Dit is niet nodig voor onze demo, dus u kunt gewoon de ‘gemiddelde’ resultaten vergelijken met de vorige mediation package resultaten.

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

Ten slotte geef ik een optie in Stata met behulp van zijn sem commando. Stata maakt het gemakkelijk om de indirecte effecten in dit voorbeeld te krijgen, maar het doet dat voor elke covariaat, dus de uitvoer is op zijn zachtst gezegd een beetje langdradig6. Wie met Stata werkt, heeft geen apart SEM-pakket nodig om dit soort resultaten te verkrijgen.

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)

Samenvatting

Modellen met indirecte effecten vereisen een zorgvuldige theoretische afweging om ze voor de gegevensanalyse te kunnen gebruiken. Als het model echter geschikt is voor uw gegevenssituatie, is het vrij eenvoudig om resultaten te verkrijgen met behulp van diverse pakketten in R. Bovendien hoeft men geen pakket voor structurele vergelijkingmodellering te gebruiken om een analyse met indirecte effecten uit te voeren, en in feite kan men ver komen met behulp van de standaard R-syntaxis. Voor strikt waargenomen, d.w.z. geen latente variabelen, is geen SEM-gereedschap nodig, of zelfs aanbevolen.

Vergelijking van pakketten samengevat

De volgende tabel kan iemand helpen beslissen welk pakket te gebruiken voor zijn behoeften, gegeven zijn theoretische overwegingen.

1. Ik heb een veel gedetailleerder document over SEM, met inbegrip van bemiddelingsanalyse.︎

2. Om de een of andere reden zie je dit in de praktijk niet veel, en je vraagt je af wat er is gedaan om de gegevens voor zo’n model geschikt te maken als het niet gerechtvaardigd was.︎

3. Imai stelt zijn artikelen beschikbaar op zijn website.︎

4. MASS wordt al meer dan tien jaar door anderen vervangen en heeft meestal alleen maar de neiging om je tidyverse en andere pakketten in de war te sturen als het wordt geladen. Het is een prima pakket (en was geweldig in de tijd), maar als je het in een pakket wilt gebruiken, zou het goed zijn om het niet te laden (of andere pakketten) in de omgeving alleen om een functie of twee te gebruiken. Ik zie het meestal alleen gebruikt voor mvrnorm (multivariate normale verdeling) en glm.nb, maar er zijn andere pakketten met die functionaliteit die extra voordelen zouden bieden, en niet dplyr functies maskeren, die tot de meest gebruikte in de R gemeenschap behoren.︎

5. brms is ermee bezig.︎

6. De opties in de code zijn er om te onderdrukken/minimaliseren wat kan.︎