Free Press

Sisällysluettelo

Päivitetty 04. elokuuta 2020. Koodi on ladattavissa täältä.

Esittely

Jossain tilanteissa saatamme tarkastella jonkin muuttujan epäsuoraa vaikutusta lopputulokseen tai tulokseen. Esimerkkinä voidaan mainita, että huonot elinolot kotona lapsuudessa voivat heikentää oppimistuloksia koulussa, mikä myöhemmin vaikuttaa kielteisesti myöhempään elämänlaatuun, esimerkiksi elinikäisiin tulotuloihin. Toisessa tapauksessa voidaan tarkastella yhtä muuttujaa, joka on kerätty useampana ajankohtana, jolloin muuttujalla on vaikutus ajankohtana 1 ajankohtana 2 ja ajankohtana 2 ajankohtana 3. Tämän muuttujan vaikutus voi olla myös suurempi. Perusajatus on jotakuinkin seuraava:

\

Toisin sanoen \(\mathcal{A}\) johtaa \(\mathcal{B}\), ja sitten \(\mathcal{B}\) johtaa \(\mathcal{C}\). Välitysmallien avulla asetamme normaalin kovariaatin \(\rightarrow\) tulospolun väliin välivaiheen muuttujan, joka meillä voisi olla tavallisessa regressioasetelmassa, ja näiden mallien avulla voimme tutkia tällaista käyttäytymistä. Edellä esitetyssä tapauksessa väliintuleva muuttuja tai välittäjä on \(\mathcal{B}\). Usein on mahdollista, että \(\mathcal{A}\) vaikuttaa edelleen suoraan \(\mathcal{C}\), mutta kuten mallissa yleensä, tämä olisi teoreettisesti perusteltua.

Välittäjäanalyysi on hyvin suosittua yhteiskuntatieteellisillä tieteenaloilla, vaikkakaan se ei suinkaan rajoitu niihin, ja sitä tehdään yleensä rakenteellisen yhtälömallinnuksen (SEM) varjolla, joka itsessäänkin on yleisemminkin grafiikkamallien erityissuuntautuneisuus1. Välitysmallin graafinen malli voi näyttää seuraavalta.

Tässä tapauksessa a ja b kuvastavat välittäjän kautta kulkevaa välillistä polkua \(\mathrm{X}\) vaikutuksesta lopputulokseen, kun taas c' on \(\mathrm{X}\) suora vaikutus lopputulokseen sen jälkeen, kun epäsuora polku on poistettu (c olisi vaikutus ennen epäsuoran vaikutuksen asettamista, ja c – c' vastaa epäsuoraa vaikutusta). Kokonaisvaikutus \(\mathrm{X}\) on yhdistetty epäsuora ja suora vaikutus.

Huomautan muutaman asian sen perusteella, mitä näen konsultoidessani kymmeniä tieteenaloja. Aluksi näyttää siltä, että hyvin harvat ihmiset, jotka luulevat tarvitsevansa välitysmallia, todella tarvitsevat sitä. Jos et esimerkiksi pysty ajattelemaan malliasi ajallisesti tai fysikaalisesti siten, että \(\mathrm{X}\) johtaa välttämättä välittäjään, joka sitten johtaa välttämättä lopputulokseen, et todennäköisesti tarvitse välitysmallia. Jos voisit nähdä nuolien kulkevan kumpaankin suuntaan, et luultavasti myöskään tarvitse tällaista mallia. Myös jos kaikki luulevat malliasi kuvattaessa, että puhut vuorovaikutuksesta (eli moderoinnista), et ehkä tarvitse tätä. Ja lopuksi, kuten voisi epäillä, jos keskeisten muuttujien (\(\mathrm{X}\)) ja välittäjän (polku a) välillä ei ole vahvaa korrelaatiota ja jos välittäjän ja lopputuloksen (lopputulosten) välillä ei ole vahvaa korrelaatiota (polku b), et luultavasti tarvitse tätä. Vaikka mikään ei estä sinua tekemästä mediaatioanalyysia, ilman näitä edellytyksiä saat lähes varmasti heikon ja luultavasti sekavamman mallin kuin mitä muuten saisit.

Lyhyesti sanottuna, mediaatio toimii parhaiten, kun muuttujien välillä on vahvasti implikoituja kausaalisia yhteyksiä. Silloinkin tällaista mallia tulisi verrata yksinkertaisempaan malliin, jossa ei ole välitystä2. Joka tapauksessa on olemassa muutamia hyvin helppoja tapoja tutkia tällaisia malleja R:ssä, ja se on tässä tavoitteena, vain demonstroida, miten voit päästä alkuun.

