Free Press

Materialförteckning

Uppdaterad 04 augusti 2020. Koden kan laddas ner här.

Inledning

I vissa situationer kan vi överväga den indirekta effekten av någon variabel på ett utfall eller resultat. Som exempel kan dåliga levnadsförhållanden hemma under barndomen försämra inlärningsresultaten i skolan, vilket senare har en negativ effekt på senare livskvalitet, t.ex. livsinkomstinkomster. I ett annat fall kan vi överväga en enda variabel som samlas in vid flera tidpunkter, så att det finns en effekt av variabeln vid tid 1 på tid 2 och tid 2 på tid 3. Den grundläggande idén är ungefär följande:

\

Med andra ord leder \(\mathcal{A}\) till \(\mathcal{B}\), och sedan leder \(\mathcal{B}\) till \(\mathcal{C}\). Med hjälp av medieringsmodeller antar vi en intervenerande variabel mellan den normala kovarianten \(\rightarrow\) och den utfallsväg som vi kan ha i den vanliga regressionsmiljön, och dessa modeller gör det möjligt för oss att undersöka sådana beteenden. I ovanstående modell är den intervenerande variabeln, eller mediatorn, \(\mathcal{B}\). Det är ofta så att vi fortfarande kan ha en direkt effekt av \(\mathcal{A}\) på \(\mathcal{C}\), men som med modellen i allmänhet skulle detta vara teoretiskt motiverat.

Mediationsanalyser är mycket populära inom samhällsvetenskapliga discipliner, även om de inte på något sätt är begränsade till dessa, och utförs vanligen under täckmanteln strukturella ekvationsmodeller (SEM), som i sig självt är en specifik inriktning av grafiska modeller mer allmänt1. Den grafiska modellen för en medieringsmodell kan se ut på följande sätt.

I detta fall återspeglar a och b den indirekta vägen för effekten av \(\mathrm{X}\) på utfallet genom mediatorn, medan c' är den direkta effekten av \(\mathrm{X}\) på utfallet efter att den indirekta vägen har tagits bort (c skulle vara effekten före den indirekta effekten, och c – c' är lika med den indirekta effekten). Den totala effekten av \(\mathrm{X}\) är den kombinerade indirekta och direkta effekten.

Jag bör notera några saker baserat på vad jag ser vid konsultationer inom dussintals discipliner. Till att börja med verkar det vara väldigt få personer som tror att de behöver en medieringsmodell som faktiskt gör det. Om du till exempel inte kan tänka på din modell i tidsmässiga eller fysiska termer, så att \(\mathrm{X}\) nödvändigtvis leder till medlaren, som sedan nödvändigtvis leder till utfallet, behöver du troligen ingen medlingsmodell. Om du kan se att pilarna går åt båda hållen behöver du troligen inte heller en sådan modell. Om alla när de beskriver din modell tror att du talar om en interaktion (även kallad moderering) behöver du kanske inte heller detta. Och slutligen, som man kan misstänka, om det inte finns någon stark korrelation mellan nyckelvariabler (\(\mathrm{X}\)) och mediator (path a), och om det inte finns någon stark korrelation mellan mediator och utfallet/utfallen (path b), behöver du förmodligen inte detta. Även om ingenting hindrar dig från att göra en medlingsanalys, kommer du utan sådana förutsättningar nästan säkert att få en svag och förmodligen mer förvirrande modell än vad du annars skulle ha fått.

Samt sett fungerar medling bäst när det finns starkt underförstådda kausala samband mellan variablerna. Även då bör en sådan modell jämföras med en enklare modell utan mediering2. I vilket fall som helst finns det några mycket enkla sätt att undersöka sådana modeller i R, och det är målet här, att bara demonstrera hur du kan komma igång.

Data

För att demonstrera medieringsmodeller med de olika paketen kommer vi att använda datasetet Jobs som följer med paketet Mediation. Här är beskrivningen.

