Free Press

Obsah

Aktualizováno 04. srpna 2020. Kód lze stáhnout zde.

Úvod

V některých situacích můžeme uvažovat o nepřímém vlivu nějaké proměnné na výsledek nebo výsledek. Jako příklad lze uvést špatné životní podmínky doma v dětství, které mohou snížit výsledky učení ve škole, což má následně negativní vliv na pozdější kvalitu života, například na celoživotní příjmy. V jiném případě můžeme uvažovat o jediné proměnné shromážděné ve více časových bodech, takže existuje vliv proměnné v čase 1 na čas 2 a v čase 2 na čas 3. To znamená, že proměnná v čase 1 má vliv na čas 3. Základní myšlenka je následující:

\

Jinými slovy, \(\mathcal{A}\) vede k \(\mathcal{B}\) a pak \(\mathcal{B}\) vede k \(\mathcal{C}\). V modelech zprostředkování předpokládáme intervenující proměnnou mezi normální kovariátou \(\rightarrow\) a výsledkem, který bychom mohli mít ve standardním regresním nastavení, a tyto modely nám umožňují zkoumat takové chování. Ve výše uvedeném případě je intervenující proměnnou nebo mediátorem \(\mathcal{B}\). Často se stává, že stále můžeme mít přímý vliv \(\mathcal{A}\) na \(\mathcal{C}\), ale stejně jako u modelu obecně by to bylo teoreticky zdůvodněno.

Analýza mediace je velmi populární ve společenskovědních disciplínách, i když se na ně v žádném případě neomezuje, a obvykle se provádí pod rouškou modelování strukturálních rovnic (SEM), což je samo o sobě specifické zaměření grafických modelů obecněji1. Grafický model mediačního modelu může vypadat následovně.

V tomto případě a a b odrážejí nepřímou cestu vlivu \(\mathrm{X}\) na výsledek prostřednictvím mediátoru, zatímco c' je přímý účinek \(\mathrm{X}\) na výsledek po odstranění nepřímé cesty (c by byl účinek před zavedením nepřímého účinku a c – c' se rovná nepřímému účinku). Celkový účinek \(\mathrm{X}\) je kombinací nepřímých a přímých účinků.

Měl bych poznamenat několik věcí na základě toho, co vidím při konzultacích v desítkách oborů. Pro začátek se zdá, že jen velmi málo lidí, kteří si myslí, že potřebují model zprostředkování, ho skutečně potřebuje. Pokud například nemůžete uvažovat o svém modelu v časových nebo fyzikálních termínech tak, že \(\mathrm{X}\) nutně vede k mediátoru, který pak nutně vede k výsledku, pravděpodobně mediační model nepotřebujete. Pokud byste viděli šipky vedoucí oběma směry, opět takový model pravděpodobně nepotřebujete. Také pokud si při popisu vašeho modelu všichni myslí, že mluvíte o interakci (neboli moderování), asi ho nepotřebujete. A konečně, jak lze předpokládat, pokud neexistuje silná korelace mezi klíčovými proměnnými (\(\mathrm{X}\)) a mediátorem (cesta a) a pokud neexistuje silná korelace mezi mediátorem a výsledkem (výsledky) (cesta b), pravděpodobně to nepotřebujete. I když vám nic nebrání provést mediační analýzu, bez těchto předpokladů budete mít téměř jistě slabý a pravděpodobně více matoucí model, než byste jinak měli.

Krátce řečeno, mediace funguje nejlépe, když mezi proměnnými existují silně implikované kauzální vazby. I v takovém případě je třeba takový model porovnat s jednodušším modelem bez mediace2. V každém případě existuje několik velmi snadných způsobů, jak takové modely v R zkoumat, a to je cílem tohoto článku, jen chceme ukázat, jak můžete začít.

Data

Pro demonstraci mediačních modelů s různými balíčky použijeme datovou sadu jobs, která je součástí balíčku mediation. Zde je jeho popis.

