Free Press

Tabella dei contenuti

Aggiornato il 04 agosto 2020. Il codice può essere scaricato qui.

Introduzione

In alcune situazioni possiamo considerare l’effetto indiretto di qualche variabile su un risultato. Per esempio, le cattive condizioni di vita a casa durante l’infanzia possono diminuire i risultati di apprendimento a scuola, che successivamente hanno un effetto negativo sulla qualità della vita successiva, per esempio, i guadagni di reddito per tutta la vita. In un altro caso potremmo considerare una singola variabile raccolta in più punti temporali, in modo che esista un effetto della variabile al tempo 1 sul tempo 2, e al tempo 2 sul tempo 3. L’idea di base è qualcosa del tipo:

In altre parole, \(\mathcal{A}}) porta a \(\mathcal{B}), e poi \(\mathcal{B}) porta a \(\mathcal{C}). Con i modelli di mediazione, poniamo una variabile interveniente tra la covariata normale \(\rightarrow\) e il percorso di risultato che potremmo avere nell’impostazione di regressione standard, e questi modelli ci permettono di indagare tali comportamenti. Nel caso precedente, la variabile interveniente, o mediatore, è \(\mathcal{B}}). Spesso potremmo ancora avere un effetto diretto di \(\mathcal{A}}} su \(\mathcal{C}}), ma come per il modello in generale, questo sarebbe teoricamente motivato.

L’analisi di mediazione è molto popolare nelle discipline delle scienze sociali, anche se non è affatto limitata a queste, e di solito viene condotta sotto il pretesto della modellazione di equazioni strutturali (SEM), che di per sé è un orientamento specifico dei modelli grafici più in generale1. Il modello grafico di un modello di mediazione potrebbe apparire come il seguente.

In questo caso, a e b riflettono il percorso indiretto dell’effetto di \(\mathrm{X}) sul risultato attraverso il mediatore, mentre c' è l’effetto diretto di \(\mathrm{X}} sull’esito dopo che il percorso indiretto è stato rimosso (c sarebbe l’effetto prima di porre l’effetto indiretto, e c – c' è uguale all’effetto indiretto). L’effetto totale di \(\mathrm{X}) è la combinazione degli effetti indiretti e diretti.

Devo notare alcune cose basate su ciò che vedo nelle consulenze in decine di discipline. Per cominciare, sembra che poche persone che pensano di aver bisogno di un modello di mediazione lo facciano davvero. Per esempio, se non potete pensare al vostro modello in termini temporali o fisici, in modo tale che \(\mathrm{X}}) porti necessariamente al mediatore, che poi porta necessariamente al risultato, probabilmente non avete bisogno di un modello di mediazione. Se potete vedere le frecce che vanno in entrambe le direzioni, di nuovo, probabilmente non avete bisogno di un tale modello. Inoltre, se quando descrivete il vostro modello, tutti pensano che stiate parlando di un’interazione (ovvero di moderazione), potreste non averne bisogno. E infine, come si potrebbe sospettare, se non c’è una forte correlazione tra le variabili chiave (\(\mathrm{X}}) e il mediatore (percorso a), e se non c’è una forte correlazione tra il mediatore e l’esito(i) (percorso b), probabilmente non ne avete bisogno. Anche se nulla vi impedisce di fare l’analisi di mediazione, senza questi prerequisiti, avrete quasi certamente un modello debole e probabilmente più confuso di quello che avreste altrimenti.

In breve, la mediazione funziona meglio quando ci sono connessioni causali fortemente implicite tra le variabili. Anche allora, un tale modello dovrebbe essere confrontato con un modello più semplice senza mediazione2. In ogni caso, ci sono alcuni modi molto semplici per studiare tali modelli in R, e questo è l’obiettivo qui, solo per dimostrare come si può iniziare.

Data

Per la dimostrazione dei modelli di mediazione con i diversi pacchetti, useremo il set di dati jobs che viene fornito con il pacchetto mediation. Ecco la descrizione.

