Free Press

Table des matières

Mise à jour le 04 août 2020. Le code peut être téléchargé ici.

Introduction

Dans certaines situations, nous pouvons considérer l’effet indirect d’une certaine variable sur un résultat ou une issue. Par exemple, de mauvaises conditions de vie à la maison dans l’enfance peuvent diminuer les résultats d’apprentissage à l’école, ce qui a ensuite un effet négatif sur la qualité de vie ultérieure, par exemple, les revenus de la vie. Dans un autre cas, nous pourrions considérer une variable unique recueillie à plusieurs points dans le temps, de sorte qu’il existe un effet de la variable au temps 1 sur le temps 2, et du temps 2 sur le temps 3. L’idée de base est quelque chose comme:

\

En d’autres termes, \(\mathcal{A}\) conduit à \(\mathcal{B}\), puis \(\mathcal{B}\) conduit à \(\mathcal{C}\). Avec les modèles de médiation, nous posons une variable intermédiaire entre la covariable normale \(\rightarrow\) et le chemin du résultat que nous pourrions avoir dans le cadre de la régression standard, et ces modèles nous permettent d’étudier de tels comportements. Dans l’exemple ci-dessus, la variable intervenante, ou médiateur, est \(\mathcal{B}\). Il est souvent le cas que nous pourrions encore avoir un effet direct de \(\mathcal{A}\) sur \(\mathcal{C}\), mais comme avec le modèle en général, ce serait théoriquement motivé.

L’analyse de médiation est très populaire dans les disciplines des sciences sociales, bien que nullement limité à ceux, et généralement menée sous le couvert de la modélisation des équations structurelles (SEM), qui est lui-même une orientation spécifique des modèles graphiques plus généralement1. Le modèle graphique d’un modèle de médiation pourrait ressembler à ce qui suit .

Dans ce cas, a et b reflètent le chemin indirect de l’effet de \(\mathrm{X}\) sur le résultat à travers le médiateur, tandis que c' est l’effet direct de \(\mathrm{X}\) sur le résultat après que le chemin indirect ait été supprimé (c serait l’effet avant de poser l’effet indirect, et c – c' est égal à l’effet indirect). L’effet total de \(\mathrm{X}\) est la combinaison des effets indirects et directs.

Je devrais noter quelques choses basées sur ce que je vois en consultant à travers des dizaines de disciplines. Pour commencer, il semble que très peu de personnes qui pensent avoir besoin d’un modèle de médiation en ont réellement besoin. Par exemple, si vous ne pouvez pas penser à votre modèle en termes temporels ou physiques, de sorte que \(\mathrm{X}\) mène nécessairement au médiateur, qui mène ensuite nécessairement au résultat, vous n’avez probablement pas besoin d’un modèle de médiation. Si vous pouvez voir les flèches aller dans un sens ou dans l’autre, là encore, vous n’avez probablement pas besoin d’un tel modèle. De même, si lors de la description de votre modèle, tout le monde pense que vous parlez d’une interaction (alias modération), vous n’en avez peut-être pas besoin. Enfin, comme on pourrait le soupçonner, s’il n’y a pas de forte corrélation entre les variables clés (\(\mathrm{X}\)) et le médiateur (path a), et s’il n’y a pas de forte corrélation entre le médiateur et le(s) résultat(s) (path b), vous n’avez probablement pas besoin de cela. Bien que rien ne vous empêchera de faire une analyse de médiation, sans ces conditions préalables, vous aurez presque certainement un modèle faible et probablement plus confus que vous auriez autrement.

En bref, la médiation fonctionne mieux lorsqu’il existe des liens causaux fortement implicites entre les variables. Même dans ce cas, un tel modèle devrait être comparé à un modèle plus simple de non médiation2. En tout cas, il y a quelques façons très faciles d’étudier de tels modèles dans R, et c’est le but ici, juste pour démontrer comment vous pouvez commencer.

Data

Pour la démonstration des modèles de médiation avec les différents paquets, nous utiliserons l’ensemble de données jobs qui est livré avec le paquet mediation. En voici la description.

