Free Press

Tartalomjegyzék

Frissítve: 2020. augusztus 04., 04. A kód letölthető innen.

Bevezetés

Az egyes helyzetekben figyelembe vehetjük valamilyen változó közvetett hatását egy kimenetelre vagy eredményre. Például a gyermekkori rossz otthoni életkörülmények csökkenthetik az iskolai tanulási eredményeket, ami később negatívan hat a későbbi életminőségre, például az élethosszig tartó jövedelemre. Egy másik esetben figyelembe vehetünk egyetlen változót, amelyet több időpontban gyűjtöttünk, úgy, hogy a változónak van hatása az 1. időpontban a 2. időpontra, és a 2. időpontban a 3. időpontra. Az alapötlet valami ilyesmi:

\

Más szóval, \(\mathcal{A}\) vezet \(\mathcal{B}\), majd \(\mathcal{B}\) vezet \(\mathcal{C}\). A mediációs modellekkel egy közbeeső változót tételezünk fel a normál kovariáns \(\rightarrow\) kimeneti útja között, amelyet a standard regressziós környezetben kaphatunk, és ezek a modellek lehetővé teszik számunkra az ilyen viselkedések vizsgálatát. A fentiekben a beavatkozó változó vagy mediátor \(\mathcal{B}\). Gyakran előfordul, hogy még mindig lehet \(\mathcal{A}\) közvetlen hatása \(\mathcal{C}\), de mint a modell esetében általában, ez is elméletileg motivált lenne.

A mediációs elemzés nagyon népszerű a társadalomtudományi tudományágakban, bár korántsem csak azokban, és általában a strukturális egyenletmodellezés (SEM) leple alatt végzik, amely maga is a grafikus modellek egy speciális irányzata általánosabban1. Egy mediációs modell grafikus modellje a következőképpen nézhet ki.

Ebben az esetben a és b a közvetítőn keresztül a \(\mathrm{X}\) eredményre gyakorolt hatásának közvetett útját tükrözik, míg c' az \(\mathrm{X}\) közvetlen hatása a kimenetelre a közvetett út eltávolítása után (c lenne a közvetett hatás megállapítása előtti hatás, és c – c' egyenlő a közvetett hatással). Az \(\mathrm{X}\) teljes hatása a közvetett és közvetlen hatások kombinációja.