Job Search Intervention Study (JOBS II). JOBS II è un esperimento sul campo randomizzato che indaga l’efficacia di un intervento di formazione al lavoro su lavoratori disoccupati. Il programma è progettato non solo per aumentare la rioccupazione tra i disoccupati, ma anche per migliorare la salute mentale di chi cerca lavoro. Nell’esperimento sul campo JOBS II, 1.801 lavoratori disoccupati hanno ricevuto un questionario di pre-screening e sono stati poi assegnati a caso a gruppi di trattamento e di controllo. Quelli del gruppo di trattamento hanno partecipato a workshop di abilità lavorative. Nei workshop, gli intervistati hanno imparato le abilità di ricerca del lavoro e le strategie di coping per affrontare le battute d’arresto nel processo di ricerca del lavoro. Quelli nella condizione di controllo hanno ricevuto un opuscolo che descriveva i consigli per la ricerca di lavoro. Nelle interviste di follow-up, le due variabili di risultato chiave erano una misura continua dei sintomi depressivi basata sulla Hopkins Symptom Checklist, e una variabile binaria, che rappresenta se l’intervistato era diventato impiegato.

Qui c’è una descrizione delle variabili in questa dimostrazione. Ce ne sono altre disponibili con cui potreste anche voler giocare.

• econ_hard: Livello di difficoltà economica prima del trattamento con valori da 1 a 5.
• sex: Variabile indicatore per il sesso. 1 = femmina
• età: Età in anni.
• educ: Fattore con cinque categorie per il livello di istruzione.
• job_seek: Una scala continua che misura il livello di autoefficacia nella ricerca di lavoro con valori da 1 a 5. La variabile mediatrice.
• depressione2: Misura dei sintomi depressivi dopo il trattamento. La variabile di risultato.
• treat: variabile indicatore per se il partecipante è stato selezionato a caso per il programma di formazione JOBS II. 1 = assegnazione alla partecipazione.

data(jobs, package = 'mediation')

Modello

Dati questi dati i modelli per il mediatore e l’esito sono i seguenti:

Così ci aspettiamo che il training sulle competenze lavorative abbia un effetto negativo sulla depressione (cioè un aumento del benessere), ma almeno una parte di questo sarebbe dovuto a un effetto positivo sulla ricerca del lavoro.

Come modello grafico, potremmo rappresentarlo succintamente come segue.

Pacchetti

Guarderemo i seguenti pacchetti per dimostrare come si può condurre un’analisi di mediazione in R:

• mediation
• lavaan
• psych
• brms

Mentre questi saranno l’obiettivo, noterò anche alcune altre alternative, inclusi Python e Stata.

mediation

Inizieremo con il pacchetto mediation, poiché fondamentalmente non richiede più capacità di programmazione di quelle già possedute dall’esecuzione di modelli di regressione standard in R. Il pacchetto fornisce l’effetto medio di mediazione causale, definito come segue dal file di aiuto e dagli articoli di Imai3:

L’effetto medio di mediazione causale (ACME) rappresenta la differenza attesa nel risultato potenziale quando il mediatore assume il valore che si realizzerebbe nella condizione di trattamento rispetto alla condizione di controllo, mentre lo stato del trattamento stesso è tenuto costante.

Nota come questa definizione sia focalizzata sui valori attesi o predetti condizionati al valore del trattamento. Questa nozione di controfattuali, o come sarebbe l’osservazione sotto l’impostazione opposta, ha una lunga storia nella modellistica a questo punto. Pensatela in questo modo: se uno è nel gruppo di trattamento, avrebbe un valore specifico per il mediatore e, dato questo, avrebbe un valore atteso specifico per il risultato. Tuttavia, potremmo porre la stessa osservazione come se fosse nel gruppo di controllo, e valutare l’effetto sul risultato attraverso il mediatore allo stesso modo. Possiamo valutare i potenziali risultati mantenendo costante il trattamento. Pensare ai cambiamenti di risultato dato il valore del mediatore non fa alcuna assunzione sul tipo di modello. È così che il pacchetto di mediazione è in grado di incorporare diversi modelli per il mediatore rispetto al risultato. Per esempio, il mediatore potrebbe essere binario, richiedendo un modello di regressione logistica, mentre il modello di risultato potrebbe essere un modello di sopravvivenza.

Nel nostro esempio, ci atterremo ai modelli lineari standard (normali). Notate anche che mentre il nostro trattamento è una variabile binaria, questo si generalizza al caso continuo, dove consideriamo il risultato di un movimento di una unità sul ‘trattamento’. Affinché il pacchetto di mediazione funzioni, dobbiamo semplicemente eseguire i nostri rispettivi modelli per il mediatore e l’esito, quindi utilizzare la funzione mediate per ottenere il risultato finale.

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)