Job Search Intervention Study (JOBS II). JOBS II är ett randomiserat fältexperiment som undersöker effekten av en intervention för arbetsträning för arbetslösa arbetstagare. Programmet är utformat för att inte bara öka återanställningen bland arbetslösa utan också förbättra den mentala hälsan hos de arbetssökande. I fältförsöket JOBS II fick 1 801 arbetslösa arbetstagare ett frågeformulär för förhandsgranskning och fördelades sedan slumpmässigt mellan behandlings- och kontrollgrupper. De som ingick i behandlingsgruppen deltog i seminarier om yrkeskunskaper. I workshoparna fick de svarande lära sig färdigheter i jobbsökning och strategier för att hantera motgångar i jobbsökningsprocessen. De som ingick i kontrollgruppen fick ett häfte med tips om jobbsökning. I uppföljningsintervjuer var de två viktigaste utfallsvariablerna ett kontinuerligt mått på depressiva symtom baserat på Hopkins Symptom Checklist och en binär variabel som representerade huruvida respondenten hade fått arbete.

Här finns en beskrivning av variablerna i denna demonstration. Det finns andra tillgängliga som du kanske också vill leka med.

• econ_hard: Nivå av ekonomisk svårighet före behandling med värden från 1 till 5.
• sex: Indikatorvariabel för kön. 1 = kvinna
• ålder: Ålder i år.
• educ: Faktor med fem kategorier för utbildningsnivå.
• job_seek: En kontinuerlig skala som mäter nivån av självförtroende för jobbsökning med värden från 1 till 5. Mediatorvariabel.
• depress2: Mått på depressiva symtom efter behandlingen. Utfallsvariabel.
• treat: Indikatorvariabel för om deltagaren slumpmässigt valdes ut till utbildningsprogrammet JOBS II. 1 = Tilldelning till deltagande.

data(jobs, package = 'mediation')

Modell

Genom dessa data är modellerna för mediatorn och utfallet följande:

\

Därmed förväntar vi oss att utbildningen i yrkesskicklighet ska ha en negativ effekt på depression (dvs. en ökning av välbefinnandet), men åtminstone en del av detta skulle bero på en positiv effekt på jobbsökning.

Som en grafisk modell kan vi kortfattat beskriva det på följande sätt.

Paket

Vi kommer att titta på följande paket för att visa hur man kan genomföra medlingsanalyser i R:

• mediation
• lavaan
• psych
• brms

Men även om dessa kommer att vara i fokus kommer jag också att nämna några andra alternativ, inklusive Python och Stata.

mediation

Vi kommer att börja med mediationspaketet, eftersom det i princip inte kräver mer programmeringsförmåga för att genomföra än vad man redan besitter från att köra vanliga regressionsmodeller i R. Paketet tillhandahåller den genomsnittliga kausala medieringseffekten, som definieras enligt följande från hjälpfilen och Imais artiklar3:

Den genomsnittliga kausala medieringseffekten (ACME) representerar den förväntade skillnaden i det potentiella utfallet när medlaren tog det värde som skulle förverkligas under behandlingstillståndet i motsats till kontrolltillståndet, samtidigt som behandlingstillståndet i sig hålls konstant.

Bemärk hur den här definitionen är inriktad på förväntade eller förutspådda värden som är betingade av behandlingsvärdet. Denna föreställning om kontrafaktiska faktorer, eller hur observationen skulle se ut under den motsatta inställningen, har en lång historia inom modellering vid denna tidpunkt. Tänk på det så här: Om man ingår i behandlingsgruppen skulle man ha ett specifikt värde för mediatorn, och med tanke på detta skulle man ha ett specifikt förväntat värde för utfallet. Vi skulle dock kunna anta att samma observatör också ingår i kontrollgruppen och bedöma effekten på utfallet genom medlaren på samma sätt. Vi kan bedöma de potentiella resultaten samtidigt som vi håller behandlingen konstant. Att tänka på förändringar i utfallet med tanke på mediatorns värde innebär inget antagande om modelltypen. Det är på detta sätt som medlingspaketet kan införliva olika modeller för mediatorn respektive utfallet. Till exempel kan mediatorn vara binär, vilket kräver en logistisk regressionsmodell, medan utfallsmodellen kan vara en överlevnadsmodell.

I vårt exempel kommer vi att hålla oss till vanliga (normala) linjära modeller. Observera också att även om vår behandling är en binär variabel, generaliseras detta till det kontinuerliga fallet, där vi betraktar resultatet av en förflyttning av en enhet på ”behandlingen”. För att mediationspaketet ska fungera kör vi helt enkelt våra respektive modeller för mediatorn och resultatet och använder sedan funktionen mediate för att få fram slutresultatet.

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)