Étude d’intervention sur la recherche d’emploi (JOBS II). JOBS II est une expérience de terrain randomisée qui étudie l’efficacité d’une intervention de formation professionnelle sur les chômeurs. Le programme est conçu non seulement pour augmenter le réemploi des chômeurs mais aussi pour améliorer la santé mentale des demandeurs d’emploi. Dans le cadre de l’expérience JOBS II, 1 801 chômeurs ont reçu un questionnaire de présélection et ont ensuite été répartis au hasard entre un groupe de traitement et un groupe de contrôle. Ceux du groupe de traitement ont participé à des ateliers de formation professionnelle. Dans ces ateliers, les personnes interrogées ont appris des techniques de recherche d’emploi et des stratégies d’adaptation pour faire face aux échecs dans le processus de recherche d’emploi. Les participants du groupe témoin ont reçu une brochure décrivant des conseils de recherche d’emploi. Lors des entretiens de suivi, les deux principales variables de résultat étaient une mesure continue des symptômes dépressifs basée sur la liste de contrôle des symptômes de Hopkins, et une variable binaire, représentant le fait que le répondant ait trouvé un emploi.

Voici une description des variables de cette démonstration. Il y en a d’autres disponibles avec lesquelles vous pourriez aussi vouloir jouer.

• econ_hard : Niveau de difficultés économiques avant traitement avec des valeurs de 1 à 5.
• sexe : variable indicatrice du sexe. 1 = femme
• âge : Âge en années.
• educ : Facteur avec cinq catégories pour le niveau d’éducation.
• job_seek : Une échelle continue mesurant le niveau d’auto-efficacité de la recherche d’emploi avec des valeurs de 1 à 5. La variable médiatrice.
• depress2 : Mesure des symptômes dépressifs après le traitement. La variable de résultat.
• treat : Variable indicatrice pour savoir si le participant a été sélectionné au hasard pour le programme de formation JOBS II. 1 = affectation à la participation.

data(jobs, package = 'mediation')

Modèle

D’après ces données, les modèles pour le médiateur et le résultat sont les suivants :

On s’attend donc à ce que la formation aux compétences professionnelles ait un effet négatif sur la dépression (c’est-à-dire une augmentation du bien-être), mais au moins une partie de cela serait due à un effet positif sur la recherche d’emploi.

Sous forme de modèle graphique, nous pourrions le représenter succinctement comme suit .

Packages

Nous allons examiner les packages suivants pour démontrer comment on peut effectuer une analyse de médiation dans R :

• médiation
• lavaan
• psych
• brms

Bien que ceux-ci seront au centre de l’attention, je noterai également d’autres alternatives, notamment Python et Stata.

médiation

Nous commencerons par le paquet de médiation, car il ne nécessite fondamentalement pas plus de capacité de programmation à mener que ce que l’on possède déjà en exécutant des modèles de régression standard dans R. Le package fournit l’effet de médiation causale moyenne, défini comme suit à partir du fichier d’aide et des articles d’Imai3:

L’effet de médiation causale moyenne (ACME) représente la différence attendue dans le résultat potentiel lorsque le médiateur a pris la valeur qui se réaliserait dans la condition de traitement par rapport à la condition de contrôle, tandis que le statut de traitement lui-même est maintenu constant.

Notez comment cette définition est axée sur les valeurs attendues ou prédites conditionnelles à la valeur de traitement. Cette notion de contrefactuel, ou ce à quoi ressemblerait l’observation dans le cadre opposé, a une longue histoire dans la modélisation à ce point. Pensez-y de cette façon : si une personne fait partie du groupe de traitement, elle aura une valeur spécifique pour le médiateur et, étant donné cela, elle aura une valeur attendue spécifique pour le résultat. Cependant, nous pourrions poser la même observation comme étant dans le groupe de contrôle, et évaluer l’effet sur le résultat par le biais du médiateur de la même manière. Nous pouvons évaluer les résultats potentiels tout en maintenant le traitement constant. Le fait de penser aux changements de résultats en fonction de la valeur du médiateur n’implique aucune hypothèse sur le type de modèle. C’est ainsi que le package de médiation est capable d’incorporer différents modèles pour le médiateur vs. l’outcome. Par exemple, le médiateur pourrait être binaire, nécessitant un modèle de régression logistique, tandis que le modèle de résultat pourrait être un modèle de survie.

Dans notre exemple, nous nous en tiendrons aux modèles linéaires standard (normaux). Notez également, que si notre traitement est une variable binaire, cela se généralise au cas continu, où nous considérons le résultat d’un mouvement d’une unité sur le  » traitement « . Pour que le package de médiation fonctionne, nous exécutons simplement nos modèles respectifs pour le médiateur et le résultat, puis nous utilisons la fonction mediate pour obtenir le résultat final.

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)

