Free Press

Table of Contents

Updated August 04, 2020. Kod można pobrać tutaj.

Wprowadzenie

W niektórych sytuacjach możemy rozważyć pośredni wpływ jakiejś zmiennej na wynik lub rezultat. Na przykład, złe warunki życia w domu w dzieciństwie mogą zmniejszyć wyniki w nauce w szkole, które następnie mają negatywny wpływ na późniejszą jakość życia, na przykład, zarobki przez całe życie. W innym przypadku możemy rozważyć pojedynczą zmienną zebraną w wielu punktach czasowych, tak że istnieje efekt zmiennej w czasie 1 na czas 2, a czas 2 na czas 3. Podstawową ideą jest coś takiego:

Innymi słowy, \(\mathcal{A}} prowadzi do \(\mathcal{B}}), a następnie \(\mathcal{B}} prowadzi do \(\mathcal{C}). W modelach mediacyjnych wprowadzamy zmienną pośredniczącą pomiędzy normalną ścieżką wyniku zmiennej, którą moglibyśmy mieć w standardowym ustawieniu regresji, a te modele pozwalają nam badać takie zachowania. W powyższym przypadku zmienna pośrednicząca, lub mediator, jest \(\). Często zdarza się, że nadal możemy mieć bezpośredni wpływ \(\mathcal{A}} na \(\mathcal{C}}, ale tak jak w przypadku modelu w ogóle, byłoby to umotywowane teoretycznie.

Analiza mediacji jest bardzo popularna w dyscyplinach nauk społecznych, choć w żaden sposób nie ogranicza się do nich, i zazwyczaj przeprowadzana jest pod przykrywką modelowania równań strukturalnych (SEM), które samo w sobie jest specyficzną orientacją modeli graficznych bardziej ogólnie1. Graficzny model mediacji może wyglądać następująco.

W tym przypadku, a i b odzwierciedlają pośrednią ścieżkę wpływu efektu \(\) na wynik poprzez mediator, podczas gdy c' jest bezpośrednim efektem \(\) na wynik po usunięciu ścieżki pośredniej (c byłby efektem przed założeniem efektu pośredniego, a c – c' równa się efektowi pośredniemu). Całkowity efekt jest połączonym efektem pośrednim i bezpośrednim.

Powinienem zauważyć kilka rzeczy opartych na tym, co widzę w konsultacjach w dziesiątkach dyscyplin. Na początek, wydaje się, że bardzo niewielu ludzi, którzy myślą, że potrzebują modelu mediacyjnego, rzeczywiście go potrzebuje. Jeśli mógłbyś zobaczyć strzałki idące w dowolnym kierunku, ponownie, prawdopodobnie nie potrzebujesz takiego modelu. Również, jeśli opisując twój model, wszyscy myślą, że mówisz o interakcji (a.k.a. moderacji), możesz tego nie potrzebować. I wreszcie, jak można podejrzewać, jeśli nie ma silnej korelacji między kluczowymi zmiennymi (\) a mediatorem (ścieżka a), i jeśli nie ma silnej korelacji między mediatorem a wynikiem(ami) (ścieżka b), prawdopodobnie nie potrzebujesz tego. Podczas gdy nic nie powstrzyma Cię przed przeprowadzeniem analizy mediacji, bez takich warunków wstępnych, prawie na pewno będziesz miał słaby i prawdopodobnie bardziej mylący model, niż miałbyś w przeciwnym razie.

W skrócie, mediacja działa najlepiej, gdy istnieją silnie implikowane związki przyczynowe między zmiennymi. Nawet wtedy, taki model powinien być porównywany z prostszym modelem bez mediacji2. W każdym razie, istnieje kilka bardzo łatwych sposobów na zbadanie takich modeli w R, i to jest cel tutaj, aby zademonstrować, jak można zacząć.

Dane

Do demonstracji modeli mediacji z różnymi pakietami, użyjemy zestawu danych jobs, który jest dostarczany z pakietem mediation. Oto opis.