Meg kell jegyeznem néhány dolgot az alapján, amit a több tucat tudományágban végzett tanácsadás során látok. Először is, úgy tűnik, hogy nagyon kevés embernek van szüksége közvetítési modellre, aki azt hiszi, hogy szüksége van rá. Például, ha nem tudja időbeli vagy fizikai értelemben elképzelni a modelljét, úgy, hogy \(\mathrm{X}\) szükségszerűen a mediátorhoz vezet, amely aztán szükségszerűen a kimenethez vezet, akkor valószínűleg nincs szüksége mediációs modellre. Ha a nyilak bármelyik irányba haladnak, megint csak valószínűleg nincs szükséged ilyen modellre. Továbbá, ha a modell leírásakor mindenki azt hiszi, hogy kölcsönhatásról (más néven moderációról) beszélsz, akkor valószínűleg nincs szükséged ilyenre. És végül, ahogyan azt gyanítani lehet, ha nincs erős korreláció a kulcsváltozók (\(\(\mathrm{X}\)) és a mediátor (út a) között, és ha nincs erős korreláció a mediátor és az eredmény(ek) között (út b), akkor valószínűleg nincs szükséged erre. Bár semmi sem akadályozza meg abban, hogy mediációs elemzést végezzen, ilyen előfeltételek nélkül szinte biztosan gyenge és valószínűleg zavarosabb modellt kap, mint amilyet egyébként kapna.

Röviden, a mediáció akkor működik a legjobban, ha a változók között erősen implikált oksági kapcsolatok vannak. Az ilyen modellt még ekkor is össze kell hasonlítani a közvetítés nélküli egyszerűbb modellel2. Mindenesetre van néhány nagyon egyszerű módja az ilyen modellek vizsgálatának az R-ben, és ez itt a cél, csak bemutatni, hogyan kezdhetjük el.

Adatok

A mediációs modellek bemutatásához a különböző csomagokkal a mediációs csomaghoz tartozó Jobs adathalmazt fogjuk használni. Íme a leírás.

Job Search Intervention Study (JOBS II). A JOBS II egy randomizált terepkísérlet, amely egy álláskeresési beavatkozás hatékonyságát vizsgálja munkanélküli munkavállalókon. A program célja, hogy ne csak a munkanélküliek újbóli elhelyezkedését növelje, hanem az álláskeresők mentális egészségét is javítsa. A JOBS II terepkísérletben 1801 munkanélküli munkavállaló kapott egy előzetes szűrő kérdőívet, majd véletlenszerűen beosztották őket kezelési és kontrollcsoportokba. A kezelési csoportba tartozók munkahelyi készségfejlesztő workshopokon vettek részt. A workshopokon a válaszadók álláskeresési készségeket és megküzdési stratégiákat tanultak az álláskeresési folyamat során felmerülő kudarcok kezelésére. A kontrollcsoportba tartozók egy füzetet kaptak, amely álláskeresési tanácsokat tartalmazott. A nyomonkövetési interjúk során a két fő kimeneti változó a depressziós tünetek folyamatos mérése volt a Hopkins tünetellenőrzési lista alapján, valamint egy bináris változó, amely azt jelezte, hogy a válaszadó elhelyezkedett-e.

A következőkben a demonstrációban szereplő változók leírása következik. Mások is rendelkezésre állnak, amelyekkel szintén érdemes lehet játszani.

• econ_hard: A kezelés előtti gazdasági nehézségek szintje 1 és 5 közötti értékekkel.
• sex: A nemre vonatkozó indikátor változó. 1 = nő
• age: Életkor években.
• educ: Az iskolai végzettségre vonatkozó öt kategóriát tartalmazó faktor.
• job_seek: Folyamatos skála, amely az álláskeresési önhatékonyság szintjét méri 1-től 5-ig terjedő értékekkel. A közvetítő változó.
• depress2: A kezelés utáni depressziós tünetek mérőszáma. A kimeneti változó.
• treat: Indikátor változó arra vonatkozóan, hogy a résztvevőt véletlenszerűen választották-e ki a JOBS II képzési programba. 1 = a részvételre való kijelölés.

data(jobs, package = 'mediation')

Modell

Az adatok ismeretében a mediátorra és a kimenetre vonatkozó modellek a következők:

\

Ezért azt várjuk, hogy a munkahelyi készségfejlesztő képzés negatív hatást gyakorol a depresszióra (azaz növeli a jóllétet), de ennek legalább egy része az álláskeresésre gyakorolt pozitív hatásnak lesz köszönhető.

Grafikus modellként ezt tömören a következőképpen ábrázolhatnánk.

Pakettek

A következő csomagokat fogjuk megvizsgálni, hogy bemutassuk, hogyan lehet mediációs elemzést végezni R-ben:

• mediáció
• lavaan
• psych
• brms

Míg ezeken lesz a hangsúly, megemlítek néhány más alternatívát is, beleértve a Pythont és a Statát.

mediáció

A mediációs csomaggal fogunk kezdeni, mivel alapvetően nem igényel több programozási képességet az elvégzéséhez, mint amennyivel az ember már rendelkezik a standard regressziós modellek R-ben történő futtatásából. A csomag biztosítja az átlagos oksági mediációs hatást, amelyet a súgófájlból és Imai cikkeiből3 a következőképpen definiálunk:

Az átlagos oksági mediációs hatás (ACME) a potenciális kimenetel várható különbségét jelenti, amikor a mediátor azt az értéket vette fel, amely a kezelési feltétel mellett realizálódna a kontrollállapottal szemben, miközben magát a kezelési állapotot állandónak tartjuk.

Megjegyezzük, hogy ez a definíció a kezelés értékétől függő várható vagy előrejelzett értékekre összpontosít. A kontrafaktumoknak ez a fogalma, vagy hogy hogyan nézne ki a megfigyelés az ellenkező beállítás mellett, hosszú múltra tekint vissza a modellezésben ezen a ponton. Gondoljon erre így: ha valaki a kezelési csoportban van, akkor egy adott értéket kapna a közvetítőre vonatkozóan, és ezt figyelembe véve egy adott várható értéket kapna az eredményre vonatkozóan. Azonban ugyanazt a megfigyelést úgy is feltehetjük, mintha a kontrollcsoportban lennénk, és ugyanúgy értékelhetjük a közvetítőn keresztül a kimenetelre gyakorolt hatást. Értékelhetjük a lehetséges kimeneteleket, miközben a kezelést állandónak tartjuk. Az eredményváltozásra való gondolkodás a mediátor értékének ismeretében nem tesz feltevést a modelltípusra vonatkozóan. Így a mediációs csomag képes különböző modelleket beépíteni a mediátorra vs. a kimenetelre vonatkozóan. Például a mediátor lehet bináris, ami logisztikus regressziós modellt igényel, míg a kimeneti modell lehet egy túlélési modell.

Példánkban maradunk a standard (normál) lineáris modelleknél. Vegyük észre azt is, hogy bár a kezelésünk bináris változó, ez általánosítható a folytonos esetre, ahol a “kezelés” egy egységnyi mozgásának eredményét tekintjük. Ahhoz, hogy a mediációs csomag működjön, egyszerűen lefuttatjuk a megfelelő modelljeinket a mediátorra és a kimenetelre, majd a mediate függvénnyel megkapjuk a végeredményt.

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)