Resultaten ovan visar att ACME inte skiljer sig statistiskt från noll, eller ingen mediering. Den genomsnittliga direkta effekten är negativ men likaså inte statistiskt anmärkningsvärd, inte heller den totala effekten (indirekt + direkt effekt). Dessutom anges soi disant ”proportion mediated”, vilket är förhållandet mellan den indirekta effekten och den totala effekten. Detta är dock inte en andel, och kan till och med vara negativt, och är därför oftast ett meningslöst tal.

Pros

• Standard R-modeller och syntax
• Multipla typer av modeller för både mediator och utfall
• Givar flera resultat samtidigt
• God dokumentation och tillhörande artiklar finns fritt tillgängliga
• Kan göra ”modererad” mediering

Begränsningar

.
• Användning av MASS4
• Enkla modeller med slumpmässiga effekter
• Funktionaliteten kan vara begränsad med vissa komplexa modeller
• Ingen kapacitet för latenta variabler

lavaan

I det specifika fallet där både medlings- och utfallsmodellerna är linjära standardmodeller med en normalfördelning för målvariabeln, är den indirekta effekten likvärdig med produkten av a och b i föregående diagram. Den direkta effekten är c'. En jämförelse mellan den fristående direkta effekten, som vi kan kalla c, och denna uppskattade direkta effekt i medlingsmodellen c', är sådan att c - c' = a*b. Det som nämndes tidigare kan nu bli tydligare, om antingen a eller b är nästan noll kan den indirekta effekten bara vara nästan noll, så det är klokt att undersöka sådana samband i förväg.

Denna produkt-of-paths (eller skillnad i koefficienter) metod är den vi kommer att använda med lavaan-paketet, och i själva verket är det, i skrivande stund, vårt enda sätt att gå till väga. lavaan är särskilt inriktad på strukturell ekvationsmodellering, såsom faktoranalyser, tillväxtmodeller och medieringsmodeller som vi genomför här, och rekommenderas starkt för sådana modeller. Även om det är begränsat till standardfallet med linjära modeller för att bedöma mediering, är det det enda av våra verktyg som lätt kan införliva latenta variabler5. Vi skulle till exempel kunna ha vårt utfall för depression som en latent variabel som ligger till grund för de enskilda enkätfrågorna. Dessutom skulle vi också kunna införliva flera mediatorer och flera utfall.

För att hålla saker och ting som vi har diskuterat kommer jag att märka a, b och c' stigarna i lavaan enligt hur de har avbildats tidigare. I övrigt är lavaan mycket lätt att använda, och när det gäller observerade variabler används standard R-formelnotation för modellerna. Utöver detta definierar vi de effekter av intresse som vi vill beräkna med operatören :=. Vi anger modellen i sin helhet som en enkel teckensträng och använder sedan funktionen sem för att göra analysen.

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

Vi ser samma utdata som tidigare och kan jämföra vår indirect-parameter med ACME som vi hade tidigare, direct-effekten jämförs med ADE och total jämförs med den tidigare totala effekten. Värdena är i huvudsak desamma.

Bemärk också att utmatningen visar \(R^2\)-värdet för båda modellerna. I fallet job_seek kan vi se att anledningen till att vi inte hittar mycket i form av mediering är att de inblandade kovariaterna inte förklarar någon variation i mediatorn till att börja med. En preliminär undersökning skulle ha besparat oss besväret i det här fallet.

Pros

• Kan hantera flera mediatorer
• Kan hantera flera ”behandlingar”
• Kan hantera flera utfall
• Kan använda latenta variabler
• En del stöd för flera nivåer
• Kan göra modererad medling och medierad. moderation (dock inte för latenta variabler)

Begränsningar

• Kräver ytterligare kodning för att skatta den indirekta effekten
• Enkla slumpmässiga effekter
• Men modellerna skulle kunna innehålla binära eller ordinalvariabler för mediatorn/utfallen, finns det inget enkelt sätt att beräkna den indirekta effekten på samma sätt som medlingspaketet i dessa fall.