Badanie Interwencyjne Poszukiwania Pracy (JOBS II). JOBS II jest randomizowanym eksperymentem terenowym, który bada skuteczność interwencji szkolenia zawodowego dla bezrobotnych pracowników. Program został zaprojektowany tak, aby nie tylko zwiększyć ponowne zatrudnienie wśród bezrobotnych, ale także poprawić zdrowie psychiczne osób poszukujących pracy. W eksperymencie terenowym JOBS II 1 801 bezrobotnych pracowników otrzymało kwestionariusz wstępny, a następnie zostało losowo przydzielonych do grupy leczonej i kontrolnej. Osoby z grupy leczonej uczestniczyły w warsztatach umiejętności zawodowych. W trakcie warsztatów respondenci uczyli się umiejętności poszukiwania pracy oraz strategii radzenia sobie z niepowodzeniami w procesie poszukiwania pracy. Osoby z grupy kontrolnej otrzymały broszurę opisującą wskazówki dotyczące poszukiwania pracy. W wywiadach kontrolnych dwiema kluczowymi zmiennymi wynikowymi były ciągła miara objawów depresji oparta na kwestionariuszu Hopkins Symptom Checklist oraz zmienna binarna określająca, czy respondent został zatrudniony.

Tutaj znajduje się opis zmiennych w tej demonstracji. Są dostępne inne zmienne, którymi również możesz chcieć się pobawić.

• econ_hard: Poziom trudności ekonomicznych przed leczeniem, z wartościami od 1 do 5.
• sex: Zmienna wskaźnikowa dla płci. 1 = kobieta
• wiek: Wiek w latach.
• educ: Czynnik z pięcioma kategoriami dla poziomu wykształcenia.
• job_seek: Ciągła skala mierząca poziom poczucia własnej skuteczności w poszukiwaniu pracy o wartościach od 1 do 5. Zmienna pośrednicząca.
• depress2: Pomiar objawów depresji po zakończeniu leczenia. Zmienna wyniku.
• treat: Zmienna wskaźnikowa określająca, czy uczestnik został losowo wybrany do programu szkoleniowego JOBS II. 1 = przypisanie do uczestnictwa.

data(jobs, package = 'mediation')

Model

Przy tych danych modele dla mediatora i wyniku są następujące:

Tak więc oczekujemy, że szkolenie w zakresie umiejętności zawodowych będzie miało negatywny wpływ na depresję (tj. wzrost dobrego samopoczucia), ale przynajmniej część z tego będzie wynikać z pozytywnego wpływu na poszukiwanie pracy.

Jako model graficzny możemy to zwięźle przedstawić w następujący sposób.

Pakiety

Przyjrzymy się następującym pakietom, aby zademonstrować, jak można przeprowadzić analizę mediacji w R:

• mediation
• lavaan
• psych
• brms

Pomimo że na nich się skupimy, zwrócę też uwagę na kilka innych alternatyw, w tym Pythona i Statę.

mediacja

Zaczniemy od pakietu mediation, ponieważ w zasadzie nie wymaga on więcej umiejętności programistycznych do przeprowadzenia niż te, które posiada się już po uruchomieniu standardowych modeli regresji w R. Pakiet zapewnia średni przyczynowy efekt mediacji, zdefiniowany w następujący sposób w pliku pomocy i artykułach Imai’a3:

Średni przyczynowy efekt mediacji (ACME) reprezentuje oczekiwaną różnicę w potencjalnym wyniku, gdy mediator przyjął wartość, która zrealizowałaby się w warunkach leczenia w przeciwieństwie do warunków kontrolnych, podczas gdy sam status leczenia jest utrzymywany na stałym poziomie.