A fenti eredmények azt mutatják, hogy az ACME statisztikailag nem különbözik a nullától vagy a közvetítés hiányától. Az átlagos közvetlen hatás negatív, de szintén nem statisztikailag figyelemre méltó, ahogy a teljes hatás (közvetett + közvetlen hatás) sem. Szintén megadva van a soi disant “közvetített arány”, amely a közvetett hatás és az összhatás aránya. Ez azonban nem arány, és akár negatív is lehet, így többnyire értelmetlen szám.

Előnyei

• Szokásos R modellek és szintaxis
• Sokféle modelltípus mind a mediátor, mind a kimenetel esetében
• Egyszerre több eredményt szolgáltat
• Jó dokumentáció és a kapcsolódó cikkek szabadon elérhetőek
• Képes a “moderált” mediációra

Korlátozások

• A MASS4 használata
• Egyszerű véletlen hatású modellek
• A funkcionalitás talán korlátozott bizonyos modellkomplexitással
• Nincs látens változó képesség

lavaan

A speciális esetben, amikor mind a mediációs, mind a kimeneti modellek standard lineáris modellek, normál eloszlású célváltozóval, a közvetett hatás egyenértékű az előző ábrán szereplő a és b ösvények szorzatával. A közvetlen hatás a c' ösvény. Az önálló közvetlen hatás összehasonlítása, amelyet c-nek nevezhetünk, vs. ez a becsült közvetlen hatás a közvetítési modellben c', olyan, hogy c - c' = a*b. Amit korábban említettünk, most talán világosabbá válik, ha akár a a, akár a b közel nulla, akkor a közvetett hatás is csak közel nulla lehet, ezért az ilyen összefüggéseket célszerű előzetesen megvizsgálni.

Ezt az utak termékének (vagy az együtthatók különbségének) megközelítését fogjuk alkalmazni a lavaan csomaggal, sőt, jelen írásunk pillanatában ez az egyetlen módszerünk. a lavaan kifejezetten a strukturális egyenletmodellezésre, például faktorelemzésre, növekedési modellekre és mediációs modellekre van kihegyezve, mint amilyeneket itt végzünk, és erősen ajánlott az ilyen modellekhez. Bár a mediáció értékelésére a standard lineáris modell esetére korlátozódik, ez az egyetlen eszközünk, amely képes látens változókat készségesen beépíteni5. Például a depresszió kimenetelét az egyes kérdőívtételek alapjául szolgáló látens változónak tekinthetjük. Ezenkívül több mediátort és több kimenetelt is beépíthetnénk.

