Free Press

data(jobs, package = 'mediation')

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)

Overstående resultater viser, at ACME ikke er statistisk adskilt fra nul, eller ingen mediation. Den gennemsnitlige direkte effekt er negativ, men ligeledes ikke statistisk bemærkelsesværdig, og det er den samlede effekt (indirekte + direkte effekt) heller ikke. Der er også angivet soi disant “proportion mediated”, som er forholdet mellem den indirekte effekt og den samlede effekt. Dette er imidlertid ikke en andel og kan endda være negativt, og er derfor for det meste et meningsløst tal.

Pros

• Standard R-modeller og syntaks
• Flere typer modeller for både mediator og resultat
• Giver flere resultater samtidigt
• God dokumentation og tilknyttede artikler er frit tilgængelige
• Kan lave ‘modereret’ mediation

Begrænsninger

• Brug af MASS4
• Enkle tilfældige effektmodeller
• Funktionalitet måske begrænset med nogle modelkompleksiteter
• Ingen latente variabelfunktioner

lavaan

I det specifikke tilfælde, hvor både medierings- og udfaldsmodeller er standard lineære modeller med en normalfordeling for målvariablen, svarer den indirekte effekt til produktet af a og b stierne i det foregående diagram. Den direkte effekt svarer til stien c'. En sammenligning af den standalone direkte effekt, som vi kan kalde c, mod denne estimerede direkte effekt i mediationsmodellen c', er sådan, at c - c' = a*b. Det, der blev nævnt tidligere, er nu måske mere klart, hvis enten a eller b er næsten nul, kan den indirekte effekt kun være næsten nul, så det er klogt at undersøge sådanne sammenhænge på forhånd.

Denne produkt-of-paths (eller forskel i koefficienter) tilgang er den vi vil benytte os af med lavaan-pakken, og faktisk er det i skrivende stund vores eneste måde at gribe det an på. lavaan er specielt rettet mod strukturel ligningsmodellering, såsom faktoranalyse, vækstmodeller og mediationsmodeller, som vi udfører her, og anbefales stærkt til sådanne modeller. Selv om det er begrænset til den lineære standardmodel til vurdering af mægling, er det det eneste af vores værktøjer, der let kan inkorporere latente variabler5. Vi kunne f.eks. have vores depressionsresultat som en latent variabel, der ligger til grund for de enkelte spørgeskemaelementer. Derudover kunne vi også inkorporere flere mediatorer og flere udfald.

For at holde tingene, som vi har diskuteret, vil jeg mærke a, b og c' stierne i lavaan i overensstemmelse med, hvordan de tidligere er blevet afbildet. Ellers er lavaan meget let at bruge, og i tilfælde af observerede variabler bruger lavaan standard R-formelnotation for modellerne. Ud over det definerer vi de effekter af interesse, som vi ønsker at beregne med :=-operatoren. Vi angiver modellen i sin helhed som en simpel tegnstreng og bruger derefter sem-funktionen til at foretage analysen.

library(lavaan)sem_model = ' job_seek ~ a*treat + econ_hard + sex + age depress2 ~ c*treat + econ_hard + sex + age + b*job_seek # direct effect direct := c # indirect effect indirect := a*b # total effect total := c + (a*b)'model_sem = sem(sem_model, data=jobs, se='boot', bootstrap=500)summary(model_sem, rsq=T) # compare with ACME in mediation

lavaan 0.6-6 ended normally after 25 iterations Estimator ML Optimization method NLMINB Number of free parameters 11 Number of observations 899 Model Test User Model: Test statistic 0.000 Degrees of freedom 0Parameter Estimates: Standard errors Bootstrap Number of requested bootstrap draws 500 Number of successful bootstrap draws 500Regressions: Estimate Std.Err z-value P(>|z|) job_seek ~ treat (a) 0.066 0.049 1.332 0.183 econ_hard 0.053 0.024 2.242 0.025 sex -0.008 0.047 -0.163 0.870 age 0.005 0.002 1.934 0.053 depress2 ~ treat (c) -0.040 0.044 -0.905 0.365 econ_hard 0.149 0.022 6.908 0.000 sex 0.107 0.038 2.831 0.005 age 0.001 0.002 0.332 0.740 job_seek (b) -0.240 0.030 -8.079 0.000Variances: Estimate Std.Err z-value P(>|z|) .job_seek 0.524 0.030 17.610 0.000 .depress2 0.373 0.022 17.178 0.000R-Square: Estimate job_seek 0.011 depress2 0.120Defined Parameters: Estimate Std.Err z-value P(>|z|) direct -0.040 0.045 -0.904 0.366 indirect -0.016 0.012 -1.324 0.185 total -0.056 0.046 -1.224 0.221