Studie intervence při hledání zaměstnání (JOBS II). JOBS II je randomizovaný terénní experiment, který zkoumá účinnost intervence v oblasti hledání zaměstnání na nezaměstnané pracovníky. Program je navržen tak, aby nejen zvýšil opětovné zaměstnávání nezaměstnaných, ale také zlepšil duševní zdraví uchazečů o zaměstnání. V rámci terénního experimentu JOBS II obdrželo 1 801 nezaměstnaných pracovníků předběžný dotazník a poté byli náhodně rozděleni do terapeutické a kontrolní skupiny. Osoby v léčebné skupině se účastnily seminářů zaměřených na pracovní dovednosti. Na seminářích se respondenti učili dovednostem při hledání zaměstnání a strategiím zvládání neúspěchů při hledání zaměstnání. Účastníci kontrolní skupiny obdrželi brožuru popisující tipy pro hledání zaměstnání. V následných rozhovorech byly dvě klíčové výsledné proměnné kontinuální míra depresivních symptomů založená na Hopkinsově kontrolním seznamu symptomů a binární proměnná, která vyjadřovala, zda se respondent stal zaměstnancem.

Popis proměnných v této ukázce. K dispozici jsou i další, se kterými si možná budete chtít pohrát.

• econ_hard: Úroveň ekonomických potíží před ošetřením s hodnotami od 1 do 5.
• sex: Indikátorová proměnná pro pohlaví. 1 = žena
• věk: Věk v letech.
• educ: Faktor s pěti kategoriemi pro dosažené vzdělání.
• job_seek: Stejnosměrná škála měřící úroveň sebepojetí při hledání zaměstnání s hodnotami od 1 do 5. Zprostředkující proměnná:
• depress2: Míra depresivních symptomů po léčbě. Výsledná proměnná.
• treat: Indikátorová proměnná pro to, zda byl účastník náhodně vybrán do školicího programu JOBS II. 1 = přiřazení k účasti.

data(jobs, package = 'mediation')

Model

Vzhledem k těmto údajům jsou modely pro mediátor a výsledek následující:

\

Očekáváme tedy, že trénink pracovních dovedností bude mít negativní vliv na depresi (tj. zvýšení pohody), ale alespoň částečně to bude způsobeno pozitivním vlivem na hledání práce.

Jako grafický model bychom to mohli stručně znázornit následovně.

Balíčky

Podíváme se na následující balíčky, abychom si ukázali, jak lze v R provádět mediační analýzu:

• mediation
• lavaan
• psych
• brms

Přestože se zaměříme na tyto balíky, upozorním i na některé další alternativy, včetně Pythonu a Staty.

mediation

Začneme balíčkem mediation, protože k jeho vedení v podstatě není zapotřebí více programátorských schopností, než které člověk již má z provádění standardních regresních modelů v R. Balík poskytuje průměrný kauzální mediační efekt, který je ze souboru nápovědy a Imaiových článků3 definován takto:

Průměrný kauzální mediační efekt (ACME) představuje očekávaný rozdíl v potenciálním výsledku, když mediátor nabyl hodnoty, která by se realizovala za léčebného stavu oproti kontrolnímu stavu, zatímco samotný stav léčby je konstantní.

Všimněte si, že tato definice je zaměřena na očekávané nebo předpokládané hodnoty podmíněné hodnotou léčby. Toto pojetí kontrafaktuálů, neboli jak by pozorování vypadalo za opačného nastavení, má v tomto bodě v modelování dlouhou historii. Představte si to tak, že pokud je člověk v léčebné skupině, měl by určitou hodnotu mediátoru a vzhledem k ní by pak měl určitou očekávanou hodnotu výsledku. Stejné pozorování bychom však mohli představit i jako pozorování v kontrolní skupině a stejně tak posoudit vliv na výsledek prostřednictvím mediátoru. Potenciální výsledky můžeme posuzovat při zachování konstantní léčby. Uvažování o změnách výsledku vzhledem k hodnotě mediátoru nevytváří žádný předpoklad o typu modelu. Takto je mediační balíček schopen zahrnout různé modely pro mediátor vs. výsledek. Například mediátor může být binární, což vyžaduje logistický regresní model, zatímco model výsledku může být model přežití.