Azért, hogy a dolgok úgy maradjanak, ahogyan eddig tárgyaltuk, a lavaanban a a, b és c' ösvényeket aszerint jelölöm, ahogyan korábban ábrázoltuk őket. Egyébként a lavaan nagyon könnyen használható, és a megfigyelt változók esetében a modellekhez a szokásos R-formulák jelölését használja. Ezen túlmenően a := operátorral definiáljuk azokat az érdekes hatásokat, amelyeket ki akarunk számítani. A modellt teljes egészében egyszerű karakterlánc formájában adjuk meg, majd a sem függvénnyel elvégezzük az elemzést.

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

A korábbi kimenetet látjuk, és a indirect paraméterünket összehasonlíthatjuk a korábbi ACME-vel, a direct hatást az ADE-vel, a total pedig a korábbi összhatással. Az értékek lényegében megegyeznek.

Megjegyezzük azt is, hogy a kimenet mindkét modell esetében a \(R^2\) értéket mutatja. A job_seek esetében láthatjuk, hogy azért nem találunk sok közvetítést, mert az érintett kovariánsok eleve nem magyarázzák a közvetítő változásait. Az előzetes vizsgálat ebben az esetben megspórolta volna nekünk a fáradságot.

Előnyei

• Kezelhet több mediátort
• Kezelhet több “kezelést”
• Kezelhet több kimenetet
• Kezelhet látens változókat
• Némi többszintű támogatás
• Moderált mediációt és mediált közvetítést is végezhet. moderációt (bár látens változók esetén nem)

Korlátozások

• Kiegészítő kódolást igényel a közvetett hatás becsléséhez
• Egyetlen véletlen hatás
• Míg a modellekbe bináris vagy ordinális változókat lehet beépíteni a mediátor/eredmények esetében, nincs egyszerű módja annak, hogy a közvetett hatást a mediációs csomag módjára kiszámítsuk ezekben a beállításokban.

piecewiseSEM

A piecewiseSEM csomag nagyon hasonlóan működik, mint a mediációs csomag. Az a szép benne a mediációs csomaghoz képest, hogy a piecewiseSEM további típusú modelleket tud kezelni, valamint további kimenetet (pl. standardizált eredmények), további opciókat (pl. többcsoportos, korrelált reziduumok) és a modell vizualizálását biztosítja.

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

A modell gyors vizualizációjának elkészítéséhez használhatjuk az ábrázolási képességeit.

plot(mediation_result)

A közvetett hatások kiszámítására jelenleg sajnos nincs automatikus mód, így az eredményeket kézzel kell bootstrapolni.

Előnyei

• Szokásos R modellek és szintaxis
• Sokféle modelltípus mind a mediátor, mind a kimenetel esetében
• Néhány SEM-stílusú eredmény (pl. illeszkedés, standardizált együtthatók, AIC)
• Az eredmények gyors ábrázolása
• Kezelhet több mediátort, “kezelést”, és kimenetekkel

Korlátozások

• Nem számítja ki automatikusan a közvetett hatásokat
• Nincs látens változó képesség

psych

A psych csomag kihasználja, hogy a standard lineáris modell esetében az eredményeket a megfelelő regressziós modelleken keresztül csak a kovariancia mátrixok alapján kaphatjuk meg. Nagyon hasonlít a lavaanhoz, bár a maximum likelihood helyett közönséges legkisebb négyzetek megközelítést használ. A szép dolog itt az a szintaxis, amely lehetővé teszi, hogy csak az érdekes hatásra összpontosítson, vagy mindent bevonjon, ami jó, ha a gazdasági nehézségek, az életkor és a nem közvetett hatásai is érdekelnék.

Ezért a demóért a tisztított változatot fogjuk használni, amely a + helyett a --t használja a nem kezelési hatásokhoz. Ez csak azt jelenti, hogy szerepelnek a modellekben, de az eredmények nem jelennek meg rájuk vonatkozóan. A mediátort a ()-vel azonosítjuk. Egy másik bónusz az eredmények gyors ábrázolása, amely megmutatja a kiigazítatlan és a kiigazított közvetlen hatások közötti különbséget, valamint a megfelelő bootstrapped intervallumot.

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