I risultati precedenti dimostrano che l’ACME non è statisticamente distinto da zero, o da nessuna mediazione. L’effetto diretto medio è negativo, ma allo stesso modo non è statisticamente notevole, e nemmeno l’effetto totale (effetto indiretto + diretto). Viene anche fornito il soi disant “proporzione mediata”, che è il rapporto tra l’effetto indiretto e il totale. Tuttavia questa non è una proporzione, e può anche essere negativa, e quindi è per lo più un numero senza senso.

Pros

• Modelli e sintassi R standard
• Tipi multipli di modelli sia per il mediatore che per il risultato
• Fornisce più risultati simultaneamente
• Buona documentazione e articoli associati sono disponibili gratuitamente
• Può fare mediazione ‘moderata’

Limitazioni

• Uso di MASS4
• Modelli semplici a effetti casuali
• Funzionalità forse limitata con alcune complessità del modello
• Nessuna capacità di variabili latenti

lavaan

Nel caso specifico in cui entrambi i modelli di mediazione e di risultato sono modelli lineari standard con una distribuzione normale per la variabile obiettivo, l’effetto indiretto è equivalente al prodotto dei percorsi a e b nel diagramma precedente. L’effetto diretto è il percorso c'. Un confronto tra l’effetto diretto indipendente, che potremmo chiamare c, e questo effetto diretto stimato nel modello di mediazione c', è tale che c - c' = a*b. Ciò che è stato menzionato prima potrebbe ora essere più chiaro, se o a o b sono quasi zero, allora l’effetto indiretto può solo essere quasi zero, quindi è prudente investigare tali relazioni in anticipo.

Questo prodotto di percorsi (o differenza di coefficienti) è l’approccio che adotteremo con il pacchetto lavaan, e infatti, al momento in cui scriviamo, questo è il nostro unico modo di procedere. lavaan è specificamente orientato alla modellazione di equazioni strutturali, come l’analisi dei fattori, i modelli di crescita e i modelli di mediazione come quelli che stiamo conducendo qui, ed è altamente consigliato per tali modelli. Mentre è limitato al caso del modello lineare standard per valutare la mediazione, è l’unico dei nostri strumenti che può incorporare prontamente le variabili latenti5. Per esempio, potremmo avere il nostro risultato della depressione come una variabile latente alla base delle singole voci del questionario. Inoltre, potremmo anche incorporare mediatori multipli e risultati multipli.

Per mantenere le cose come le abbiamo discusse, etichetterò i percorsi a, b e c' in lavaan secondo come sono stati rappresentati in precedenza. Altrimenti lavaan è molto facile da usare, e nel caso di variabili osservate, usa la notazione standard della formula di R per i modelli. Oltre a questo definiamo gli effetti di interesse che vogliamo calcolare con l’operatore :=. Specifichiamo il modello nella sua interezza come una semplice stringa di caratteri, poi usiamo la funzione sem per fare l’analisi.

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

Vediamo lo stesso output di prima e possiamo confrontare il nostro parametro indirect con l’ACME che avevamo prima, l’effetto direct viene confrontato con l’ADE, e il total si confronta con il precedente effetto totale. I valori sono essenzialmente gli stessi.

Nota anche che l’output mostra il valore \(R^2\) per entrambi i modelli. Nel caso di job_seek, possiamo vedere che la ragione per cui non stiamo trovando molta mediazione è perché le covariate coinvolte non spiegano alcuna variazione nel mediatore per cominciare. Un’indagine preliminare ci avrebbe risparmiato la fatica in questo caso.

Pro

• Può gestire mediatori multipli
• Può gestire “trattamenti” multipli
• Può gestire esiti multipli
• Può usare variabili latenti
• Alcuni supporti multilivello
• Può fare mediazione moderata e moderata (anche se non per le variabili latenti)

Limitazioni

• Richiede una codifica aggiuntiva per stimare l’effetto indiretto
• Solo effetti casuali
• Mentre i modelli potrebbero incorporare variabili binarie o ordinali per il mediatore/esiti, non c’è un modo semplice per calcolare l’effetto indiretto alla maniera del pacchetto mediazione in questi contesti.

piecewiseSEM