Data

Demonstraatioksi mediatiomalleista eri paketeilla käytämme jobs-dataa, joka tulee mediation-paketin mukana. Tässä on kuvaus.

Job Search Intervention Study (JOBS II). JOBS II on satunnaistettu kenttäkoe, jossa tutkitaan työhönvalmennusintervention tehoa työttömiin työntekijöihin. Ohjelman tarkoituksena on paitsi lisätä työttömien uudelleen työllistymistä myös parantaa työnhakijoiden mielenterveyttä. JOBS II -kenttäkokeilussa 1 801 työtöntä työntekijää sai esitutkimuskyselyn, minkä jälkeen heidät jaettiin satunnaisesti hoito- ja kontrolliryhmiin. Hoitoryhmään kuuluvat osallistuivat työelämätaitoja käsitteleviin työpajoihin. Työpajoissa vastaajat oppivat työnhakutaitoja ja selviytymisstrategioita työnhakuprosessin takaiskujen käsittelemiseksi. Kontrolliryhmään kuuluvat saivat kirjasen, jossa kuvattiin työnhakuvinkkejä. Seurantahaastatteluissa kaksi keskeistä tulosmuuttujaa olivat Hopkinsin oireiden tarkistuslistaan perustuva jatkuva masennusoireiden mittari ja binäärimuuttuja, joka kuvasi sitä, oliko vastaaja työllistynyt.

Tässä on kuvaus tämän demonstraation muuttujista. Käytettävissä on muitakin, joilla voit myös leikkiä.

• econ_hard: Taloudellisten vaikeuksien taso ennen käsittelyä, arvot 1-5.
• sex: Indikaattorimuuttuja sukupuolelle. 1 = nainen
• age: Ikä vuosina.
• educ: Koulutustasoa kuvaava tekijä, jossa on viisi luokkaa.
• job_seek: Jatkuva asteikko, joka mittaa työnhaun itsetehokkuuden tasoa ja jonka arvot ovat 1-5. Välittävä muuttuja.
• depress2: Masennusoireiden mittaaminen hoidon jälkeen. Tulosmuuttuja.
• treat: Indikaattorimuuttuja sille, valittiinko osallistuja satunnaisesti JOBS II -koulutusohjelmaan. 1 = osoittaminen osallistumiseen.

data(jobs, package = 'mediation')

Malli

Tämän aineiston perusteella välittäjä- ja lopputulosmallit ovat seuraavat:

\

Edellytys on siis, että työelämätaitojen harjoittelulla olisi negatiivinen vaikutus masennukseen (eli hyvinvoinnin lisääntyminen), mutta ainakin osa tästä johtuisi positiivisesta vaikutuksesta työnhakuun.

Graafisena mallina voisimme kuvata asian ytimekkäästi seuraavasti.

Paketit

Tarkastelemme seuraavia paketteja havainnollistaaksemme, miten välitysanalyysia voidaan tehdä R:ssä:

• mediation
• lavaan
• psych
• brms

Näihin keskitytään, mutta huomioin myös joitakin muita vaihtoehtoja, kuten Python ja Stata.

mediation

Aloitamme mediation-paketilla, koska sen suorittaminen ei periaatteessa vaadi enempää ohjelmointitaitoja kuin mitä omistaa jo tavallisten regressiomallien suorittamisesta R:ssä. Paketti tarjoaa keskimääräisen kausaalisen välitysvaikutuksen (average causal mediation effect), joka määritellään seuraavasti ohjetiedostosta ja Imain artikkeleista3:

Keskimääräinen kausaalinen välitysvaikutus (average causal mediation effect, ACME) edustaa odotettua eroa potentiaalisessa lopputuloksessa silloin, kun välittäjä ottaisi arvon, joka toteutuisi hoitotilanteen vallitessa vastakohtana kontrollointitilanteelle, kun taas hoitotilanne sinällään pidetään muuttumattomana.

>