Az eredmények ugyanazok, más csomagolásban, de talán az eddigi legegyszerűbb út, mivel csak egy függvényhívást igényelt. A pszichológiai csomag bónuszként több közvetítőt és eredményt is kezel.

Pros

• A legegyszerűbb szintaxis, alapvetően egysoros modell
• Gyorsan ábrázolja az eredményeket
• Kezelhet több mediátort, ‘kezelést’, és kimenetekkel
• Meg tudja csinálni a “moderált” mediációt

Korlátozások

• Korlátozott a standard lineáris modellre (lm)
• MASS használata

brms

A következő demónkban a szerintem legerősebb csomaghoz, a brms-hez érkezünk. A név a Bayesian Regression Modeling with Stan rövidítése, a Stan pedig egy hatékony valószínűségi programozási nyelv a Bayes-analízishez. Nem megyek bele a Bayes-analízis részleteibe, de nézze meg nyugodtan az erről szóló dokumentumomat.

Általában úgy járunk el, mint eddig, megadjuk a mediátor modellt és a kimeneti modellt. A brms nem csinál semmi különlegeset a mediációs elemzéshez, de a hipotézisfüggvénye lehetővé teheti számunkra, hogy teszteljük az útvonalak terméke megközelítést. Ezenkívül az sjstats csomag lényegében ugyanúgy szolgáltatja az eredményeket, ahogyan a mediációs csomag teszi ezt számunkra, és ha már itt tartunk, a mediációs csomag alapvetően amúgy is egy bayesi megoldási kísérlet a frequentista módszerekkel. Ha valóban különböző eloszlásokkal rendelkeznénk a kimenetelre és a közvetítőre vonatkozóan, akkor viszonylag könnyen megkaphatnánk ezeket az átlagos előrejelzési értékeket és azok különbségeit, mivel a bayesi megközelítések mindig utólagos előrejelzési eloszlásokban gondolkodnak. Mindenesetre itt van a kód.

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

A kimeneten minden, amiben jobseek_* szerepel, az a mediátor modell eredménye, míg depress2_* az eredményé. Ezen a ponton ugyanaz a régi történet áll előttünk, de a Bayes-megközelítéssel több szórakoztató dolgot nézhetünk meg. Például láthatjuk, hogy valójában nem jól ragadjuk meg a depresszió kimenetelének ferdeségét. A megjósolt értékeink és a megfigyelt értékek nem teljesen egyeznek. A mediátor esetében egy kicsit jobban állunk, de talán még mindig egy kicsit magasan vagyunk néhány modellalapú előrejelzésünkkel.

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

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

Pros

• egyszerű szintaxis
• Extrem erős- A modellek többnyire a képzeletre korlátozódnak
• Gyakorlatilag azt teszi, amit a mediációs csomag közelít
• A bayesi következtetés összes előnye: Diagnosztika, utólagos előrejelző ellenőrzések, modell-összehasonlítás stb.

Korlátozások

• Kisebb becslés
• ‘Kézzel’ végzett számítások szükségesek a standard lineáris modellen való túllépéshez, de ez már általános megközelítés a bayesi szemszögből
• Némi kényelem szükséges a bayesi megközelítéssel kapcsolatban

Nagyobb komplexitás

A felsorolt csomagok némelyike képes bonyolultabb modelleket kezelni, vagy további megközelítéseket nyújt a közvetett hatások vizsgálatához.

Interakciók

Egyes modellek interakciókat tartalmaznak akár a közvetítő modell, akár a kimenet tekintetében, és ezt sajnos gyakran mediált moderációnak vagy moderált közvetítésnek nevezik. Én személy szerint nem látom az előnyét annak, hogy kétértelmű neveket adjunk egy egyébként egyszerű koncepciónak (ha még mindig nem is olyan egyértelmű modellnek), de ez a hajó már régen elment. Nem fogok belemenni a részletekbe, de az ötlet az, hogy valahol a modellben interakciós kifejezés lehet, és a kölcsönhatás magában foglalhatja a kezelési változót, a mediátort vagy mindkettőt.