Il pacchetto piecewiseSEM funziona molto simile al pacchetto mediation. La cosa bella di questo rispetto al pacchetto mediazione è che piecewiseSEM può gestire tipi aggiuntivi di modelli, così come fornire output aggiuntivi (ad esempio risultati standardizzati), opzioni aggiuntive (ad esempio multigruppo, residui correlati) e visualizzazione del modello.

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

Possiamo usare le sue capacità di plotting per creare una rapida visualizzazione del modello.

plot(mediation_result)

Purtroppo, non c’è un modo automatico per calcolare gli effetti indiretti al momento, quindi si dovrebbe fare il bootstrap dei risultati a mano.

Pro

• Modelli e sintassi standard di R
• Tipi multipli di modelli sia per il mediatore che per il risultato
• Alcuni risultati in stile SEM (es. fit, coefficienti standardizzati, AIC)
• Partitura rapida dei risultati
• Può gestire mediatori multipli, ‘trattamenti’, e risultati

Limitazioni

• Non calcola automaticamente gli effetti indiretti
• Nessuna capacità di variabili latenti

psych

Il pacchetto psych sfrutta il fatto che nel caso del modello lineare standard, si possono ottenere i risultati attraverso i modelli di regressione appropriati basati solo sulle matrici di covarianza. È molto simile a lavaan, anche se utilizza un approccio ai minimi quadrati ordinari in opposizione alla massima verosimiglianza. La cosa bella qui è una sintassi che permette di concentrarsi solo sull’effetto di interesse, o di includere tutto, il che è bello se si è interessati anche agli effetti indiretti per le difficoltà economiche, l’età e il sesso.

Per questa demo useremo la versione ripulita usando il -, invece del +, per gli effetti non di trattamento. Questo significa semplicemente che sono inclusi nei modelli, ma i risultati non sono mostrati riguardo ad essi. Il mediatore è identificato con (). Un altro bonus è un rapido grafico dei risultati, che mostra la differenza tra gli effetti diretti non aggiustati e aggiustati, e l’appropriato intervallo bootstrapped.

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

Stessi risultati, confezione diversa, ma probabilmente il percorso più semplice, dato che richiede solo una chiamata alla funzione. Il pacchetto psych gestisce anche mediatori e risultati multipli come bonus.

Pro

• Sintassi più semplice, fondamentalmente un modello di una sola linea
• Trama veloce dei risultati
• Può gestire mediatori multipli, ‘trattamenti’, e risultati
• Può fare mediazione ‘moderata’

Limitazioni

• Limitato al modello lineare standard (lm)
• Uso di MASS

brms

Per la nostra prossima demo arriviamo a quello che credo sia il pacchetto più potente, brms. Il nome sta per Bayesian Regression Modeling with Stan, e Stan è un potente linguaggio di programmazione probabilistico per l’analisi bayesiana. Non entrerò nei dettagli dell’analisi bayesiana, ma sentitevi liberi di vedere il mio documento che lo fa.

In genere facciamo come abbiamo fatto prima, specificando il modello del mediatore e il modello dell’esito. brms non fa nulla di speciale per l’analisi di mediazione, ma la sua funzione di ipotesi può permetterci di testare l’approccio del prodotto dei percorsi. Inoltre, il pacchetto sjstats fornirà essenzialmente i risultati nello stesso modo in cui il pacchetto mediazione lo fa per noi, e per questo, il pacchetto mediazione è fondamentalmente un tentativo di soluzione bayesiana usando comunque metodi frequentisti. Se avessimo distribuzioni diverse per l’esito e il mediatore, avremmo un tempo relativamente facile per ottenere questi valori medi di previsione e le loro differenze, poiché gli approcci bayesiani pensano sempre alle distribuzioni predittive posteriori. In ogni caso, ecco il codice.

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

Nell’output, qualsiasi cosa con jobseek_* è un risultato per il modello mediatore, mentre depress2_* è per il risultato. A questo punto abbiamo la stessa vecchia storia, ma con l’approccio bayesiano abbiamo più cose divertenti da guardare. Per esempio, possiamo vedere che non stiamo effettivamente catturando bene l’asimmetria del risultato della depressione. I nostri valori previsti rispetto a quelli osservati non corrispondono del tutto. Siamo un po’ meglio per il mediatore, ma forse ancora un po’ alti con alcune delle nostre previsioni basate sul modello.

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

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