Huomaa, kuinka tässä määritelmässä keskitytään odotettuihin tai ennustettuihin arvoihin, jotka ovat riippuvaisia hoitotilan arvosta. Tällä käsitteellä kontrafaktuaaleista eli siitä, miltä havainto näyttäisi päinvastaisessa asetelmassa, on tässä vaiheessa pitkä historia mallintamisessa. Ajattele asiaa näin: jos joku on hoitoryhmässä, hänellä olisi tietty arvo välittäjälle, ja tämän perusteella hänellä olisi tietty odotettu arvo lopputulokselle. Voisimme kuitenkin olettaa, että sama havainto on myös kontrolliryhmässä, ja arvioida välittäjän vaikutusta lopputulokseen aivan samalla tavalla. Voimme arvioida mahdollisia tuloksia pitäen hoidon vakiona. Tulosmuutosten tarkastelu välittäjän arvon perusteella ei tee oletusta mallityypistä. Näin välityspaketti pystyy sisällyttämään erilaisia malleja välittäjälle ja lopputulokselle. Välittäjä voi esimerkiksi olla binäärinen, jolloin tarvitaan logistinen regressiomalli, kun taas lopputulosmalli voi olla eloonjäämismalli.

Esimerkissämme pysytään tavanomaisissa (normaaleissa) lineaarisissa malleissa. Huomaa myös, että vaikka käsittelymme on binäärinen muuttuja, tämä yleistyy jatkuvaan tapaukseen, jossa tarkastelemme tulosta, joka syntyy yhden yksikön liikkeestä ”käsittelyssä”. Jotta välityspaketti toimisi, suoritamme yksinkertaisesti mallit välittäjälle ja lopputulokselle ja käytämme sitten mediate-funktiota saadaksemme lopputuloksen.

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)

Ylläolevat tulokset osoittavat, että ACME ei eroa tilastollisesti nollasta tai välittämisestä. Keskimääräinen suora vaikutus on negatiivinen, mutta sekään ei ole tilastollisesti merkittävä, kuten ei myöskään kokonaisvaikutus (epäsuora + suora vaikutus). Lisäksi on annettu soi disant ”välitetyn vaikutuksen osuus”, joka on epäsuoran vaikutuksen suhde kokonaisvaikutukseen. Tämä ei kuitenkaan ole suhdeluku, ja se voi olla jopa negatiivinen, joten se on enimmäkseen merkityksetön luku.

Pros

• Standardit R-mallit ja syntaksi
• Monentyyppisiä malleja sekä välittäjälle että lopputulokselle
• Tarjoaa useita tuloksia yhtäaikaisesti
• Hyvä dokumentaatio ja siihen liittyvät artikkelit ovat vapaasti saatavissa
• Voi tehdä ’moderoitua’ välittämistä

Limiitit

• MASS4:n käyttö
• Yksinkertaiset satunnaisefektimallit
• Toiminnallisuus ehkä rajoitettu joidenkin mallien monimutkaisuuden vuoksi
• Ei latenttimuuttujaominaisuuksia

lavaan

Erikoistapauksessa, jossa sekä välitystapa- että lopputulosmallit ovat tavanomaisia lineaarisia malleja, joissa on normaali jakauma kohdemuuttujan osalta, epäsuora vaikutus vastaa edellisen kaavion a– ja b-polkujen tuloa. Suora vaikutus on polku c'. Vertailu itsenäisen suoran vaikutuksen, jota voisimme kutsua c, ja tämän estimoidun suoran vaikutuksen välillä välitysmallissa c', on sellainen, että c - c' = a*b. Se, mitä aiemmin mainittiin, saattaa nyt olla selkeämpää, jos joko a tai b on lähes nolla, niin epäsuora vaikutus voi olla vain lähes nolla, joten on järkevää tutkia tällaisia suhteita etukäteen.

Tämä polkujen tuote (tai kertoimien erotus) -lähestymistapa on se, jota tulemme käyttämään lavaan-paketin avulla, ja itse asiassa tätä kirjoittaessamme se on ainoa tapamme toimia. lavaan on suunnattu erityisesti rakenteelliseen yhtälömallinnukseen, kuten faktorianalyyseihin, kasvumalleihin ja välitysmalleihin, jollaisia teemme tässä, ja sitä suositellaan lämpimästi tällaisia malleja varten. Vaikka se rajoittuu tavanomaiseen lineaarisen mallin tapaukseen arvioida välitystä, se on ainoa työkaluistamme, joka voi sisällyttää latentteja muuttujia helposti5. Esimerkiksi masennustuloksemme voisi olla yksittäisten kyselylomakkeen kohtien taustalla oleva latentti muuttuja. Lisäksi voisimme myös sisällyttää useita välittäjiä ja useita lopputuloksia.