Les résultats ci-dessus démontrent que l’ACME n’est pas statistiquement distinct de zéro, ou de l’absence de médiation. L’effet direct moyen est négatif mais de même, il n’est pas statistiquement notable, pas plus que l’effet total (effet indirect + direct). Le soi disant ‘proportion médiée’, qui est le rapport entre l’effet indirect et le total, est également fourni. Cependant, ce n’est pas une proportion, et peut même être négatif, et donc est surtout un nombre sans signification.

Pros

• Modèles et syntaxe R standard
• Multiples types de modèles pour le médiateur et le résultat
• Fournit plusieurs résultats simultanément
• Bonne documentation et les articles associés sont librement disponibles
• Peut faire de la médiation « modérée »

Limitations

.
• Utilisation de MASS4
• Modèles simples à effets aléatoires
• Fonctionnalité peut-être limitée avec certaines complexités de modèle
• Pas de capacités de variables latentes

lavaan

Dans le cas spécifique où les deux modèles de médiation et de résultat sont des modèles linéaires standard avec une distribution normale pour la variable cible, l’effet indirect est équivalent au produit des chemins a et b dans le diagramme précédent. L’effet direct est le chemin c'. Une comparaison de l’effet direct autonome, que nous pourrions appeler c, contre cet effet direct estimé dans le modèle de médiation c', est telle que c - c' = a*b. Ce qui a été mentionné plus tôt pourrait maintenant être plus clair, si soit a ou b sont presque zéro, alors l’effet indirect ne peut être que presque zéro, il est donc prudent d’étudier de telles relations au préalable.

Cette approche par produit de chemins (ou différence de coefficients) est celle que nous adopterons avec le paquet lavaan, et en fait, à l’heure où nous écrivons ces lignes, c’est notre seule façon de procéder. lavaan est spécifiquement orienté vers la modélisation par équations structurelles, comme l’analyse factorielle, les modèles de croissance et les modèles de médiation comme nous le faisons ici, et il est fortement recommandé pour de tels modèles. Bien qu’il soit limité au cas standard du modèle linéaire pour évaluer la médiation, c’est le seul de nos outils qui peut incorporer facilement des variables latentes5. Par exemple, nous pourrions avoir notre résultat de dépression comme une variable latente sous-jacente aux items du questionnaire individuel. En outre, nous pourrions également incorporer des médiateurs multiples et des résultats multiples.

Pour garder les choses comme nous en avons discuté, je vais étiqueter les chemins a, b et c' dans lavaan selon la façon dont ils ont été représentés précédemment. Sinon, lavaan est très facile à utiliser, et dans le cas des variables observées, utilise la notation standard des formules R pour les modèles. Au-delà, nous définissons les effets d’intérêt que nous voulons calculer avec l’opérateur :=. Nous spécifions le modèle dans son intégralité sous la forme d’une simple chaîne de caractères, puis nous utilisons la fonction sem pour effectuer l’analyse.

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

Nous voyons la même sortie que précédemment et pouvons comparer notre paramètre indirect à l’ACME que nous avions auparavant, l’effet direct est comparé à l’ADE, et le total se compare à l’effet total précédent. Les valeurs sont essentiellement les mêmes.

Notez également que la sortie montre la valeur \(R^2\) pour les deux modèles. Dans le cas de job_seek, nous pouvons voir que la raison pour laquelle nous ne trouvons pas beaucoup de médiation est que les covariables impliquées n’expliquent aucune variation dans le médiateur pour commencer. Une enquête préliminaire nous aurait épargné la peine dans ce cas.

Pros

• Peut traiter de multiples médiateurs
• Peut traiter de multiples ‘traitements’
• Peut traiter de multiples résultats
• Peut utiliser des variables latentes
• Un certain soutien multi-niveaux
• Peut faire de la médiation modérée et de la modération médiée (mais pas pour les variables latentes). modération (mais pas pour les variables latentes)

Limitations

• Requiert un codage supplémentaire pour estimer l’effet indirect
• Effets aléatoires simples
• Alors que les modèles pourraient incorporer des variables binaires ou ordinales pour le médiateur/les résultats, il n’y a pas de moyen direct de calculer l’effet indirect à la manière du paquet de médiation dans ces contextes.

piecewiseSEM