V našem příkladu se budeme držet standardních (normálních) lineárních modelů. Všimněte si také, že i když je naše léčba binární proměnnou, zobecňuje se to na spojitý případ, kdy uvažujeme výsledek pohybu o jednu jednotku na „léčbě“. Aby balík pro mediaci fungoval, jednoduše spustíme naše příslušné modely pro mediátor a výsledek a poté použijeme funkci mediate, abychom získali konečný výsledek.

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)

Výše uvedené výsledky ukazují, že ACME se statisticky neliší od nuly, resp. není zprostředkován. Průměrný přímý účinek je záporný, ale stejně tak není statisticky významný, stejně jako celkový účinek (nepřímý + přímý účinek). Rovněž je uveden soi disant „proportion mediated“, což je poměr nepřímého účinku k celkovému. Nejedná se však o podíl, může být dokonce záporný, a tak se většinou jedná o bezvýznamné číslo.

Pros

• Standardní modely a syntaxe R
• Více typů modelů pro mediátor i výsledek
• Poskytuje více výsledků současně
• Dobrá dokumentace a související články jsou volně dostupné
• Může provádět ‚moderovanou‘ mediaci

Omezení

.
• Použití MASS4
• Jednoduché modely náhodných efektů
• Funkčnost je možná omezena některými složitostmi modelu
• Není možné použít latentní proměnnou

lavaan

Ve specifickém případě, kdy mediační i výsledkové modely jsou standardní lineární modely s normálním rozdělením cílové proměnné, nepřímý účinek odpovídá součinu cest a a b v předchozím diagramu. Přímý účinek je cesta c'. Porovnání samostatného přímého účinku, který bychom mohli nazvat c, s tímto odhadnutým přímým účinkem v mediačním modelu c', je takové, že c - c' = a*b. To, co bylo uvedeno dříve, by nyní mohlo být jasnější, pokud je buď a, nebo b téměř nulový, pak nepřímý účinek může být pouze téměř nulový, takže je rozumné takové vztahy předem prozkoumat.

Tento přístup založený na součinu cest (nebo rozdílu koeficientů) použijeme s balíčkem lavaan a v podstatě je to v době psaní tohoto článku náš jediný způsob, jak na to jít. lavaan je speciálně zaměřen na modelování strukturálních rovnic, jako je faktorová analýza, růstové modely a mediační modely, jaké zde provádíme, a je pro takové modely velmi doporučován. Je sice omezen na standardní případ lineárního modelu pro posouzení mediace, ale jako jediný z našich nástrojů dokáže snadno začlenit latentní proměnné5. Například bychom mohli mít náš výsledek deprese jako latentní proměnnou, která je základem jednotlivých položek dotazníku. Kromě toho bychom mohli zahrnout i více mediátorů a více výsledků.

Abychom zachovali to, co jsme probírali, budu cesty a, b a c' v laváži označovat podle toho, jak byly zobrazeny dříve. Jinak se lavaan používá velmi snadno a v případě pozorovaných proměnných používá standardní zápis vzorců R pro modely. Kromě toho definujeme zájmové efekty, které chceme vypočítat, pomocí operátoru :=. Celý model zadáme jako jednoduchý řetězec znaků a poté pomocí funkce sem provedeme analýzu.

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

Vidíme stejný výstup jako předtím a můžeme porovnat náš parametr indirect s ACME, který jsme měli předtím, efekt direct se porovnává s ADE a total se porovnává s předchozím celkovým efektem. Hodnoty jsou v podstatě stejné.

Všimněte si také, že výstup ukazuje hodnotu \(R^2\) pro oba modely. V případě job_seek vidíme, že důvodem, proč nenacházíme mnoho mediátorů, je to, že zúčastněné kovariáty nevysvětlují na začátku žádnou variabilitu mediátoru. Předběžné šetření by nám v tomto případě ušetřilo práci.