Pitäytyäkseni siinä, mistä olemme keskustelleet, merkitsen lavaan a, b ja c' polut sen mukaan, miten ne on aiemmin kuvattu. Muuten lavaan on erittäin helppokäyttöinen, ja havainnoitujen muuttujien tapauksessa se käyttää malleissa tavallista R-kaavojen merkintätapaa. Sen lisäksi määrittelemme kiinnostavat vaikutukset, jotka haluamme laskea :=-operaattorilla. Määrittelemme mallin kokonaisuudessaan yksinkertaisena merkkijonona ja käytämme sitten sem-funktiota analyysin tekemiseen.

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

Näemme saman tulosteen kuin ennenkin ja voimme verrata indirect-parametria ACME:hen, joka meillä oli ennenkin, direct-vaikutusta verrataan ADE:hen ja total verrataan edelliseen kokonaisvaikutukseen. Arvot ovat pääosin samat.

Huomaa myös, että tulosteessa näkyy \(R^2\) arvo molemmille malleille. job_seek:n tapauksessa voimme nähdä, että syy siihen, että emme löydä juuri mitään sovittelua, on se, että mukana olevat kovariaatit eivät alun perin selitä mitään vaihtelua sovittelijassa. Ennakkotutkimus olisi tässä tapauksessa säästänyt meidät vaivalta.

Pros

• Voi käsitellä useita välittäjiä
• Voi käsitellä useita ”hoitoja”
• Voi käsitellä useita lopputuloksia
• Voi käyttää latentteja muuttujia
• Joitakin monitasoisia tukia
• Voi tehdä moderoitua välittämistä ja medioituja moderoinnin (ei kuitenkaan latentteja muuttujia)

Limitations

• Tarvitsee lisäkoodausta epäsuoran vaikutuksen estimoimiseksi
• Yksittäisiä satunnaisvaikutuksia
• Malleihin voisi sisällyttää binääri- tai ordinaalimuuttujia välittäjälle/tuloksille, ei ole suoraviivaista tapaa laskea välillistä vaikutusta välityspaketin tapaan näissä asetuksissa.

piecewiseSEM

P piecewiseSEM-paketti toimii hyvin samankaltaisesti kuin mediation-paketti. Hienoa tässä suhteessa mediation-pakettiin on se, että piecewiseSEM voi käsitellä lisätyyppisiä malleja sekä tarjota lisätulosteita (esim. standardoidut tulokset), lisäoptioita (esim. moniryhmä, korreloidut residuaalit) ja mallin visualisointia.

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

Voidaan käyttää sen piirto-ominaisuuksia mallin nopean visualisoinnin luomiseen.

plot(mediation_result)

Epäsuorien vaikutusten laskemiseen ei valitettavasti ole tällä hetkellä automaattista tapaa, joten tulokset pitäisi bootstrapata käsin.

Pros

• Standardit R-mallit ja syntaksi
• Monentyyppisiä malleja sekä välittäjälle että lopputulokselle
• Joitakin SEM-tyylisiä tuloksia (esim. fit, standardoidut kertoimet, AIC)
• Tulosten nopea piirtäminen
• Voi käsitellä useita välittäjiä, ’treatments’, and outcomes

Limitations

• Ei laske automaattisesti epäsuoria vaikutuksia
• Ei latenttimuuttujaominaisuuksia

psych

Psykopaketti hyödyntää sitä, että tavanomaisen lineaarisen mallin tapauksessa tulokset saadaan sopivien regressiomallien kautta pelkkien kovarianssimatriisien perusteella. Se on hyvin samankaltainen kuin lavaan, vaikkakin käyttää tavallista pienimmän neliösumman lähestymistapaa maksimiluotettavuuden sijasta. Hienoa tässä on syntaksi, jonka avulla voit keskittyä vain kiinnostavaan vaikutukseen tai sisällyttää kaiken, mikä on mukavaa, jos olisit kiinnostunut myös taloudellisten vaikeuksien, iän ja sukupuolen epäsuorista vaikutuksista.

Tässä demossa käytämme puhdistettua versiota, jossa käytetään --merkintää +:n sijaan muille kuin hoitovaikutuksille. Tämä tarkoittaa vain sitä, että ne sisältyvät malleihin, mutta niitä koskevia tuloksia ei näytetä. Välittäjä tunnistetaan ():llä. Toinen bonus on tulosten nopea kuvaaja, jossa näkyy sopeuttamattomien ja sopeutettujen suorien vaikutusten välinen ero ja sopiva bootstrapped-väli.

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