Legyen elég annyi, hogy mivel olyan standard modellezési eszközöket használunk, mint a lm és annak kiterjesztései, a kölcsönhatások beépítése triviális az összes fenti csomag esetében, de a product-of-paths típusú megközelítés nem érvényesül (a*b != c').

Generalizált lineáris modellek

Egyes esetekben a mediátorunk vagy kimenetünk lehet bináris, számláló, vagy valami, ahol a normál eloszlás feltételezése nem feltétlenül a legjobb ötlet. Vagy esetleg nemlineáris kapcsolatokat szeretnénk vizsgálni a kezelés/közvetítő/eredmény között. Vagy lehetnek olyan adataink, amelyek korrelált megfigyeléseket tartalmaznak, mint például ismételt mérések vagy hasonló. A mediációs csomag különösen büszke erre, de a brms bármit meg tud csinálni, amit csak tud, és még többet is, bár lehet, hogy egy kicsit több munkát kell végeznünk, hogy ténylegesen kiszámítsuk az eredményt. lavaan valóban képes egy korlátozott számú modellt elvégezni bináris és ordinális változókra, de a megfelelő közvetett becslés megszerzése nagyon fárasztó kézi megközelítést igényelne.

Missing data

Az ilyen adatok kezelése során, különösen a társadalomtudományokban, gyakran hiányoznak adatok bármelyik kovariánsról. Néha ezeket elhagyhatjuk, ha nincs túl sok, de más esetekben valamit tenni akarunk vele. A lavaan, psych és brms csomagok egy vagy több lehetőséget biztosítanak a helyzet kezelésére (pl. többszörös imputálás).

Alternatívák

A modelleket eddig csomópontok hálózataként ábrázoltuk, amelyeket ívek/élek/útvonalak kötnek össze. A vitánk az úgynevezett irányított aciklikus gráfok (DAG) körül forog, ahol a nyilak csak egy irányba mehetnek, visszacsatolási hurok nélkül. Bármely kimeneti változó eredménye az azt megelőző nyilak függvénye, és feltételesen független a többitől. Egyes elméleti modellek ezt lazíthatják, mások pedig egyáltalán nem tartalmaznak nyilakat, azaz irányítatlanok, így minket csak a kapcsolatok érdekelnek (pl. egyes társadalmi hálózatoknál).

bnlearn

A bnlearn csomag lehetővé teszi irányított, részben irányított és irányítatlan gráfok vizsgálatát. A DAG-ok szempontjából lényegében az eddig tárgyalt közvetítési modellek megkettőzésére használhatjuk. A jó dolog azonban az, hogy ez a csomag hatékonyan teszteli az útvonalakat a felvételre, ahelyett, hogy feltételezné őket, de szükség szerint elméleti megkötéseket továbbra is előírhatunk. Ezután nem csak a Bayes-hálózatok és Pearl oksági gráfelméletének alapjául szolgáló Bayes-hálózatokkal és Pearl oksági gráfelméletével elvi módon kereshetjük meg az érdekes útvonalakat, hanem olyan eszközökkel is rendelkezni fogunk, amelyekkel keresztvalidálással tovább kerülhetjük a túlillesztést.

A kezdeti modell esetében biztosítjuk, hogy a kezelés – mediátor, kezelés – eredmény és mediátor – eredmény (a fehér lista) között léteznek útvonalak. Az olyan értelmetlen útvonalakat, mint a kezelésre (amelyet véletlenszerűen osztottunk ki), a nemre, a gazdasági nehézségekre és az életkorra mutató nyilak (a feketelista), letiltjuk. Egyébként megnézzük, mit sugallnak az adatok.

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)

A grafikonon látjuk, hogy a dolgok egy kicsit megváltoztak. Például az életkor most már csak az álláskeresési önhatékonysággal van kapcsolatban, a nem pedig csak a depresszióra van hatással.

Ha az útvonalakat úgy korlátozzuk, hogy csak azok legyenek, mint a korábbi példáinkban, akkor ugyanazokat az eredményeket kapjuk.

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 