Pros

• Může zpracovávat více mediátorů
• Může zpracovávat více „ošetření“
• Může zpracovávat více výsledků
• Může používat latentní proměnné
• Nějaká víceúrovňová podpora
• Může provádět moderovanou mediaci a mediovanou moderace (i když ne pro latentní proměnné)

Omezení

• Vyžaduje další kódování pro odhad nepřímého účinku
• Jednotlivé náhodné efekty
• Mohly by modely zahrnovat binární nebo ordinální proměnné pro mediátor/výsledky, neexistuje žádný přímočarý způsob, jak v těchto nastaveních vypočítat nepřímý účinek způsobem balíčku mediace.

piecewiseSEM

Balíček piecewiseSEM pracuje velmi podobně jako balíček mediation. Příjemné na něm oproti balíčku mediation je to, že piecewiseSEM umí pracovat s dalšími typy modelů a také poskytuje další výstupy (např. standardizované výsledky), další možnosti (např. více skupin, korelovaná rezidua) a vizualizaci modelu.

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

Můžeme využít jeho možnosti vykreslování a vytvořit rychlou vizualizaci modelu.

plot(mediation_result)

Naneštěstí v současné době neexistuje žádný automatický způsob výpočtu nepřímých účinků, takže by bylo nutné provést bootstrapování výsledků ručně.

Pros

• Standardní modely a syntaxe R
• Více typů modelů pro mediátor i výsledek
• Některé výsledky ve stylu SEM (např. fit, standardizované koeficienty, AIC)
• Rychlý graf výsledků
• Může zpracovávat více mediátorů, „ošetření“, a výsledky

Omezení

• Nepočítá automaticky nepřímé účinky
• Neumožňuje výpočet latentních proměnných

psych

Balíček psych využívá toho, že v případě standardního lineárního modelu lze získat výsledky prostřednictvím příslušných regresních modelů pouze na základě kovariančních matic. Je velmi podobný balíku lavaan, i když používá přístup obyčejných nejmenších čtverců na rozdíl od přístupu maximální věrohodnosti. Příjemná je zde syntaxe, která umožňuje zaměřit se pouze na efekt, který vás zajímá, nebo zahrnout vše, což je příjemné, pokud by vás zajímaly i nepřímé efekty pro ekonomické potíže, věk a pohlaví.

Pro tuto ukázku použijeme vyčištěnou verzi s použitím - namísto + pro efekty bez léčby. To pouze znamená, že jsou zahrnuty do modelů, ale výsledky týkající se jich se nezobrazují. Mediátor je identifikován pomocí (). Dalším bonusem je rychlý graf výsledků, který ukazuje rozdíl mezi neupravenými a upravenými přímými účinky a příslušný bootstrapovaný interval.

library(psych)mediation_psych = mediate( depress2 ~ treat + (job_seek) - econ_hard - sex - age, data = jobs, n.iter = 500)

mediation_psych