Samoja tuloksia, erilainen pakkaus, mutta kenties helpoin reitti tähän asti, koska se vaati vain yhden funktiokutsun. Psykopaketti käsittelee bonuksena myös useita välittäjiä ja tuloksia.

Pros

• Leppoisin syntaksi, periaatteessa yhden rivin malli
• Nopea tulosten piirtäminen
• Voi käsitellä useita välittäjiä, ’hoitoja’, ja tuloksia
• Voi tehdä ’moderoitua’ sovittelua

Limitations

• Limitoituu tavalliseen lineaariseen malliin (lm)
• Use of MASS

brms

Seuraavaksi esittelyssämme käsittelemme mielestäni mahtavinta pakettia brms. Nimi tulee sanoista Bayesian Regression Modeling with Stan, ja Stan on tehokas todennäköisyysohjelmointikieli bayesilaiseen analyysiin. En aio mennä yksityiskohtiin Bayesian analyysistä, mutta tutustukaa vapaasti dokumenttiini, joka tekee niin.

Toimimme yleensä kuten aiemminkin määrittelemällä välittäjämallin ja lopputulosmallin. brms ei tee mitään erityistä välitysanalyysiä varten, mutta sen hypoteesifunktio voi antaa meille mahdollisuuden testata product-of-paths -lähestymistapaa. Lisäksi sjstats-paketti antaa tulokset periaatteessa samalla tavalla kuin mediaatiopaketti, ja itse asiassa mediaatiopaketti on pohjimmiltaan joka tapauksessa yritys Bayesin ratkaisuun frequentistisia menetelmiä käyttäen. Jos meillä olisi erilaiset jakaumat lopputulokselle ja välittäjälle, meidän olisi suhteellisen helppo saada nämä keskimääräiset ennustearvot ja niiden erot, koska bayesiläiset lähestymistavat ajattelevat aina posteriorisia ennustejakaumia. Joka tapauksessa tässä on koodi.

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

Tuloksessa kaikki, missä on jobseek_*, on välittäjämallin tulos, kun taas depress2_* on tuloksen tulos. Tässä vaiheessa meillä on sama vanha tarina, mutta Bayesin lähestymistavan avulla meillä on hauskempia asioita tarkasteltavana. Voimme esimerkiksi nähdä, että emme itse asiassa kuvaa hyvin masennustuloksen vinoutta. Ennustetut arvot ja havaitut arvot eivät aivan vastaa toisiaan. Olemme hieman parempia välittäjän osalta, mutta ehkä edelleen hieman korkealla joidenkin mallipohjaisten ennusteidemme kanssa.

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

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

Pros

• Suoraviivainen syntaksi
• Erittäin tehokas- Mallit rajoittuvat enimmäkseen omaan mielikuvitukseen
• Periaatteessa tekee sen, mitä välityspaketti approksimoi
• Kaikki Bayesin päättelyn edut: Diagnostiikka, posterioriset ennustetarkastelut, mallien vertailu jne.

Limitations

• Halvempi estimoida
• ’By-hand’-laskelmia tarvitaan, kun mennään pidemmälle kuin tavallinen lineaarinen malli, mutta tämä on jo yleinen lähestymistapa bayesiläisestä näkökulmasta
• Tarvitaan jonkin verran mukavuutta bayesiläisen lähestymistavan kanssa

Lisää monimutkaisuutta

Jotkut mainituista paketeista pystyvät käsittelemään monimutkaisempia malleja tai tarjoavat lisälähestymistapoja epäsuorien vaikutusten tutkimiseen.

Vuorovaikutukset

Joissakin malleissa on vuorovaikutuksia joko välitysmallin tai lopputuloksen osalta, ja valitettavasti tähän viitataan usein nimellä mediated moderation tai moderated mediation. Henkilökohtaisesti en näe etua siitä, että annetaan moniselitteisiä nimiä sille, mikä muuten voisi olla suoraviivainen käsite (joskaan ei silti niin suoraviivainen malli), mutta se laiva purjehti jo kauan sitten. En aio mennä yksityiskohtiin, mutta ajatuksena on, että mallissa saattaa olla vuorovaikutustermi jossain, ja vuorovaikutukseen saattaa liittyä hoitomuuttuja, välittäjä tai molemmat.