Zauważ, jak ta definicja skupia się na oczekiwanych lub przewidywanych wartościach warunkowych w stosunku do wartości leczenia. To pojęcie kontrfaktycznych, lub jak wyglądałaby obserwacja w przeciwnym ustawieniu, ma długą historię w modelowaniu w tym momencie. Pomyśl o tym w ten sposób, jeśli ktoś jest w grupie leczonej, miałby określoną wartość dla mediatora, i biorąc to pod uwagę, miałby wtedy określoną oczekiwaną wartość dla wyniku. Jednakże, moglibyśmy założyć tę samą obserwację jako bycie w grupie kontrolnej, jak również, i ocenić wpływ na wynik przez mediatora tak samo. Możemy ocenić potencjalne wyniki przy zachowaniu stałej wartości leczenia. Myślenie o zmianach wyników biorąc pod uwagę wartość mediatora nie czyni żadnych założeń co do typu modelu. W ten sposób pakiet mediacyjny jest w stanie włączyć różne modele dla mediatora i wyniku. Na przykład, mediator może być binarny, co wymaga modelu regresji logistycznej, podczas gdy model wyniku może być modelem przeżycia.

W naszym przykładzie, będziemy trzymać się standardowych (normalnych) modeli liniowych. Zauważ również, że podczas gdy nasze leczenie jest zmienną binarną, generalizuje to do przypadku ciągłego, gdzie rozważamy wynik ruchu o jedną jednostkę na „leczeniu”. Aby pakiet mediacji działał, po prostu uruchamiamy nasze odpowiednie modele dla mediatora i wyniku, a następnie używamy funkcji mediate, aby uzyskać ostateczny wynik.

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)

Powyższe wyniki pokazują, że ACME nie jest statystycznie różny od zera, czyli braku mediacji. Średni efekt bezpośredni jest ujemny, ale również nie jest statystycznie istotny, podobnie jak efekt całkowity (efekt pośredni + bezpośredni). Podano również soi disant „proporcję pośrednictwa”, która jest stosunkiem efektu pośredniego do całkowitego. Jednak nie jest to proporcja, a nawet może być ujemna, a więc jest to w większości przypadków liczba bez znaczenia.

Pros

• Standardowe modele i składnia R
• Wielokrotne typy modeli zarówno dla mediatora, jak i wyniku
• Dostarcza wiele wyników jednocześnie
• Dobra dokumentacja i powiązane artykuły są swobodnie dostępne
• Może wykonywać „moderowane” mediacje

Ograniczenia

.
• Użycie MASS4
• Proste modele efektów losowych
• Funkcjonalność może być ograniczona przy pewnych złożonościach modelu
• Brak możliwości zmiennych ukrytych

lavaan

W szczególnym przypadku, gdy zarówno modele mediacji, jak i wyniku są standardowymi modelami liniowymi z rozkładem normalnym dla zmiennej docelowej, efekt pośredni jest równoważny iloczynowi ścieżek a i b na poprzednim wykresie. Efektem bezpośrednim jest ścieżka c'. Porównanie samodzielnego efektu bezpośredniego, który moglibyśmy nazwać c, z tym oszacowanym efektem bezpośrednim w modelu mediacyjnym c', jest takie, że c - c' = a*b. Jeśli a lub b są bliskie zeru, to efekt pośredni może być tylko bliski zeru, więc rozsądnie jest zbadać takie zależności wcześniej.

To podejście oparte na iloczynie ścieżek (lub różnicy współczynników) jest tym, które zastosujemy z pakietem lavaan, i w rzeczywistości, od momentu napisania tego tekstu, jest to nasz jedyny sposób postępowania. lavaan jest szczególnie ukierunkowany na modelowanie równań strukturalnych, takie jak analiza czynnikowa, modele wzrostu i modele mediacji, takie jak te, które tutaj przeprowadzamy, i jest wysoce zalecany dla takich modeli. Chociaż jest on ograniczony do standardowego przypadku modelu liniowego do oceny mediacji, jest to jedyne z naszych narzędzi, które może łatwo włączyć zmienne ukryte5. Na przykład, moglibyśmy mieć nasz wynik depresji jako zmienną latentną leżącą u podstaw poszczególnych pozycji kwestionariusza. Ponadto, moglibyśmy również włączyć wiele mediatorów i wiele wyników.