Vi ser det samme output som før og kan sammenligne vores indirect parameter med den ACME vi havde før, direct effekten sammenlignes med ADE, og total sammenlignes med den tidligere samlede effekt. Værdierne er stort set de samme.

Bemærk også, at outputtet viser \(R^2\)-værdien for begge modeller. I tilfældet job_seek kan vi se, at grunden til, at vi ikke finder meget i form af mediation, er, at de involverede kovariater ikke forklarer nogen variation i mediatoren til at begynde med. En forudgående undersøgelse ville have sparet os for besværet i dette tilfælde.

Pros

• Kan håndtere flere mediatorer
• Kan håndtere flere ‘behandlinger’
• Kan håndtere flere udfald
• Kan bruge latente variabler
• En vis støtte til flere niveauer
• Kan lave modereret mediation og medieret moderation (dog ikke for latente variabler)

Begrænsninger

• Kræver yderligere kodning for at estimere den indirekte effekt
• Enkle tilfældige effekter
• Mens modellerne kan inkorporere binære eller ordinale variabler for mediator/udfald, er der ingen direkte måde at beregne den indirekte effekt på som i mediationspakken i disse indstillinger.

piecewiseSEM

Den piecewiseSEM-pakke fungerer meget lig med mediations-pakken. Det gode ved denne i forhold til mediationspakken er, at piecewiseSEM kan håndtere yderligere typer af modeller, samt give yderligere output (f.eks. standardiserede resultater), yderligere indstillinger (f.eks. multigruppe, korrelerede residualer) og visualisering af modellen.

library(piecewiseSEM)model_mediator <- lm(job_seek ~ treat + econ_hard + sex + age, data = jobs)model_outcome <- lm(depress2 ~ treat + econ_hard + sex + age + job_seek, data = jobs)mediation_result <- psem(model_mediator, model_outcome, data = jobs)summary(mediation_result)

Structural Equation Model of mediation_result Call: job_seek ~ treat + econ_hard + sex + age depress2 ~ treat + econ_hard + sex + age + job_seek AIC BIC 26.000 88.417---Tests of directed separation: No independence claims present. Tests of directed separation not possible.Global goodness-of-fit: Fisher's C = 0 with P-value = 1 and on 0 degrees of freedom---Coefficients: Response Predictor Estimate Std.Error DF Crit.Value P.Value Std.Estimate job_seek treat 0.0656 0.0515 894 1.2748 0.2027 0.0425 job_seek econ_hard 0.0532 0.0246 894 2.1612 0.0309 0.0720 * job_seek sex -0.0076 0.0487 894 -0.1567 0.8755 -0.0052 job_seek age 0.0046 0.0023 894 1.9779 0.0482 0.0658 * depress2 treat -0.0403 0.0435 893 -0.9255 0.3550 -0.0291 depress2 econ_hard 0.1485 0.0208 893 7.1323 0.0000 0.2248 *** depress2 sex 0.1068 0.0411 893 2.5957 0.0096 0.0818 ** depress2 age 0.0006 0.0020 893 0.3306 0.7410 0.0104 depress2 job_seek -0.2400 0.0282 893 -8.4960 0.0000 -0.2682 *** Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05---Individual R-squared: Response method R.squared job_seek none 0.01 depress2 none 0.12

Vi kan bruge dens plotteringsmuligheder til at skabe en hurtig visualisering af modellen.

plot(mediation_result)

Der er desværre ingen automatisk måde at beregne de indirekte effekter på i øjeblikket, så man skal bootstrape resultaterne i hånden.

Pros

• Standard R-modeller og syntaks
• Flere typer modeller for både mediator og resultat
• En del resultater i SEM-stil (f.eks. fit, standardiserede koefficienter, AIC)
• Snapt plot af resultater
• Kan håndtere flere mediatorer, ‘behandlinger’, og resultater

Begrænsninger

• Beregner ikke automatisk indirekte effekter
• Ingen muligheder for latente variabler

psych