Se riittää, että koska käytämme vakiomallinnustyökaluja, kuten lm ja sen laajennuksia, vuorovaikutusten sisällyttäminen on triviaalia kaikille edellä mainituille paketeille, mutta polkujen tuote -tyyppinen lähestymistapa ei päde (a*b != c').

Yleistetyt lineaariset mallit

Jossain tapauksissa mediaattorimme tai lopputuloksemme voi olla binäärinen, laskentatulos tai jotain muuta sellaista, jossa normaalin jakauman olettamana olettamukset eivät ehkä olekaan paras mahdollinen idea. Tai saatamme haluta tutkia epälineaarisia suhteita hoidon/välittäjän/tuloksen välillä. Tai meillä voi olla tietoja, joissa on korreloituneita havaintoja, kuten toistuvia mittauksia tai vastaavaa. Mediation-paketti on ylpeä erityisesti tästä, mutta brms voi tehdä mitä tahansa, mitä se voi tehdä, ja enemmänkin, vaikka saatat joutua tekemään hieman enemmän työtä saadaksesi tuloksen oikeasti laskettua. lavaan voi itse asiassa tehdä rajoitetun joukon malleja binääri- ja ordinaalimuuttujille, mutta sopivan epäsuoran estimaatin saaminen vaatisi hyvin työlästä käsityötä.

Puuttuvat tiedot

Usein tällaisten aineistojen kanssa työskennellessämme erityisesti yhteiskuntatieteiden alalla tiedot puuttuvat usein jostakin kovariaatista. Joskus voimme jättää nämä pois, jos niitä ei ole liikaa, mutta toisissa tapauksissa haluamme tehdä asialle jotain. Paketit lavaan, psych ja brms tarjoavat yhden tai useamman tavan käsitellä tilannetta (esim. moninkertainen imputaatio).

Vaihtoehdot

Olemme kuvanneet mallit solmujen verkostoina, joita yhdistävät kaaret/särmät/polut. Keskustelumme pyörii niin sanottujen suunnattujen asyklisten graafien (Directed Acyclic Graphs, DAG) ympärillä, joissa nuolet voivat kulkea vain yhteen suuntaan ilman palautesilmukoita. Minkä tahansa tulosmuuttujan tulos on sitä edeltävien nuolten funktio ja ehdollisesti riippumaton muista. Joissakin teoreettisissa malleissa tätä voidaan höllentää, ja toisissa malleissa ei ole lainkaan nuolia, eli ne ovat suuntaamattomia, jolloin olemme kiinnostuneita vain yhteyksistä (esim. joidenkin sosiaalisten verkostojen kohdalla).

bnlearn

Bnlearn-paketti mahdollistaa suunnattujen, osittain suunnattujen ja suuntaamattomien graafien tutkimisen. DAG:ien osalta voimme käyttää sitä lähinnä monistamaan välitysmalleja, joista olemme keskustelleet. Hienoa on kuitenkin se, että tämä paketti testaa tehokkaasti polkujen sisällyttämistä sen sijaan, että olettaisi niitä, mutta voimme silti asettaa teoreettisia rajoituksia tarpeen mukaan. Sen lisäksi, että voimme sitten etsiä kiinnostavia polkuja periaatteellisella tavalla bayesilaisten verkkojen ja Pearlin kausaaligraafiteorian pohjalta, meillä on myös työkaluja, joilla voimme edelleen välttää ylisovittamista ristiinvalidoinnin avulla.

Alustavan mallin osalta varmistamme, että polkuja on olemassa hoidon – välittäjän, hoidon – lopputuloksen ja välittäjän – lopputuloksen välillä (valkoinen lista). Kiellämme järjettömät polut, kuten nuolet hoitoon (joka osoitettiin satunnaisesti), sukupuoleen, taloudellisiin vaikeuksiin ja ikään (musta lista). Muuten katsomme, mitä aineisto antaa ymmärtää.

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)

Havainnosta näemme, että asiat ovat hieman muuttuneet. Esimerkiksi ikä liittyy nyt vain työnhaun itsetehokkuuteen, ja sukupuolella on vaikutus vain masennukseen.

Jos rajoitamme polut vain sellaisiksi kuin ne ovat aiemmissa esimerkeissämme, saamme samat tulokset.

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 