Le package piecewiseSEM fonctionne de manière très similaire au package mediation. Ce qui est bien par rapport au package mediation, c’est que piecewiseSEM peut gérer des types de modèles supplémentaires, ainsi que fournir des sorties supplémentaires (par exemple, des résultats normalisés), des options supplémentaires (par exemple, multigroupe, résidus corrélés) et une visualisation du modèle.

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

Nous pouvons utiliser ses capacités de traçage pour créer une visualisation rapide du modèle.

plot(mediation_result)

Malheureusement, il n’y a pas de moyen automatique de calculer les effets indirects à l’heure actuelle, il faudrait donc bootstrap les résultats à la main.

Pros

• Modèles et syntaxe R standard
• Multiples types de modèles pour le médiateur et le résultat
• Certains résultats de style SEM (par ex. ajustement, coefficients standardisés, AIC)
• Tracé rapide des résultats
• Peut gérer plusieurs médiateurs, « traitements », et résultats

Limitations

• Ne calcule pas automatiquement les effets indirects
• Pas de capacités de variables latentes

psych

Le paquet psych tire avantage du fait que dans le cas du modèle linéaire standard, on peut obtenir les résultats via les modèles de régression appropriés basés sur les matrices de covariance seules. Il est très similaire à lavaan, bien qu’il utilise une approche des moindres carrés ordinaires par opposition au maximum de vraisemblance. Ce qui est bien ici, c’est une syntaxe qui vous permet de vous concentrer uniquement sur l’effet qui vous intéresse, ou de tout inclure, ce qui est bien si vous vous intéressez également aux effets indirects des difficultés économiques, de l’âge et du sexe.

Pour cette démo, nous utiliserons la version nettoyée en utilisant le -, au lieu du +, pour les effets de non-traitement. Cela signifie simplement qu’ils sont inclus avec les modèles, mais que les résultats ne sont pas montrés les concernant. Le médiateur est identifié par (). Un autre bonus est un graphique rapide des résultats, montrant la différence entre les effets directs non ajustés et ajustés, et l’intervalle bootstrapped approprié.

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

Mêmes résultats, emballage différent, mais probablement la voie la plus facile jusqu’à présent car elle ne nécessite qu’un appel de fonction. Le paquet psych gère également plusieurs médiateurs et résultats en prime.

Pros

• Syntaxe la plus simple, essentiellement un modèle d’une ligne
• Tracé rapide des résultats
• Peut gérer de multiples médiateurs, ‘traitements’, et résultats
• Peut faire de la médiation ‘modérée’

Limitations

• Limité au modèle linéaire standard (lm)
• Utilisation de MASS

brms

Pour notre prochaine démonstration, nous arrivons à ce que je pense être le paquet le plus puissant, brms. Le nom signifie Bayesian Regression Modeling with Stan, et Stan est un puissant langage de programmation probabiliste pour l’analyse bayésienne. Je n’entrerai pas dans les détails de l’analyse bayésienne, mais n’hésitez pas à consulter mon document qui le fait.

Nous faisons généralement comme précédemment, en spécifiant le modèle du médiateur et le modèle du résultat. brms ne fait rien de spécial pour l’analyse de la médiation, mais sa fonction d’hypothèse peut nous permettre de tester l’approche du produit des chemins. De plus, le paquet sjstats fournira essentiellement les résultats de la même manière que le paquet mediation le fait pour nous, et d’ailleurs, le paquet mediation est essentiellement une tentative de solution bayésienne utilisant des méthodes fréquentistes de toute façon. Si nous avions des distributions différentes pour le résultat et le médiateur, il nous serait relativement facile d’obtenir ces valeurs de prédiction moyennes et leurs différences, car les approches bayésiennes pensent toujours aux distributions prédictives postérieures. En tout cas, voici le code.

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

Dans la sortie, tout ce qui comporte jobseek_* est un résultat pour le modèle médiateur, tandis que depress2_* est pour le résultat. Nous avons la même vieille histoire à ce stade, mais avec l’approche bayésienne, nous avons des choses plus amusantes à regarder. Par exemple, nous pouvons voir que nous ne capturons pas bien l’asymétrie du résultat de la dépression. Nos valeurs prédites par rapport aux valeurs observées ne correspondent pas tout à fait. Nous sommes un peu mieux pour le médiateur, mais peut-être encore un peu élevé avec certaines de nos prédictions basées sur le modèle.

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

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