Pro

• Sintassi semplice
• Estremamente potente- I modelli sono per lo più limitati alla propria immaginazione
• Fa praticamente ciò che il pacchetto mediazione approssima
• Tutti i vantaggi dell’inferenza bayesiana: diagnostica, controlli predittivi posteriori, confronto dei modelli, ecc.

Limitazioni

• Più basso da stimare
• Calcoli ‘a mano’ necessari per andare oltre il modello lineare standard, ma questo è già un approccio comune dalla prospettiva bayesiana
• È richiesto un certo comfort con l’approccio bayesiano

Più complessità

Alcuni dei pacchetti menzionati possono gestire modelli più complessi o fornire approcci aggiuntivi per studiare gli effetti indiretti.

Interazioni

Alcuni modelli coinvolgono interazioni sia per il modello di mediazione che per il risultato, e sfortunatamente questo viene spesso chiamato moderazione mediata o mediazione moderata. Personalmente non vedo il vantaggio di dare nomi ambigui a quello che altrimenti potrebbe essere un concetto semplice (anche se il modello non è così semplice), ma quella nave è salpata molto tempo fa. Non entrerò nei dettagli, ma l’idea è che si potrebbe avere un termine di interazione da qualche parte nel modello, e l’interazione potrebbe coinvolgere la variabile di trattamento, il mediatore, o entrambi.

Basta dire che, poiché stiamo usando strumenti di modellazione standard come lm e le sue estensioni, incorporare le interazioni è banale per tutti i pacchetti di cui sopra, ma il tipo di approccio del prodotto dei percorsi non è valido (a*b != c').

Modelli lineari generalizzati

In alcuni casi il nostro mediatore o risultato può essere binario, conteggio, o qualcosa dove assumere una distribuzione normale potrebbe non essere la migliore idea. Oppure potremmo voler studiare relazioni non lineari tra trattamento/mediatore/risultato. Oppure potremmo avere dati che hanno osservazioni correlate come misure ripetute o simili. Il pacchetto mediazione si vanta in particolare di questo, ma brms può fare tutto ciò che può fare e anche di più, anche se potreste dover fare un po’ più di lavoro per calcolare effettivamente il risultato. lavaan può effettivamente fare un insieme limitato di modelli per variabili binarie e ordinali, ma ottenere la stima indiretta appropriata richiederebbe un approccio molto noioso a mano.

Dati mancanti

Spesso quando si ha a che fare con tali dati, specialmente nelle scienze sociali, i dati mancano spesso su qualsiasi covariata. A volte possiamo lasciarle perdere se non sono troppe, ma in altri casi vorremo fare qualcosa al riguardo. I pacchetti lavaan, psych, e brms forniscono uno o più modi per affrontare la situazione (per esempio l’imputazione multipla).

Alternative

Abbiamo rappresentato i modelli come reti di nodi, con archi/spigoli/percorsi che li collegano. La nostra discussione ruota intorno a quelli che sono chiamati Directed Acyclic Graphs (DAG) dove le frecce possono andare solo in una direzione senza loop di feedback. Il risultato di qualsiasi variabile di risultato è una funzione delle frecce che la precedono, e condizionatamente indipendente dalle altre. Alcuni modelli teorici possono rilassare questo, e altri possono non avere alcuna freccia, cioè non essere diretti, in modo tale che siamo interessati solo alle connessioni (ad esempio con alcune reti sociali).

bnlearn

Il pacchetto bnlearn permette di studiare grafi diretti, parzialmente diretti e non diretti. In termini di DAG, possiamo usarlo per duplicare essenzialmente i modelli di mediazione che abbiamo discusso. La cosa bella però è che questo pacchetto testerà efficientemente i percorsi per l’inclusione piuttosto che assumerli, ma possiamo ancora imporre vincoli teorici come necessario. Non solo potremo quindi cercare i percorsi di interesse in un modo basato su principi con le reti bayesiane e la teoria dei grafi causali di Pearl come base, ma avremo anche strumenti per evitare ulteriormente l’overfitting attraverso la cross-validazione.

Per il modello iniziale, ci assicureremo che esistano percorsi tra trattamento – mediatore, trattamento – risultato e mediatore – risultato (la whitelist). Non permetteremo percorsi insensati come avere frecce verso il trattamento (che è stato assegnato casualmente), il sesso, le difficoltà economiche e l’età (la lista nera). Altrimenti, vedremo cosa suggeriscono i dati.

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)