A legfontosabb, hogy a becsült paraméterek megegyeznek azzal, amit a korábbi csomagokkal kaptunk. Ez lényegében egyenértékű a lavaan használatával az alapértelmezett maximum likelihood becslővel.

Ha a kezelést és a nemet használjuk faktorokként, a bnlearn olyan feltételes modelleket fog produkálni, amelyek különbözőek attól függően, hogy milyen faktorértéket veszünk. Más szóval, külön modellünk lenne arra az esetre, ha treatment == 'treatment', és külön modellünk arra az esetre, ha treatment == control. A mi esetünkben ez azonos lenne azzal, ha mindent kölcsönhatásba engednénk a kezeléssel, pl. lm( job_seek ~ treat * (econ_hard + sex + age)), és ugyanígy a depressziós modell esetében is. Ez kiterjeszthető lenne potenciálisan bármilyen bináris változóra (pl. a nemet is beleértve). Ha a mediátor bináris változó, akkor valószínűleg ezt szeretnénk csinálni.

Python

CSCAR igazgatója, Kerby Shedden tartott egy Python workshopot a mediációs modellekről, ezért itt mutatom be a statsmodels implementációt. Ez az Imai megközelítést követi, így a mediációs csomag Python verziójának tekinthető. A kimenet lényegében ugyanaz, mint amit a kezelést faktorváltozóként használva kapnánk, ahol minden egyes kezelési kategóriára külön eredményeket kapunk. Erre a mi demónkban nincs szükség, így csak az “átlagos” eredményeket lehet összehasonlítani az előző mediációs csomag eredményeivel.

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

Végül a Stata sem parancsával biztosítok egy lehetőséget. A Stata megkönnyíti a közvetett hatások kinyerését ebben a példában, de ezt minden kovariánsra vonatkozóan megteszi, így a kimenet enyhén szólva is kissé bőbeszédű6. A Statával dolgozóknak nincs szükségük külön SEM csomagra ahhoz, hogy ilyen jellegű eredményeket kapjanak.

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)

Összefoglaló

A közvetett hatásokkal rendelkező modellek gondos elméleti megfontolást igényelnek az adatelemzéshez való alkalmazáshoz. Ha azonban a modell megfelel az adathelyzetnek, akkor az R különböző csomagjaival meglehetősen könnyű eredményeket kapni. Továbbá nem kell strukturális egyenletmodellező csomagot használni a közvetett hatásokkal végzett elemzéshez, sőt, a standard R szintaxis használatával is messzire juthatunk. Szigorúan megfigyelt, azaz nem látens változók esetén nincs szükség SEM eszközre, sőt nem is ajánlott.

\

Package comparison summarized

A következő táblázat segíthet eldönteni, hogy elméleti megfontolások alapján melyik csomagot használjuk az igényeinkhez.

1. A SEM-ről van egy sokkal részletesebb dokumentumom, beleértve a mediációs elemzést is.︎

2. Valamiért a gyakorlatban nem nagyon látni ilyet, és az ember elgondolkodik, hogy mit tettek azért, hogy az adatokat alkalmassá tegyék egy ilyen modellre, ha az nem volt indokolt.︎

3. Imai elérhetővé teszi cikkeit a honlapján.︎

4. A MASS-t már több mint egy évtizede felváltották mások, és többnyire csak arra hajlamos, hogy összezavarja a tidyverse-t és más csomagokat, amikor betöltik. Ez egy jó csomag (és nagyszerű volt annak idején), de ha egy csomagban akarod használni, akkor jó lenne, ha nem töltenéd be (vagy más csomagokat) a környezetbe csak azért, hogy egy-két függvényt használj. Leginkább csak az mvrnorm (többváltozós normális eloszlás) és a glm.nb esetében látom használni, de vannak más csomagok is ezzel a funkcionalitással, amelyek további előnyöket biztosítanának, és nem a dplyr függvényeket takarják, amelyek az R közösségben a leggyakrabban használtak közé tartoznak.︎

5. brms dolgozik rajta.︎

6. A kódban lévő opciók azért vannak, hogy elnyomják/csökkentik, ami lehet.︎