Pros

• Syntaxe simple
• Extrêmement puissant-. Les modèles sont principalement limités à l’imagination de chacun
• Fait essentiellement ce que le package de médiation approxime
• Tous les avantages de l’inférence bayésienne : diagnostics, contrôles prédictifs postérieurs, comparaison de modèles, etc.

Limitations

• Plus faible pour estimer
• Calculs « à la main » nécessaires pour aller au-delà du modèle linéaire standard, mais c’est déjà une approche commune du point de vue bayésien
• Un certain confort avec l’approche bayésienne est requis

Plus de complexité

Certains des paquets mentionnés peuvent gérer des modèles plus complexes ou fournir des approches supplémentaires pour étudier les effets indirects.

Interactions

Certains modèles impliquent des interactions soit pour le modèle de médiation, soit pour le résultat, et malheureusement on parle souvent de modération médiate ou de médiation modérée. Personnellement, je ne vois pas l’avantage de donner des noms ambigus à ce qui pourrait être un concept simple (même si le modèle n’est pas si simple), mais ce navire a navigué il y a longtemps. Je ne vais pas entrer dans les détails, mais l’idée est que vous pourriez avoir un terme d’interaction quelque part dans le modèle, et l’interaction pourrait impliquer la variable de traitement, le médiateur, ou les deux.

Il suffit de dire que, puisque nous utilisons des outils de modélisation standard comme lm et ses extensions, l’incorporation d’interactions est triviale pour tous les paquets ci-dessus, mais le type d’approche product-of-paths ne tient pas (a*b != c').

Modèles linéaires généralisés

Dans certains cas, notre médiateur ou notre résultat peut être binaire, compter, ou quelque chose où supposer une distribution normale pourrait ne pas être la meilleure idée. Ou nous pourrions vouloir étudier les relations non linéaires entre le traitement/médiateur/résultat. Ou nous pouvons avoir des données qui ont des observations corrélées comme des mesures répétées ou similaires. Le paquet de médiation se targue de cela en particulier, mais brms peut faire tout ce qu’il peut faire et plus, bien que vous pourriez avoir à faire un peu plus de travail pour calculer réellement le résultat. lavaan peut effectivement faire un ensemble limité de modèles pour les variables binaires et ordinales, mais obtenir l’estimation indirecte appropriée nécessiterait une approche manuelle très fastidieuse.

Données manquantes

Souvent, lorsque nous traitons avec de telles données, en particulier dans les sciences sociales, les données sont souvent manquantes sur l’une des covariables. Parfois, nous pouvons les laisser tomber s’il n’y en a pas trop, mais dans d’autres cas, nous voudrons faire quelque chose à ce sujet. Les paquets lavaan, psych, et brms fournissent une ou plusieurs façons de gérer la situation (par exemple, l’imputation multiple).

Alternatives

Nous avons représenté les modèles comme des réseaux de nœuds, avec des arcs/flèches/chemins les reliant. Notre discussion tourne autour de ce que l’on appelle des graphes acycliques dirigés (DAG) où les flèches ne peuvent aller que dans une seule direction, sans boucle de rétroaction. Le résultat de toute variable de résultat est une fonction des flèches qui la précèdent, et conditionnellement indépendant des autres. Certains modèles théoriques peuvent relaxer cela, et d’autres peuvent n’avoir aucune flèche du tout, c’est-à-dire être non dirigés, de telle sorte que nous sommes intéressés par les seules connexions (par exemple avec certains réseaux sociaux).

bnlearn

Le paquet bnlearn permet l’investigation de graphes dirigés, partiellement dirigés et non dirigés. En termes de DAGs, nous pouvons l’utiliser pour essentiellement dupliquer les modèles de médiation dont nous avons discuté. Ce qui est bien, c’est que ce paquet testera efficacement les chemins pour l’inclusion plutôt que de les supposer, mais nous pouvons toujours imposer des contraintes théoriques si nécessaire. Non seulement nous pouvons alors rechercher les chemins d’intérêt d’une manière fondée sur des principes avec les réseaux bayésiens et la théorie des graphes causaux de Pearl comme base, mais nous aurons également des outils pour éviter davantage le surajustement via la validation croisée.

Pour le modèle initial, nous nous assurerons que des chemins existent entre traitement – médiateur, traitement – résultat, et médiateur – résultat (la liste blanche). Nous interdirons les chemins absurdes comme les flèches vers le traitement (qui a été assigné au hasard), le sexe, les difficultés économiques et l’âge (la liste noire). Sinon, nous verrons ce que les données suggèrent.

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)