Aby zachować rzeczy tak, jak rozmawialiśmy, będę oznaczał ścieżki a, b i c' w lavaan zgodnie z tym, jak zostały one przedstawione wcześniej. Poza tym lavaan jest bardzo łatwy w użyciu, a w przypadku zmiennych obserwowanych używa standardowej notacji wzorów R dla modeli. Poza tym definiujemy interesujące nas efekty, które chcemy obliczyć za pomocą operatora :=. Określamy model w całości jako prosty ciąg znaków, a następnie używamy funkcji sem do przeprowadzenia analizy.

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

Widzimy to samo wyjście co poprzednio i możemy porównać nasz indirect parametr z ACME, który mieliśmy wcześniej, direct efekt jest porównywany z ADE, a total porównuje się z poprzednim całkowitym efektem. Wartości są zasadniczo takie same.

Zauważ również, że dane wyjściowe pokazują wartość ^2 R dla obu modeli. W przypadku job_seek widzimy, że powodem, dla którego nie znajdujemy wiele w sposobie mediacji, jest to, że zaangażowane zmienne nie wyjaśniają żadnej zmienności mediatora na początku. Wstępne dochodzenie oszczędziłoby nam kłopotów w tym przypadku.

Pros

• Można obsługiwać wiele mediatorów
• Można obsługiwać wiele 'zabiegów’
• Można obsługiwać wiele wyników
• Można używać zmiennych ukrytych
• Jakieś wsparcie wielopoziomowe
• Można robić moderowane pośrednictwo i moderowaną moderację (choć nie dla zmiennych ukrytych). (choć nie dla zmiennych ukrytych)

Ograniczenia

• Wymaga dodatkowego kodowania, aby oszacować efekt pośredni
• Pojedyncze efekty losowe
• Mimo, że modele mogą zawierać zmienne binarne lub porządkowe dla mediatora/wyników, nie ma prostego sposobu na obliczenie efektu pośredniego w taki sposób, jak w przypadku pakietu mediacyjnego w tych warunkach.

piecewiseSEM

Pakiet piecewiseSEM działa bardzo podobnie do pakietu mediacji. W porównaniu z pakietem mediation zaletą pakietu piecewiseSEM jest to, że może on obsługiwać dodatkowe typy modeli, a także zapewniać dodatkowe dane wyjściowe (np. wyniki znormalizowane), dodatkowe opcje (np. multigrupa, skorelowane resztki) i wizualizację 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

Możemy wykorzystać jego możliwości plotowania, aby stworzyć szybką wizualizację modelu.

plot(mediation_result)

Niestety, obecnie nie ma automatycznego sposobu obliczania efektów pośrednich, więc trzeba by ręcznie bootstrapować wyniki.

Pros

• Standardowe modele R i składnia
• Wiele typów modeli zarówno dla mediatora, jak i wyniku
• Pewne wyniki w stylu SEM (np. dopasowanie, standaryzowane współczynniki, AIC)
• Szybki wykres wyników
• Możliwość obsługi wielu mediatorów, „zabiegów i wynikami

Ograniczenia

• Nie oblicza automatycznie efektów pośrednich
• Brak możliwości stosowania zmiennych ukrytych

psych

Pakiet psych wykorzystuje fakt, że w przypadku standardowego modelu liniowego można uzyskać wyniki za pomocą odpowiednich modeli regresji opartych tylko na macierzach kowariancji. Jest on bardzo podobny do lavaan, choć używa zwykłego podejścia najmniejszych kwadratów w przeciwieństwie do maksymalnego prawdopodobieństwa. Fajną rzeczą jest tutaj składnia, która pozwala skupić się tylko na interesującym nas efekcie lub uwzględnić wszystko, co jest fajne, jeśli interesowały nas również efekty pośrednie dla trudności ekonomicznych, wieku i płci.

Dla tego demo użyjemy wyczyszczonej wersji używając -, zamiast +, dla efektów nie związanych z leczeniem. Oznacza to tylko, że są one uwzględnione w modelach, ale wyniki ich dotyczące nie są pokazywane. Mediator jest identyfikowany za pomocą (). Kolejnym bonusem jest szybki wykres wyników, pokazujący różnicę między nieskorygowanymi i skorygowanymi efektami bezpośrednimi oraz odpowiedni przedział bootstrapowany.

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