Mediation/Moderation Analysis Call: mediate(y = depress2 ~ treat + (job_seek) - econ_hard - sex - age, data = jobs, n.iter = 500)The DV (Y) was depress2* . The IV (X) was treat* . The mediating variable(s) = job_seek* . Variable(s) partialled out were econ_hard sex ageTotal effect(c) of treat* on depress2* = -0.06 S.E. = 0.05 t = -1.24 df= 895 with p = 0.21Direct effect (c') of treat* on depress2* removing job_seek* = -0.04 S.E. = 0.15 t = 14.91 df= 893 with p = 4.6e-45Indirect effect (ab) of treat* on depress2* through job_seek* = -0.02 Mean bootstrapped indirect effect = -0.02 with standard error = 0.01 Lower CI = -0.04 Upper CI = 0.01R = 1.07 R2 = 1.15 F = -3510.4 on 2 and 893 DF p-value: 1 To see the longer output, specify short = FALSE in the print statement or ask for the summary

summary(mediation_psych)

Call: mediate(y = depress2 ~ treat + (job_seek) - econ_hard - sex - age, data = jobs, n.iter = 500)Direct effect estimates (traditional regression) (c') depress2* se t df ProbIntercept 2.21 0.15 14.91 893 4.60e-45treat -0.04 0.04 -0.93 893 3.55e-01job_seek -0.24 0.03 -8.50 893 8.14e-17R = 1.07 R2 = 1.15 F = -3510.4 on 2 and 893 DF p-value: 1 Total effect estimates (c) depress2* se t df Probtreat -0.06 0.05 -1.24 895 0.215 'a' effect estimates job_seek se t df ProbIntercept 3.67 0.13 29.33 894 5.65e-133treat 0.07 0.05 1.27 894 2.03e-01 'b' effect estimates depress2* se t df Probjob_seek -0.24 0.03 -8.5 894 7.83e-17 'ab' effect estimates (through mediators) depress2* boot sd lower uppertreat -0.02 -0.02 0.01 -0.04 0.01

Stejné výsledky, jiný obal, ale možná zatím nejjednodušší cesta, protože vyžadovala pouze jedno volání funkce. Balíček Psych si jako bonus poradí i s více prostředníky a výsledky.

Pros

• Nejjednodušší syntaxe, v podstatě jednořádkový model
• Rychlý graf výsledků
• Může zpracovávat více mediátorů, „ošetření“, a výsledky
• Může provádět ‚moderovanou‘ mediaci

Omezení

• Omezeno na standardní lineární model (lm)
• Použití MASS

brms

Pro naši další ukázku se dostáváme k podle mého názoru nejvýkonnějšímu balíku, brms. Název znamená Bayesian Regression Modeling with Stan a Stan je výkonný pravděpodobnostní programovací jazyk pro Bayesovskou analýzu. Nebudu zabíhat do podrobností o bayesovské analýze, ale klidně se podívejte na můj dokument, který to dělá.

Obvykle postupujeme stejně jako předtím, tedy zadáváme model mediátoru a model výsledku. brms neumí pro mediační analýzu nic zvláštního, ale jeho funkce hypotézy nám může umožnit testovat přístup založený na součinu cest. Navíc balík sjstats nám v podstatě poskytne výsledky stejným způsobem, jakým to dělá balík mediace, a když na to přijde, balík mediace je stejně v podstatě pokusem o bayesovské řešení pomocí frekvenčních metod. Pokud bychom skutečně měli různá rozdělení pro výsledek a mediátor, poměrně snadno bychom získali tyto průměrné hodnoty predikce a jejich rozdíly, protože bayesovské přístupy vždy uvažují o posteriorních predikčních rozděleních. V každém případě je zde 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()

Ve výstupu je cokoli s jobseek_* výsledkem pro model mediátoru, zatímco depress2_* pro výsledek. V tomto bodě máme stále stejný příběh, ale díky bayesovskému přístupu se můžeme podívat na zábavnější věci. Například vidíme, že ve skutečnosti dobře nezachycujeme zkreslení výsledku deprese. Naše předpovídané hodnoty oproti pozorovaným se tak docela neshodují. U mediátoru jsme na tom o něco lépe, ale u některých našich předpovědí založených na modelu jsme možná stále trochu vysoko.

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

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

Pros

• Přímočará syntaxe
• Extrémně silná-. Modely jsou většinou omezeny na vlastní představivost
• V podstatě dělá to, co aproximuje mediační balík
• Všechny výhody bayesovské inference: diagnostiku, posteriorní prediktivní kontroly, porovnávání modelů atd.

Omezení

• Nejjednodušší odhad
• „Ruční“ výpočty potřebné pro překročení standardního lineárního modelu, ale to je z bayesovského hlediska již běžný přístup
• Vyžaduje se určitý komfort s bayesovským přístupem

Větší složitost

Některé z uvedených balíků si poradí se složitějšími modely nebo poskytují další přístupy ke zkoumání nepřímých účinků.

Interakce

Některé modely zahrnují interakce buď pro mediační model, nebo pro výsledek, a to se bohužel často označuje jako zprostředkovaná moderace nebo moderovaná mediace. Osobně nevidím výhodu v dávání nejednoznačných názvů něčemu, co by jinak mohlo být přímočarým konceptem (i když stále ne tak přímočarým modelem), ale tato loď už dávno odplula. Nebudu zabíhat do podrobností, ale jde o to, že někde v modelu můžete mít interakční člen a interakce se může týkat proměnné léčby, mediátoru nebo obou.

Stačí říct, že vzhledem k tomu, že používáme standardní modelovací nástroje jako lm a jejich rozšíření, je začlenění interakcí triviální pro všechny výše uvedené balíčky, ale přístup typu součin cest neplatí (a*b != c').

Generalizované lineární modely

V některých případech může být náš mediátor nebo výsledek binární, počítaný nebo něco, kde předpokládat normální rozdělení nemusí být nejlepší nápad. Nebo můžeme chtít zkoumat nelineární vztahy mezi léčbou/mediátorem/výsledkem. Nebo můžeme mít data, která mají korelovaná pozorování, jako jsou opakovaná měření apod. Zejména balík mediation se tím pyšní, ale brms umí všechno, co umí, a ještě víc, i když si možná budete muset dát trochu víc práce, abyste výsledek skutečně vypočítali. lavaan skutečně umí omezenou sadu modelů pro binární a ordinální proměnné, ale získání vhodného nepřímého odhadu by vyžadovalo velmi zdlouhavý ruční přístup.

Chybějící data

Při práci s takovými daty, zejména ve společenských vědách, často chybí údaje o některé z kovariát. Někdy je můžeme vypustit, pokud jich není příliš mnoho, ale v jiných případech s tím budeme chtít něco udělat. Balíčky lavaan, psych a brms poskytují jeden nebo více způsobů, jak se s touto situací vypořádat (např. vícenásobná imputace).

Alternativy

Modely jsme zobrazovali jako sítě uzlů, které spojují oblouky/hrany/cesty. Naše diskuse se točí kolem tzv. směrovaných acyklických grafů (DAG), kde šipky mohou směřovat pouze jedním směrem bez zpětných vazeb. Výsledek jakékoli výsledné proměnné je funkcí šipek, které jí předcházejí, a je podmíněně nezávislý na ostatních. V některých teoretických modelech může být tento požadavek zmírněn a v jiných nemusí být šipky vůbec žádné, tj. jsou neusměrněné, takže nás zajímají pouze spojení (např. u některých sociálních sítí).

bnlearn

Balík bnlearn umožňuje zkoumat směrované, částečně směrované a neusměrněné grafy. Pokud jde o DAGy, můžeme jej použít v podstatě ke zopakování modelů zprostředkování, o kterých jsme hovořili. Příjemné však je, že tento balíček bude efektivně testovat cesty pro zařazení, nikoli je předpokládat, ale stále můžeme podle potřeby ukládat teoretická omezení. Nejenže pak můžeme hledat cesty, které nás zajímají, principiálním způsobem pomocí bayesovských sítí a Pearlovy teorie kauzálních grafů jako základu, ale také budeme mít k dispozici nástroje, jak se dále vyhnout nadměrnému fitování pomocí křížové validace.

Pro počáteční model se ujistíme, že existují cesty mezi léčbou – mediátorem, léčbou – výsledkem a mediátorem – výsledkem (whitelist). Zakážeme nesmyslné cesty, jako například mít šipky k léčbě (která byla náhodně přiřazena), pohlaví, ekonomické potíže a věk (černá listina). Jinak se podíváme, co naznačují data.

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)

V grafu vidíme, že se situace trochu změnila. Například věk nyní souvisí pouze se soběstačností při hledání zaměstnání a pohlaví má vliv pouze na depresi.

Kdybychom omezili cesty pouze na ty, které jsou v našich předchozích příkladech, dostali bychom stejné výsledky.

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 

Hlavní věc, které si musíme všimnout, je, že odhadované parametry se rovnají tomu, co jsme dostali s předchozími balíčky. Je to v podstatě ekvivalentní použití lavaan s výchozím odhadem maximální věrohodnosti.

Použijeme-li jako faktory léčbu a pohlaví, bnlearn vytvoří podmíněné modely, které se liší v závislosti na přijaté hodnotě faktoru. Jinými slovy, jeden model by byl zvlášť pro případ, kdy treatment == 'treatment', a druhý pro případ, kdy treatment == control. V našem případě by to bylo totožné s tím, kdybychom vše nechali interagovat s léčbou, např. lm( job_seek ~ treat * (econ_hard + sex + age)), a stejně tak pro model deprese. To by se rozšířilo na potenciálně jakoukoli binární proměnnou (např. včetně pohlaví). Pokud je mediátor binární proměnnou, je to pravděpodobně to, co bychom chtěli udělat.

Python

Ředitel CSCAR Kerby Shedden uspořádal workshop o mediačních modelech v jazyce Python, takže zde ukazuji implementaci statsmodels. Vychází z přístupu Imai, a tak ji lze považovat za pythonovskou verzi mediačního balíčku. Výstup je v podstatě stejný jako při použití léčby jako faktorové proměnné, kdy získáte samostatné výsledky pro každou kategorii léčby. To je pro naši ukázku zbytečné, takže stačí porovnat „průměrné“ výsledky s výsledky předchozího mediačního balíčku.

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

Nakonec uvádím možnost ve Stata pomocí jejího příkazu sem. Stata v tomto příkladu umožňuje snadno získat nepřímé efekty, ale dělá to pro každý kovariát, takže výstup je přinejmenším poněkud mnohomluvný6. Ti, kdo pracují se systémem Stata, nepotřebují k získání těchto druhů výsledků samostatný balík SEM.

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)

Shrnutí

Modely s nepřímými účinky vyžadují pečlivé teoretické zvážení, aby je bylo možné použít pro analýzu dat. Pokud je však model vhodný pro vaši datovou situaci, je poměrně snadné získat výsledky pomocí různých balíčků v R. Navíc k provedení analýzy s nepřímými účinky není nutné použít balíček pro modelování strukturálních rovnic a ve skutečnosti se lze dostat daleko pomocí standardní syntaxe R. V případě, že je model vhodný pro vaši datovou situaci, je možné získat výsledky pomocí různých balíčků. Pro striktně pozorované, tj. žádné latentní, proměnné není žádný nástroj SEM nutný, nebo dokonce doporučený.

\

Souhrnné srovnání balíčků

Následující tabulka může pomoci při rozhodování, který balíček použít pro své potřeby vzhledem k teoretickým úvahám.

1. Mám mnohem podrobnější dokument o SEM, včetně mediační analýzy.︎

2. Z nějakého důvodu se s tím v praxi příliš nesetkáte a člověk si klade otázku, co bylo uděláno pro to, aby data takový model umožňovala, když to nebylo opodstatněné.︎

3. Imai své články zpřístupňuje na svých webových stránkách.︎

4. MASS je v tuto chvíli již více než deset let nahrazen jinými a většinou má tendenci vám při načítání zaneřádit tidyverse a další balíky. Je to fajn balík (a ve své době byl skvělý), ale pokud ho chcete používat v balíku, bylo by dobré ho (nebo jiné balíky) nenačítat do prostředí jen kvůli použití jedné nebo dvou funkcí. Většinou ho vidím používat jen pro mvrnorm (vícerozměrné normální rozdělení) a glm.nb, ale existují i jiné balíky s touto funkcí, které by poskytly další výhody, a nemaskují funkce dplyr, které patří v komunitě R k nejpoužívanějším.︎

5. brms na tom pracuje.︎

6. Volby v kódu jsou tam proto, aby potlačily/minimalizovaly to, co může být.︎