Huomattava on se, että estimoidut muuttujat ovat yhtä suuret kuin aiemmilla paketeillamme. Se vastaa olennaisesti lavaanin käyttöä oletusarvoisella Maximum Likelihood -estimaattorilla.

Jos käytämme tekijöinä hoitoa ja sukupuolta, bnlearn tuottaa ehdollisia malleja, jotka ovat erilaisia riippuen otetusta tekijäarvosta. Toisin sanoen saataisiin erillinen malli sille, kun treatment == 'treatment', ja toinen sille, kun treatment == control. Meidän tapauksessamme tämä olisi identtistä sen kanssa, että annettaisiin kaiken olla vuorovaikutuksessa hoidon kanssa, esim. lm( job_seek ~ treat * (econ_hard + sex + age)), ja samoin masennuksen mallissa. Tämä ulottuisi mahdollisesti mihin tahansa binääriseen muuttujaan (esim. sukupuoli mukaan lukien). Jos välittäjä on binäärimuuttuja, tämä on todennäköisesti se, mitä haluaisimme tehdä.

Python

CSCAR:n johtaja Kerby Shedden on pitänyt Python-työpajan välitysmalleista, joten esittelen tässä statsmodels-toteutuksen. Se noudattaa Imai:n lähestymistapaa, joten sitä voidaan pitää Python-versiona välityspaketista. Tulos on periaatteessa sama kuin se, joka saataisiin käyttämällä hoitoa tekijämuuttujana, jolloin saat erilliset tulokset kullekin hoitoluokalle. Tämä on tarpeetonta demomme kannalta, joten voit vain verrata ”keskimääräisiä” tuloksia edellisen mediatiopaketin tuloksiin.

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

Viimeiseksi annan Statassa vaihtoehdon sen sem-komennon avulla. Stata helpottaa epäsuorien vaikutusten saamista tässä esimerkissä, mutta se tekee sen jokaiselle kovariaatille, joten tuloste on lievästi sanottuna hieman runsas6. Statan kanssa työskentelevät eivät tarvitse erillistä SEM-pakettia tällaisten tulosten saamiseksi.

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)

Yhteenveto

Mallit, joissa on epäsuoria vaikutuksia, vaativat huolellista teoreettista pohdintaa, jotta ne voidaan ottaa käyttöön data-analyysissä. Jos malli on kuitenkin sopiva aineistotilanteeseen, on melko helppoa saada tuloksia erilaisista R-paketeista. Lisäksi ei tarvitse käyttää rakenneyhtälömallinnuspakettia suorittaakseen analyysin epäsuorilla vaikutuksilla, ja itse asiassa voi päästä pitkälle käyttämällä tavallista R-syntaksia. Puhtaasti havaittujen, eli ei latenttien, muuttujien osalta ei SEM-työkalua tarvita tai edes suositella.

\

Pakettien vertailu tiivistettynä

Seuraava taulukko voi auttaa päättämään, mitä pakettia käyttää tarpeisiinsa ottaen huomioon teoreettiset näkökohdat.

1. Minulla on paljon yksityiskohtaisempi dokumentti SEM:stä, mukaan luettuna välitysanalyysi.︎

2. Jostain syystä tätä ei juurikaan näe käytännössä, ja ihmettelee, mitä on tehty, jotta aineisto soveltuisi tällaiseen malliin, jos se ei ole perusteltua.︎

3. Imai asettaa artikkelinsa saataville verkkosivuillaan.︎

4. MASS on korvattu muilla jo yli vuosikymmenen ajan tässä vaiheessa, ja sillä on enimmäkseen vain taipumus sotkea tidyverseä ja muita paketteja ladattaessa. Se on hieno paketti (ja oli aikoinaan loistava), mutta jos sitä haluaa käyttää paketissa, olisi hyvä olla lataamatta sitä (tai muita paketteja) ympäristöön vain yhden tai kahden funktion käyttöä varten. Enimmäkseen näen sitä käytettävän vain mvrnormiin (monimuuttujainen normaalijakauma) ja glm.nb:hen, mutta on muitakin paketteja, joissa on kyseisiä toimintoja, jotka tarjoaisivat lisähyötyjä, eikä naamioida dplyr-funktioita, jotka ovat R-yhteisössä yleisimmin käytettyjen joukossa.︎

5. brms työskentelee sen parissa.︎

6. Koodissa olevat vaihtoehdot ovat olemassa tukahduttamaan/minimoimaan sitä, mitä voi olla.︎