Te same wyniki, inne opakowanie, ale prawdopodobnie najłatwiejsza droga, ponieważ wymagała tylko jednego wywołania funkcji. Pakiet psych obsługuje również wielu mediatorów i wyniki jako bonus.

Pros

• Najprostsza składnia, w zasadzie model jednowierszowy
• Szybki wykres wyników
• Możliwość obsługi wielu mediatorów, 'zabiegów’, i wyniki
• Możliwość 'moderowanej’ mediacji

Ograniczenia

• Ograniczony do standardowego modelu liniowego (lm)
• Użycie MASS

brms

Do naszego następnego dema dochodzimy do tego, co uważam za najpotężniejszy pakiet, brms. Nazwa oznacza Bayesian Regression Modeling with Stan, a Stan jest potężnym probabilistycznym językiem programowania do analizy bayesowskiej. Nie będę zagłębiał się w szczegóły analizy bayesowskiej, ale zapraszam do zapoznania się z moim dokumentem, który to robi.

Zazwyczaj postępujemy tak, jak wcześniej, określając model mediatora i model wyniku. brms nie robi nic specjalnego dla analizy mediacji, ale jego funkcja hipotezy może pozwolić nam na przetestowanie podejścia product-of-paths. Co więcej, pakiet sjstats zasadniczo dostarczy wyniki w ten sam sposób, w jaki robi to za nas pakiet mediatorów, a pakiet mediatorów jest w zasadzie próbą rozwiązania bayesowskiego przy użyciu metod frequentystycznych. Gdybyśmy mieli różne rozkłady dla wyniku i mediatora, mielibyśmy stosunkowo łatwy czas na uzyskanie tych średnich wartości predykcji i ich różnic, jako że podejścia bayesowskie zawsze myślą o rozkładach predykcyjnych posterior. W każdym razie, oto kod.

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

W danych wyjściowych, wszystko z jobseek_* jest wynikiem dla modelu mediatora, podczas gdy depress2_* jest dla wyniku. W tym momencie mamy tę samą starą historię, ale dzięki podejściu bayesowskiemu mamy więcej zabawnych rzeczy do oglądania. Na przykład, możemy zobaczyć, że tak naprawdę nie uchwyciliśmy dobrze skośności wyniku depresji. Nasze przewidywane wartości vs. obserwowane nie do końca się zgadzają. Jesteśmy trochę lepsi dla mediatora, ale być może nadal trochę za wysoko z niektórymi z naszych przewidywań opartych na modelu.

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

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

Pros

• Prosta składnia
• Niezwykle potężne-. Modele są w większości ograniczone do naszej wyobraźni
• Podstawowo robi to, co pakiet mediation aproksymuje
• Wszystkie zalety wnioskowania bayesowskiego: diagnostyka, posterior predictive checks, porównanie modeli itp.

Limitations

• Slower to estimate
• ’By-hand’ calculations needed for going beyond the standard linear model, ale jest to już powszechne podejście z perspektywy bayesowskiej
• Wymagana pewna wygoda z podejściem bayesowskim

Większa złożoność

Niektóre z wymienionych pakietów mogą obsługiwać bardziej złożone modele lub zapewnić dodatkowe podejścia do badania efektów pośrednich.

Interakcje

Niektóre modele zawierają interakcje albo dla modelu pośredniczącego albo wyniku, i niestety jest to często określane jako moderacja pośrednicząca lub moderowana mediacja. Osobiście nie widzę korzyści w nadawaniu niejednoznacznych nazw temu, co w przeciwnym razie może być prostą koncepcją (jeśli wciąż nie jest tak prostym modelem), ale ten statek odpłynął dawno temu. Nie będę się zagłębiał w szczegóły, ale chodzi o to, że możesz mieć termin interakcji gdzieś w modelu, a interakcja może dotyczyć zmiennej leczenia, mediatora lub obu.