Nous voyons dans le graphique que les choses ont un peu changé. Par exemple, l’âge n’est maintenant lié qu’à l’auto-efficacité de la recherche d’emploi, et le sexe n’a un effet que sur la dépression.

Si nous restreignons les chemins pour qu’ils soient uniquement ce qu’ils sont dans nos exemples précédents, nous obtiendrons les mêmes résultats.

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 principale chose à noter est que les paramètres estimés sont égaux à ce que nous avons obtenu avec les paquets précédents. C’est essentiellement équivalent à l’utilisation de lavaan avec l’estimateur de vraisemblance maximale par défaut.

Si nous utilisons le traitement et le sexe comme facteurs, bnlearn produira des modèles conditionnels qui sont différents selon la valeur du facteur prise. En d’autres termes, on aurait un modèle distinct pour quand treatment == 'treatment' et un pour quand treatment == control. Dans notre cas, cela serait identique au fait de permettre à tout d’interagir avec le traitement, par exemple lm( job_seek ~ treat * (econ_hard + sex + age)), et de même pour le modèle de dépression. Cela s’étendrait potentiellement à toute variable binaire (par exemple, y compris le sexe). Si le médiateur est une variable binaire, c’est probablement ce que nous voudrions faire.

Python

Le directeur du CCSRA, Kerby Shedden, a donné un atelier Python sur les modèles de médiation, donc je montre la mise en œuvre de statsmodels ici. Elle suit l’approche Imai et peut donc être considérée comme la version Python du package mediation. La sortie est essentiellement la même que celle que vous obtiendriez en utilisant le traitement comme variable factorielle, où vous obtenez des résultats séparés pour chaque catégorie de traitement. Ceci n’est pas nécessaire pour notre démo, donc vous pouvez simplement comparer les résultats ‘moyens’ aux résultats du package de médiation précédent.

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

Enfin, je fournis une option dans Stata en utilisant sa commande sem. Stata permet d’obtenir facilement les effets indirects dans cet exemple, mais il le fait pour chaque covariable, de sorte que la sortie est un peu verbeuse pour le moins6. Pour ceux qui travaillent avec Stata, ils n’ont pas besoin d’un paquet SEM séparé pour obtenir ce genre de résultats.

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)

Résumé

Les modèles avec effets indirects nécessitent une considération théorique minutieuse à employer pour l’analyse des données. Cependant, si le modèle est approprié à votre situation de données, il est assez facile d’obtenir des résultats à partir d’une variété de paquets dans R. En outre, il n’est pas nécessaire d’utiliser un paquet de modélisation d’équation structurelle pour effectuer une analyse avec des effets indirects, et en fait, on peut aller loin en utilisant la syntaxe R standard. Pour des variables strictement observées, c’est-à-dire non latentes, aucun outil SEM n’est nécessaire, ni même recommandé.

Comparaison de packages résumée

Le tableau suivant peut aider une personne à décider quel package utiliser pour ses besoins compte tenu de ses considérations théoriques.

1. J’ai un document beaucoup plus détaillé sur le SEM, y compris l’analyse de médiation.︎

2. Pour une raison quelconque, on ne voit pas beaucoup cela dans la pratique, et on se demande ce qui a été fait pour que les données se prêtent à un tel modèle si ce n’était pas justifié.︎

3. Imai met ses articles à disposition sur son site web.︎

4. MASS a été supplanté par d’autres depuis plus d’une décennie à ce stade, et il a surtout tendance à salir votre tidyverse et d’autres paquets quand il est chargé. C’est un bon paquet (et il était génial à l’époque), mais si vous voulez l’utiliser dans un paquet, il serait bon de ne pas le charger (ou d’autres paquets) dans l’environnement juste pour utiliser une fonction ou deux. Je le vois surtout juste utilisé pour mvrnorm (distribution normale multivariée) et glm.nb, mais il y a d’autres paquets avec cette fonctionnalité qui apporteraient des avantages supplémentaires, et ne pas masquer les fonctions dplyr, qui sont parmi les plus utilisées dans la communauté R.︎

5. brms y travaille.︎

6. Les options dans le code sont là pour supprimer/minimiser ce qui peut l’être.︎