Pakken psych udnytter det faktum, at man i tilfældet med den lineære standardmodel kan få resultaterne via de passende regressionsmodeller baseret på kovariansmatricerne alene. Den minder meget om lavaan, om end den anvender en almindelig mindste kvadraters metode i modsætning til maximum likelihood. Det gode her er en syntaks, der gør det muligt kun at fokusere på den ønskede effekt eller at inkludere alt, hvilket er rart, hvis man også var interesseret i de indirekte effekter for økonomiske vanskeligheder, alder og køn.

I denne demo bruger vi den rensede version med - i stedet for + for de ikke-behandlingseffekter. Det betyder blot, at de er medtaget i modellerne, men at resultaterne ikke vises vedrørende dem. Mediatoren er identificeret med (). En anden bonus er et hurtigt plot af resultaterne, der viser forskellen mellem de ujusterede og justerede direkte effekter og det passende bootstrappede 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

Samme resultater, anden indpakning, men muligvis den nemmeste vej endnu, da den kun krævede ét funktionskald. Psyk-pakken håndterer også flere mediatorer og resultater som en bonus.

Pros

• Nemmeste syntaks, dybest set en model på én linje
• Snapt plot af resultater
• Kan håndtere flere mediatorer, ‘behandlinger’, og resultater
• Kan lave ‘modereret’ mediation

Begrænsninger

• Begrænset til standard lineær model (lm)
• Brug af MASS

brms

For vores næste demo kommer vi til det, som jeg mener er den mest kraftfulde pakke, brms. Navnet står for Bayesian Regression Modeling with Stan, og Stan er et kraftfuldt probabilistisk programmeringssprog til Bayesian-analyse. Jeg vil ikke gå i detaljer om bayesiansk analyse, men du er velkommen til at se mit dokument, der gør det.

Vi gør generelt som vi har gjort før, idet vi specificerer mediatormodellen og udfaldsmodellen. brms gør ikke noget specielt for mediationsanalyse, men dens hypotesefunktion kan give os mulighed for at teste produkt-af-stier-tilgangen. Desuden vil sjstats-pakken i det væsentlige levere resultaterne på samme måde som mediationspakken gør for os, og for den sags skyld er mediationspakken i bund og grund alligevel et forsøg på en bayesiansk løsning ved hjælp af frekventistiske metoder. Hvis vi havde forskellige fordelinger for udfaldet og mediatoren, ville vi have relativt let ved at få disse gennemsnitlige forudsigelsesværdier og deres forskelle, da bayesianske tilgange altid tænker på posteriore forudsigelsesfordelinger. Under alle omstændigheder er her koden:

library(brms)model_mediator <- bf(job_seek ~ treat + econ_hard + sex + age)model_outcome <- bf(depress2 ~ treat + job_seek + econ_hard + sex + age)med_result = brm( model_mediator + model_outcome + set_rescor(FALSE), data = jobs)save(med_result, file = 'data/mediation_brms.RData')

load('data/mediation_brms.RData')summary(med_result)

 Family: MV(gaussian, gaussian) Links: mu = identity; sigma = identity mu = identity; sigma = identity Formula: job_seek ~ treat + econ_hard + sex + age depress2 ~ treat + job_seek + econ_hard + sex + age Data: jobs (Number of observations: 899) Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1; total post-warmup samples = 4000Population-Level Effects: Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESSjobseek_Intercept 3.67 0.12 3.43 3.91 1.00 6699 3749depress2_Intercept 2.21 0.15 1.92 2.50 1.00 6174 3091jobseek_treat 0.07 0.05 -0.03 0.17 1.00 6322 2709jobseek_econ_hard 0.05 0.02 0.00 0.10 1.00 6266 2656jobseek_sex -0.01 0.05 -0.10 0.09 1.00 5741 2655jobseek_age 0.00 0.00 0.00 0.01 1.00 6539 2846depress2_treat -0.04 0.04 -0.12 0.04 1.00 5458 3102depress2_job_seek -0.24 0.03 -0.30 -0.18 1.00 5950 2938depress2_econ_hard 0.15 0.02 0.11 0.19 1.00 7543 3102depress2_sex 0.11 0.04 0.03 0.19 1.00 5599 2699depress2_age 0.00 0.00 -0.00 0.00 1.00 4555 2887Family Specific Parameters: Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESSsigma_jobseek 0.73 0.02 0.69 0.76 1.00 6639 3276sigma_depress2 0.61 0.01 0.59 0.64 1.00 6145 2987Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESSand Tail_ESS are effective sample size measures, and Rhat is the potentialscale reduction factor on split chains (at convergence, Rhat = 1).