Wystarczy powiedzieć, że skoro używamy standardowych narzędzi modelowania, takich jak lm i jego rozszerzeń, włączenie interakcji jest trywialne dla wszystkich powyższych pakietów, ale podejście typu iloczyn ścieżek nie ma zastosowania (a*b != c').

Uogólnione modele liniowe

W niektórych przypadkach nasz mediator lub wynik może być binarny, zliczany lub coś, gdzie założenie rozkładu normalnego może nie być najlepszym pomysłem. Albo możemy chcieć zbadać nieliniowe zależności pomiędzy leczeniem/mediatorem/wynikiem. Albo możemy mieć dane, które mają skorelowane obserwacje, takie jak powtarzane pomiary lub podobne. Pakiet mediation szczyci się tym w szczególności, ale brms może zrobić wszystko, co może zrobić i więcej, choć może trzeba będzie wykonać trochę więcej pracy, aby faktycznie obliczyć wynik. lavaan może faktycznie wykonać ograniczony zestaw modeli dla zmiennych binarnych i porządkowych, ale uzyskanie odpowiedniego oszacowania pośredniego wymagałoby bardzo żmudnego podejścia ręcznego.

Brakujące dane

Często, gdy mamy do czynienia z takimi danymi, szczególnie w naukach społecznych, często brakuje danych na temat któregoś ze zmiennych. Czasami możemy je pominąć, jeśli nie jest ich zbyt wiele, ale w innych przypadkach będziemy chcieli coś z tym zrobić. Pakiety lavaan, psych i brms zapewniają jeden lub więcej sposobów radzenia sobie z tą sytuacją (np. wielokrotna imputacja).

Alternatywy

Obrazowaliśmy modele jako sieci węzłów, z łukami/krawędziami/ścieżkami łączącymi je. Nasza dyskusja obraca się wokół tego, co nazywamy grafami skierowanymi (Directed Acyclic Graphs – DAG), gdzie strzałki mogą iść tylko w jednym kierunku, bez pętli sprzężenia zwrotnego. Wynik każdej zmiennej wynikowej jest funkcją poprzedzających ją strzałek i jest warunkowo niezależny od innych. Niektóre modele teoretyczne mogą to rozluźnić, a inne mogą w ogóle nie mieć strzałek, tj. być nieukierunkowane, tak, że interesują nas tylko połączenia (np. w niektórych sieciach społecznych).

bnlearn

Pakiet bnlearn pozwala na badanie grafów skierowanych, częściowo skierowanych i nieukierunkowanych. Jeśli chodzi o DAG, możemy go użyć do zduplikowania modeli mediacji, o których już mówiliśmy. Fajną rzeczą jest to, że ten pakiet będzie efektywnie testował ścieżki do włączenia, zamiast je zakładać, ale nadal możemy narzucać teoretyczne ograniczenia w razie potrzeby. Nie tylko możemy wtedy szukać interesujących nas ścieżek w sposób zasadniczy z sieciami bayesowskimi i teorią grafu przyczynowego Pearla jako podstawą, ale również będziemy mieli narzędzia do dalszego unikania przepasowania poprzez walidację krzyżową.

Dla początkowego modelu, upewnimy się, że ścieżki istnieją pomiędzy leczeniem – mediatorem, leczeniem – wynikiem i mediatorem – wynikiem (biała lista). Będziemy wyłączać nonsensowne ścieżki, takie jak posiadanie strzałek do leczenia (które zostało losowo przypisane), płeć, trudności ekonomiczne i wiek (czarna lista). W przeciwnym razie zobaczymy, co sugerują dane.

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)

Widzimy na wykresie, że sprawy trochę się zmieniły. Na przykład, wiek ma teraz związek tylko z poczuciem własnej skuteczności w poszukiwaniu pracy, a płeć ma wpływ tylko na depresję.

Jeśli ograniczymy ścieżki, aby były tylko takie, jakie są w naszych poprzednich przykładach, otrzymamy te same wyniki.

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 