Vediamo nel grafico che le cose sono un po’ cambiate. Per esempio, l’età ora si riferisce solo all’autoefficacia nella ricerca di lavoro, e il sesso ha un effetto solo sulla depressione.

Se restringiamo i percorsi per essere solo quello che sono nei nostri esempi precedenti, otterremmo gli stessi risultati.

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 

La cosa principale da notare è che i parametri stimati sono gli stessi che abbiamo ottenuto con i pacchetti precedenti. E’ essenzialmente equivalente all’uso di lavaan con lo stimatore di massima verosimiglianza predefinito.

Se usiamo il trattamento e il sesso come fattori, bnlearn produrrà modelli condizionali che sono diversi a seconda del valore del fattore preso. In altre parole, si avrebbe un modello separato per quando treatment == 'treatment' e uno per quando treatment == control. Nel nostro caso, questo sarebbe identico a permettere che tutto interagisca con il trattamento, per esempio lm( job_seek ~ treat * (econ_hard + sex + age)), e lo stesso per il modello della depressione. Questo si estenderebbe potenzialmente a qualsiasi variabile binaria (per esempio includendo il sesso). Se il mediatore è una variabile binaria, questo è probabilmente ciò che vorremmo fare.

Python

Il direttore del CSCAR Kerby Shedden ha tenuto un workshop Python sui modelli di mediazione, quindi mostro qui l’implementazione di statsmodels. Segue l’approccio di Imai e quindi può essere vista come la versione Python del pacchetto di mediazione. L’output è essenzialmente lo stesso di quello che si avrebbe usando il trattamento come variabile fattore, dove si ottengono risultati separati per ogni categoria di trattamento. Questo non è necessario per la nostra dimostrazione, quindi potete semplicemente confrontare i risultati ‘medi’ con i risultati del precedente pacchetto di mediazione.

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

Infine, fornisco un’opzione in Stata usando il suo comando sem. Stata rende facile ottenere gli effetti indiretti in questo esempio, ma lo fa per ogni covariata, quindi l’output è un po’ prolisso a dir poco6. Per coloro che lavorano con Stata, non hanno bisogno di un pacchetto SEM separato per ottenere questo tipo di risultati.

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)

Sommario

I modelli con effetti indiretti richiedono un’attenta considerazione teorica da utilizzare per l’analisi dei dati. Tuttavia, se il modello è appropriato per la vostra situazione di dati, è abbastanza facile ottenere risultati da una varietà di pacchetti in R. Inoltre, non è necessario utilizzare un pacchetto di modellazione di equazioni strutturali per condurre un’analisi con effetti indiretti, e in effetti, si può arrivare lontano utilizzando la sintassi standard di R. Per le variabili strettamente osservate, cioè non latenti, nessuno strumento SEM è necessario o addirittura raccomandato.

Confronto dei pacchetti riassunto

La seguente tabella può aiutare a decidere quale pacchetto usare per i propri bisogni, date le proprie considerazioni teoriche.

1. Ho un documento molto più dettagliato sul SEM, inclusa l’analisi di mediazione.︎

2. Per qualche motivo non si vede molto nella pratica, e ci si chiede cosa sia stato fatto per rendere i dati adatti a un tale modello se non era garantito.︎

3. Imai rende i suoi articoli disponibili sul suo sito web.︎

4. MASS è stato sostituito da altri da oltre un decennio a questo punto, e per lo più tende a incasinare il tuo tidyverse e altri pacchetti quando viene caricato. È un buon pacchetto (ed era ottimo in passato), ma se volete usarlo in un pacchetto, sarebbe bene non caricarlo (o altri pacchetti) nell’ambiente solo per usare una o due funzioni. Per lo più lo vedo solo usato per mvrnorm (distribuzione normale multivariata) e glm.nb, ma ci sono altri pacchetti con quella funzionalità che fornirebbero ulteriori benefici, e non mascherare le funzioni dplyr, che sono tra le più usate nella comunità R.︎

5. brms ci sta lavorando.︎

6. Le opzioni nel codice sono lì per sopprimere/minimizzare ciò che può essere.︎