piecewiseSEM

Paketet piecewiseSEM fungerar mycket likt mediationspaketet. Det fina med detta i förhållande till mediationspaketet är att piecewiseSEM kan hantera ytterligare typer av modeller, samt ge ytterligare utdata (t.ex. standardiserade resultat), ytterligare alternativ (t.ex. multigrupp, korrelerade residualer) och visualisering av modellen.

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

Vi kan använda dess plottningsmöjligheter för att skapa en snabb visualisering av modellen.

plot(mediation_result)

Det finns tyvärr inget automatiskt sätt att beräkna indirekta effekter i dagsläget, så man måste bootstrappa resultaten för hand.

Pros

• Standard R-modeller och syntax
• Multipla typer av modeller för både mediator och utfall
• En del resultat i SEM-stil (t.ex. fit, standardiserade koefficienter, AIC)
• Snabb plottning av resultat
• Kan hantera flera mediatorer, ”behandlingar”, och utfall

Begränsningar

• Beräknar inte automatiskt indirekta effekter
• Ingen kapacitet för latenta variabler

psyk

Paketet psyk drar nytta av det faktum att i fallet med den linjära standardmodellen kan man få fram resultaten via lämpliga regressionsmodeller baserade enbart på kovariansmatriserna. Det är mycket likt lavaan, även om det använder en vanlig minsta kvadratmetoden i motsats till maximal sannolikhet. Det fina här är en syntax som gör att man kan fokusera endast på den effekt man är intresserad av eller inkludera allt, vilket är trevligt om man är intresserad av de indirekta effekterna för ekonomiska svårigheter, ålder och kön också.

För den här demonstrationen använder vi den renodlade versionen med -, istället för +, för icke-behandlingseffekter. Detta innebär bara att de ingår i modellerna, men att resultaten inte visas när det gäller dem. Mediatorn identifieras med (). En annan bonus är en snabb plott av resultaten som visar skillnaden mellan de ojusterade och justerade direkta effekterna och det lämpliga bootstrappade intervallet.

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

Samma resultat, annan förpackning, men möjligen den enklaste vägen hittills eftersom den bara krävde ett funktionsanrop. Psykpaketet hanterar också flera medlare och resultat som en bonus.

Pros

• Enklaste syntaxen, i princip en enradsmodell
• Snabb plottning av resultat
• Kan hantera flera mediatorer, ”behandlingar”, och utfall
• Kan göra ’modererad’ mediering

Begränsningar

• Begränsad till linjär standardmodell (lm)
• Användning av MASS

brms

För vår nästa demo kommer vi till det som jag tycker är det mest kraftfulla paketet, brms. Namnet står för Bayesian Regression Modeling with Stan, och Stan är ett kraftfullt probabilistiskt programmeringsspråk för bayesiansk analys. Jag kommer inte att gå in på detaljer om Bayesiansk analys, men se gärna mitt dokument som gör det.

Vi gör i allmänhet som vi har gjort tidigare och specificerar mediatormodellen och utfallsmodellen. brms gör inget speciellt för mediationsanalys, men dess hypotesfunktion kan göra det möjligt för oss att testa produkt-av-stigar-metoden. Dessutom kommer sjstats-paketet i huvudsak att ge resultaten på samma sätt som mediationspaketet gör för oss, och för den delen är mediationspaketet i princip ett försök till en bayesiansk lösning med hjälp av frekventistiska metoder ändå. Om vi verkligen hade olika fördelningar för utfallet och medlaren skulle vi ha relativt lätt att få fram dessa genomsnittliga prediktionsvärden och deras skillnader, eftersom bayesianska metoder alltid tänker på efterföljande prediktionsfördelningar. Här är i alla fall koden.

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

I resultatet är allt med jobseek_* ett resultat för mediatormodellen, medan depress2_* är för resultatet. Vi har samma gamla historia vid det här laget, men med den bayesianska metoden har vi roligare saker att titta på. Vi kan till exempel se att vi faktiskt inte fångar skevheten i depressionsutfallet på ett bra sätt. Våra förutspådda värden jämfört med de observerade stämmer inte riktigt överens. Vi är lite bättre för mediatorn, men kanske fortfarande lite för högt med vissa av våra modellbaserade förutsägelser.

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

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