Główną rzeczą do zauważenia jest to, że szacowane parametry są takie same, jak te, które otrzymaliśmy w poprzednich pakietach. Jest to w zasadzie równoważne użyciu lavaan z domyślnym estymatorem maksymalnej wiarygodności.

Jeśli użyjemy leczenia i płci jako czynników, bnlearn wyprodukuje modele warunkowe, które są różne w zależności od wartości czynnika. Innymi słowy, można by mieć osobny model dla gdy treatment == 'treatment' i jeden dla gdy treatment == control. W naszym przypadku, byłoby to identyczne z dopuszczeniem wszystkiego do interakcji z leczeniem, np. lm( job_seek ~ treat * (econ_hard + sex + age)), i podobnie dla modelu depresji. To rozszerzyłoby się na potencjalnie każdą zmienną binarną (np. włączając płeć). Jeśli mediator jest zmienną binarną, to jest to prawdopodobnie to, co chcielibyśmy zrobić.

Python

DyrektorCSCAR Kerby Shedden poprowadził warsztaty Pythona na temat modeli mediacji, więc pokazuję tutaj implementację statsmodels. Jest ona zgodna z podejściem Imai i może być postrzegana jako wersja pakietu mediacyjnego w Pythonie. Dane wyjściowe są zasadniczo takie same jak te, które otrzymalibyśmy używając leczenia jako zmiennej czynnikowej, gdzie otrzymujemy oddzielne wyniki dla każdej kategorii leczenia. Jest to niepotrzebne dla naszego demo, więc możesz po prostu porównać 'średnie’ wyniki z poprzednimi wynikami pakietu mediacyjnego.

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

Na koniec, dostarczam opcję w Stata używając jej polecenia sem. Stata ułatwia uzyskanie efektów pośrednich w tym przykładzie, ale robi to dla każdego zmiennego, więc dane wyjściowe są nieco zbyt obszerne6. Osoby pracujące z programem Stata nie potrzebują oddzielnego pakietu SEM, aby uzyskać tego rodzaju wyniki.

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)

Podsumowanie

Modele z efektami pośrednimi wymagają starannego rozważenia teoretycznego w celu zastosowania ich do analizy danych. Jednakże, jeśli model jest odpowiedni do sytuacji danych, dość łatwo jest uzyskać wyniki z różnych pakietów w R. Co więcej, nie trzeba używać pakietu do modelowania równań strukturalnych, aby przeprowadzić analizę z efektami pośrednimi, a w rzeczywistości można zajść daleko używając standardowej składni R. Dla ściśle obserwowanych, tj. nie ukrytych, zmiennych, żadne narzędzie SEM nie jest konieczne, a nawet zalecane.

Porównanie pakietów podsumowane

Poniższa tabela może pomóc zdecydować, który pakiet wykorzystać do swoich potrzeb, biorąc pod uwagę rozważania teoretyczne.

1. Mam dużo bardziej szczegółowy dokument na temat SEM, w tym analizy mediacji.︎

2. Z jakiegoś powodu nie widać tego zbyt często w praktyce i można się zastanawiać, co zrobiono, aby dane nadawały się do takiego modelu, jeśli nie było to uzasadnione.︎

3. Imai udostępnia swoje artykuły na swojej stronie internetowej.︎

4. MASS został zastąpiony przez innych przez ponad dekadę w tym momencie, i przeważnie ma tendencję do zamulania twojego tidyverse i innych pakietów kiedy jest załadowany. To dobry pakiet (i był świetny w swoim czasie), ale jeśli chcesz go użyć w pakiecie, dobrze byłoby nie ładować go (lub innych pakietów) w środowisku tylko po to, by użyć funkcji lub dwóch. Najczęściej widzę go po prostu używanego do mvrnorm (wielowymiarowa dystrybucja normalna) i glm.nb, ale są inne pakiety z tą funkcjonalnością, które zapewniłyby dodatkowe korzyści, a nie maskują funkcje dplyr, które są jednymi z najczęściej używanych w społeczności R.︎

5. brms pracuje nad tym.︎

6. Opcje w kodzie są po to, aby tłumić/minimalizować to, co może być.︎