# using brms we can calculate the indirect effect as follows# hypothesis(med_result, 'jobseek_treat*depress2_job_seek = 0')# sjstats provides similar printing as the mediation package# print(sjstats::mediation(med_result), digits=4)sjstats::mediation(med_result) %>% kable_df()

I outputtet er alt med jobseek_* et resultat for mediatormodellen, mens depress2_* er et resultat for resultatet. Vi har den samme gamle historie på dette punkt, men med den bayesianske tilgang har vi flere sjove ting at kigge på. For eksempel kan vi se, at vi faktisk ikke fanger skævheden i depressionens udfald godt. Vores forudsagte værdier i forhold til de observerede stemmer ikke helt overens. Vi er en smule bedre for mediatoren, men måske stadig lidt for høje med nogle af vores modelbaserede forudsigelser.

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

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

Pros

• Straightforward syntax
• Extremely powerful- Modellerne er for det meste begrænset til ens fantasi
• Gør stort set det, som mediationspakken tilnærmer sig
• Alle fordelene ved Bayesisk inferens: diagnostik, posterior predictive checks, modelsammenligning osv.

Begrænsninger

• Sværere at estimere
• ‘By-hand’-beregninger, der er nødvendige for at gå ud over den lineære standardmodel, men dette er allerede en almindelig tilgang fra det bayesianske perspektiv
• En vis komfort med den bayesianske tilgang er påkrævet

Mere kompleksitet

En del af de nævnte pakker kan håndtere mere komplekse modeller eller giver yderligere tilgange til at undersøge indirekte effekter.

Interaktioner

Nogle modeller involverer interaktioner enten for mediationsmodellen eller resultatet, og desværre omtales dette ofte som medieret moderation eller modereret mediation. Personligt kan jeg ikke se fordelen ved at give tvetydige navne til noget, der ellers kunne være et ligetil koncept (om end stadig ikke så ligetil model), men det skib er sejlet for længe siden. Jeg vil ikke gå i detaljer, men ideen er, at du måske har et interaktionsterm et sted i modellen, og at interaktionen kan involvere behandlingsvariablen, mediatoren eller begge dele.