Pros

• Straightforward syntax
• Extremt kraftfull- Modellerna begränsas mestadels av ens fantasi
• Gör i princip det som mediationspaketet approximerar
• Alla fördelar med Bayesiansk inferens: Diagnostik, efterföljande prediktiva kontroller, jämförelse av modeller osv.

Begränsningar

• Snålare skattning
• ”By hand”-beräkningar som behövs för att gå utöver den linjära standardmodellen, men detta är redan ett vanligt tillvägagångssätt från det bayesianska perspektivet
• En viss bekvämlighet med det bayesianska tillvägagångssättet krävs

Mer komplexitet

Vissa av de nämnda paketen kan hantera mer komplexa modeller eller tillhandahålla ytterligare tillvägagångssätt för att undersöka indirekta effekter.

Interaktioner

Vissa modeller inbegriper interaktioner antingen för medieringsmodellen eller utfallet, och tyvärr kallas detta ofta för medierad moderering eller modererad mediering. Personligen ser jag inte fördelen med att ge tvetydiga namn åt något som annars skulle kunna vara ett okomplicerat begrepp (även om det fortfarande inte är en lika okomplicerad modell), men det skeppet har seglat för länge sedan. Jag tänker inte gå in på detaljerna, men tanken är att du kan ha en interaktionsterm någonstans i modellen, och interaktionen kan inbegripa behandlingsvariabeln, mediatorn eller båda.

Det räcker med att säga att eftersom vi använder standardmodelleringsverktyg som lm och utvidgningar av det, är det trivialt att införliva interaktioner för alla ovanstående paket, men att den typ av tillvägagångssätt som bygger på produkten av vägar inte gäller (a*b != c').

Generaliserade linjära modeller

I vissa fall kan vår mediator eller vårt utfall vara binär, räkna, eller något annat där det kanske inte är den bästa idén att anta en normalfördelning. Eller så kanske vi vill undersöka icke-linjära samband mellan behandling/mediator/utfall. Eller så kanske vi har data med korrelerade observationer, t.ex. upprepade mätningar eller liknande. Mediationspaketet är stolt över just detta, men brms kan göra allt det kan göra och mer, även om du kanske måste göra lite mer arbete för att faktiskt beräkna resultatet. lavaan kan faktiskt göra en begränsad uppsättning modeller för binära och ordinala variabler, men för att få fram den lämpliga indirekta skattningen skulle det krävas ett mycket tråkigt tillvägagångssätt för hand.

Missing data

Ofta när vi har att göra med sådana data, särskilt inom samhällsvetenskaperna, saknas det ofta data om någon av kovariaterna. Ibland kan vi släppa dessa om det inte är för många, men i andra fall vill vi göra något åt saken. Paketen lavaan, psych och brms tillhandahåller ett eller flera sätt att hantera situationen (t.ex. multipel imputering).

Alternativ

Vi har avbildat modellerna som nätverk av noder, med bågar/kanter/stigar som förbinder dem. Vår diskussion kretsar kring vad som kallas riktade acykliska grafer (Directed Acyclic Graphs, DAG) där pilarna endast kan gå i en riktning utan återkopplingsslingor. Resultatet av en utfallsvariabel är en funktion av de pilar som föregår den och är villkorligt oberoende av andra. Vissa teoretiska modeller kan lätta på detta, och andra kan ha inga pilar alls, dvs. vara oriktade, så att vi bara är intresserade av kopplingarna (t.ex. med vissa sociala nätverk).

bnlearn

Paketet bnlearn gör det möjligt att undersöka riktade, delvis riktade och oriktade grafer. När det gäller DAGs kan vi använda det för att i huvudsak duplicera de förmedlingsmodeller som vi har diskuterat. Det fina är dock att det här paketet effektivt kommer att testa vägar för inkludering i stället för att anta dem, men vi kan fortfarande införa teoretiska begränsningar vid behov. Vi kan då inte bara söka efter de intressanta vägarna på ett principiellt sätt med bayesianska nätverk och Pearls kausala grafteori som grund, vi kommer också att ha verktyg för att ytterligare undvika överanpassning via korsvalidering.

För den första modellen kommer vi att se till att det finns vägar mellan behandling – mediator, behandling – utfall och mediator – utfall (whitelistan). Vi kommer att förbjuda nonsensiska vägar som att ha pilar till behandlingen (som tilldelades slumpmässigt), kön, ekonomiska svårigheter och ålder (den svarta listan). Annars får vi se vad data tyder på.

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)