Det er nok at sige, da vi bruger standardmodelleringsværktøjer som lm og udvidelser af det, er det trivielt at inkorporere interaktioner for alle de ovennævnte pakker, men produkt-of-paths-typen af tilgang holder ikke (a*b != c').

Generaliserede lineære modeller

I nogle tilfælde kan vores mediator eller resultat være binær, tælle eller noget, hvor det måske ikke er den bedste idé at antage en normalfordeling. Eller vi ønsker måske at undersøge ikke-lineære sammenhænge mellem behandling/mediator/udbytte. Eller vi har måske data med korrelerede observationer som f.eks. gentagne målinger eller lignende. Mediations-pakken er især stolt af dette, men brms kan gøre alt, hvad den kan, og mere til, selv om du måske skal gøre lidt mere arbejde for rent faktisk at beregne resultatet. lavaan kan faktisk lave et begrænset sæt modeller for binære og ordinale variabler, men at få det passende indirekte estimat ville kræve en meget kedelig tilgang med hånden.

Manglende data

Ofte når man har med sådanne data at gøre, især inden for samfundsvidenskab, mangler der ofte data om nogen af kovariaterne. Nogle gange kan vi droppe disse, hvis der ikke er for mange, men i andre tilfælde vil vi ønske at gøre noget ved det. Pakkerne lavaan, psych og brms giver en eller flere måder at håndtere situationen på (f.eks. multiple imputation).

Alternativer

Vi har afbildet modellerne som netværk af knudepunkter, med buer/kanter/stier, der forbinder dem. Vores diskussion drejer sig om det, der kaldes Directed Acyclic Graphs (DAG), hvor pilene kun kan gå i én retning uden feedback loops. Resultatet af enhver udfaldsvariabel er en funktion af de pilene, der går forud for den, og er betinget uafhængigt af andre. Nogle teoretiske modeller kan lempe dette, og andre kan have slet ingen pile, dvs. være udirigerede, således at vi kun er interesseret i forbindelserne (f.eks. med nogle sociale netværk).

bnlearn

Den bnlearn-pakke giver mulighed for undersøgelse af dirigerede, delvist dirigerede og udirigerede grafer. Med hensyn til DAG’er kan vi bruge den til stort set at duplikere de formidlingsmodeller, som vi har diskuteret. Det gode er dog, at denne pakke effektivt vil teste stier for inklusion i stedet for at antage dem, men vi kan stadig pålægge teoretiske begrænsninger efter behov. Ikke alene kan vi så søge efter de stier af interesse på en principiel måde med bayesianske netværk og Pearls kausalgrafteori som grundlag, vi vil også have værktøjer til yderligere at undgå overfitting via krydsvalidering.

For den indledende model vil vi sikre os, at der findes stier mellem behandling – mediator, behandling – resultat og mediator – resultat (whitelisten). Vi vil udelukke nonsens-pædagogiske stier som f.eks. at have pile til behandlingen (som blev tildelt tilfældigt), køn, økonomiske vanskeligheder og alder (sortlisten). Ellers vil vi se, hvad dataene antyder.

whitelist = data.frame( from = c('treat', 'treat', 'job_seek'), to = c('job_seek', 'depress2', 'depress2'))blacklist = expand.grid( from = colnames(mediation_result$model.y$model), to = c('treat', 'sex', 'age', 'econ_hard'))# For simpler output we'll use treatment and sex as numeric (explained later)library(dplyr)jobs_trim = jobs %>% select(depress2, treat, econ_hard, sex, age, job_seek) %>% mutate( treat = as.numeric(jobs$treat), sex = as.numeric(jobs$sex) )# extract path coefficients if desired# parameters = bn.fit(model, jobs_trim)# parameters$job_seek# parameters$econ_hard# parameters$depress2

library(bnlearn)model = gs(jobs_trim, whitelist = whitelist, blacklist = blacklist)plot(model)

Vi kan se i plottet, at tingene har ændret sig en smule. For eksempel hænger alder nu kun sammen med jobsøgende self-efficacy, og køn har kun en effekt på depression.

Hvis vi begrænser stierne til kun at være det, de er i vores tidligere eksempler, får vi de samme resultater.

library(bnlearn)whitelist = data.frame( from = c('treat', 'age', 'sex', 'econ_hard', 'treat', 'job_seek', 'age', 'sex', 'econ_hard'), to = c('job_seek', 'job_seek','job_seek','job_seek', 'depress2', 'depress2', 'depress2', 'depress2', 'depress2'))blacklist = expand.grid( from = colnames(mediation_result$model.y$model), to = c('treat', 'sex', 'age', 'econ_hard')) model = gs(jobs_trim, whitelist = whitelist, blacklist = blacklist)plot(model)

parameters = bn.fit(model, jobs_trim)parameters$depress2$coefficients

 (Intercept) treat econ_hard sex age job_seek 2.2076414333 -0.0402647000 0.1485433818 0.1068048699 0.0006488642 -0.2399549527 

parameters$job_seek$coefficients

 (Intercept) treat econ_hard sex age 3.670584908 0.065615003 0.053162413 -0.007637336 0.004586492 

Det vigtigste at bemærke er, at de estimerede parametre svarer til det samme, som vi fik med de tidligere pakker. Det svarer stort set til at bruge lavaan med standard maximum likelihood estimator.

Hvis vi bruger behandling og køn som faktorer, vil bnlearn producere betingede modeller, der er forskellige afhængig af den faktorværdi, der tages. Med andre ord vil man have en separat model for når treatment == 'treatment' og en for når treatment == control. I vores tilfælde ville dette være identisk med at lade alt interagere med behandling, f.eks. lm( job_seek ~ treat * (econ_hard + sex + age)), og det samme ville gælde for depressionsmodellen. Dette ville blive udvidet til potentielt enhver binær variabel (f.eks. inklusive køn). Hvis mediatoren er en binær variabel, er det sandsynligvis det, vi ønsker at gøre.

Python

CSCAR-direktør Kerby Shedden har givet en Python-workshop om mediationsmodeller, så jeg viser statsmodels-implementeringen her. Den følger Imai-tilgangen og kan således ses som Python-versionen af mediationspakken. Outputtet er stort set det samme som det, man ville få ved at bruge behandling som en faktorvariabel, hvor man får separate resultater for hver behandlingskategori. Dette er unødvendigt for vores demo, så du kan bare sammenligne de “gennemsnitlige” resultater med de tidligere resultater fra mediationspakken.

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

Sidst giver jeg en mulighed i Stata ved hjælp af dens sem-kommando. Stata gør det nemt at få de indirekte effekter i dette eksempel, men den gør det for hver enkelt kovariat, så outputtet er mildest talstærkt for at sige det mildt6. For dem, der arbejder med Stata, er det ikke nødvendigt med en separat SEM-pakke for at få denne slags resultater.

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)