Vi ser i diagrammet att saker och ting har förändrats en aning. Till exempel relaterar ålder nu bara till jobbsökande self-efficacy och kön har bara en effekt på depression.

Om vi begränsar stigarna till att bara vara vad de är i våra tidigare exempel får vi samma resultat.

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 

Det viktigaste att notera är att de skattade parametrarna är lika med samma sak som vi fick med tidigare paket. Det är i huvudsak likvärdigt med att använda lavaan med standard maximum likelihood estimatorn.

Om vi använder behandling och kön som faktorer kommer bnlearn att producera villkorliga modeller som är olika beroende på vilket faktorvärde som tas. Med andra ord skulle man ha en separat modell för när treatment == 'treatment' och en för när treatment == control. I vårt fall skulle detta vara identiskt med att låta allt interagera med behandling, t.ex. lm( job_seek ~ treat * (econ_hard + sex + age)), och likaså för depressionsmodellen. Detta skulle utvidgas till potentiellt alla binära variabler (t.ex. inklusive kön). Om mediatorn är en binär variabel är detta sannolikt vad vi skulle vilja göra.

Python

CSCAR-chefen Kerby Shedden har gett en Python-workshop om medieringsmodeller, så jag visar statsmodels implementering här. Den följer Imai-metoden och kan därför ses som Python-versionen av mediationspaketet. Resultatet är i huvudsak detsamma som det du skulle få om du använde behandling som en faktorvariabel, där du får separata resultat för varje behandlingskategori. Detta är onödigt för vår demo, så du kan bara jämföra de ”genomsnittliga” resultaten med resultaten från det tidigare mediationspaketet.

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

Slutligt tillhandahåller jag ett alternativ i Stata med hjälp av dess sem-kommando. Stata gör det enkelt att få fram de indirekta effekterna i det här exemplet, men det gör det för varje kovariat, så utmatningen är minst sagt lite spretig6. De som arbetar med Stata behöver inte ett separat SEM-paket för att få den här typen av resultat.

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)

Sammanfattning

Modeller med indirekta effekter kräver noggranna teoretiska överväganden för att användas för dataanalys. Om modellen är lämplig för din datasituation är det dock ganska enkelt att få fram resultat från en mängd olika paket i R. Dessutom behöver man inte använda ett paket för strukturell ekvationsmodellering för att genomföra en analys med indirekta effekter, och i själva verket kan man komma långt med hjälp av R:s standardsyntax. För strikt observerade, dvs. inga latenta, variabler är inget SEM-verktyg nödvändigt, eller till och med rekommenderat.

\

Paketjämförelse sammanfattat

Följande tabell kan hjälpa en att avgöra vilket paket man ska använda för sina behov med tanke på sina teoretiska överväganden.

1. Jag har ett mycket mer detaljerat dokument om SEM, inklusive medlingsanalys.︎

2. För någon anledning ser man inte detta i praktiken så mycket, och man undrar vad som gjordes för att göra data mottagliga för en sådan modell om det inte var befogat. ︎

3. Imai gör sina artiklar tillgängliga på sin webbplats.︎

4. MASS har ersatts av andra i över ett decennium vid det här laget, och det tenderar för det mesta bara att förstöra ditt tidyverse och andra paket när det laddas. Det är ett bra paket (och var bra på den tiden), men om du vill använda det i ett paket skulle det vara bra att inte ladda det (eller andra paket) i miljön bara för att använda en funktion eller två. Jag ser oftast bara att det används för mvrnorm (multivariat normalfördelning) och glm.nb, men det finns andra paket med den funktionen som skulle ge ytterligare fördelar, och inte maskera dplyr-funktioner, som är bland de mest använda i R-communityt.︎

5. brms arbetar på det. ︎

6. Optionerna i koden är till för att undertrycka/minimera vad som kan vara. ︎