Free Press

Tabla de contenidos

Actualizado el 04 de agosto de 2020. El código se puede descargar aquí.

Introducción

En algunas situaciones podemos considerar el efecto indirecto de alguna variable sobre un resultado o desenlace. Por ejemplo, las malas condiciones de vida en el hogar durante la infancia pueden disminuir los resultados del aprendizaje en la escuela, lo que posteriormente tiene un efecto negativo en la calidad de vida posterior, por ejemplo, los ingresos a lo largo de la vida. En otro caso, podríamos considerar una única variable recogida en múltiples puntos temporales, de forma que exista un efecto de la variable en el momento 1 sobre el momento 2, y del momento 2 sobre el momento 3. La idea básica es algo así como:

\Ndecir que \N(\mathcal{A}) lleva a \N(\mathcal{B}\N), y luego \N(\mathcal{B}\Nlleva a \N(\mathcal{C}\N). Con los modelos de mediación, planteamos una variable intermedia entre la covariable normal \(\rightarrow\) y la trayectoria de los resultados que podríamos tener en el marco de la regresión estándar, y estos modelos nos permiten investigar tales comportamientos. En el caso anterior, la variable interviniente, o mediadora, es \mathcal{B}. A menudo se da el caso de que todavía podríamos tener un efecto directo de \(\mathcal{A}) sobre \(\mathcal{C}), pero al igual que con el modelo en general, esto estaría motivado teóricamente.

El análisis de la mediación es muy popular en las disciplinas de las ciencias sociales, aunque de ninguna manera se limita a ellas, y por lo general se lleva a cabo bajo la apariencia de modelado de ecuaciones estructurales (SEM), que a su vez es una orientación específica de los modelos gráficos más generalmente1. El modelo gráfico de un modelo de mediación podría tener el siguiente aspecto.

En este caso, a y b reflejan el camino indirecto del efecto de \(\mathrm{{X}\) sobre el resultado a través del mediador, mientras que c' es el efecto directo de \(\mathrm{X}\) sobre el resultado después de la ruta indirecta se ha eliminado (c sería el efecto antes de plantear el efecto indirecto, y c – c' es igual al efecto indirecto). El efecto total de \N(\Nmathrm{X}\Nes la combinación de los efectos indirectos y directos.

Debo señalar algunas cosas basadas en lo que veo en la consultoría a través de docenas de disciplinas. Para empezar, parece que muy pocas personas que piensan que necesitan un modelo de mediación realmente lo hacen. Por ejemplo, si usted no puede pensar en su modelo en términos temporales o físicos, de tal manera que \ (\mathrm{X}\) necesariamente conduce al mediador, que luego necesariamente conduce al resultado, es probable que no necesite un modelo de mediación. Si puede ver que las flechas van en cualquier dirección, de nuevo, probablemente no necesite dicho modelo. Además, si al describir su modelo, todo el mundo piensa que está hablando de una interacción (también conocida como moderación), es posible que no lo necesite. Y, por último, como se podría sospechar, si no hay una fuerte correlación entre las variables clave (\mathrm{X}\) y el mediador (camino a), y si no hay una fuerte correlación entre el mediador y el resultado(s) (camino b), probablemente no lo necesite. Aunque nada le impedirá hacer un análisis de mediación, sin estos requisitos previos, es casi seguro que tendrá un modelo débil y probablemente más confuso que el que tendría de otro modo.

En resumen, la mediación funciona mejor cuando hay conexiones causales fuertemente implicadas entre las variables. Incluso entonces, tal modelo debe ser comparado con un modelo más simple de no mediación2. En cualquier caso, hay algunas maneras muy fáciles de investigar tales modelos en R, y ese es el objetivo aquí, sólo para demostrar cómo se puede empezar.

Datos

Para la demostración de los modelos de mediación con los diferentes paquetes, utilizaremos el conjunto de datos de trabajos que viene con el paquete de mediación. Aquí está la descripción.

Estudio de Intervención de Búsqueda de Empleo (JOBS II). JOBS II es un experimento de campo aleatorizado que investiga la eficacia de una intervención de formación laboral en trabajadores desempleados. El programa está diseñado no sólo para aumentar el reempleo entre los desempleados, sino también para mejorar la salud mental de los solicitantes de empleo. En el experimento de campo JOBS II, 1.801 trabajadores desempleados recibieron un cuestionario de preselección y fueron asignados aleatoriamente a grupos de tratamiento y de control. Los del grupo de tratamiento participaron en talleres de capacitación laboral. En los talleres, los encuestados aprendieron técnicas de búsqueda de empleo y estrategias para hacer frente a los contratiempos en el proceso de búsqueda de empleo. Los del grupo de control recibieron un folleto con consejos para la búsqueda de empleo. En las entrevistas de seguimiento, las dos variables clave de los resultados fueron una medida continua de los síntomas depresivos basada en la lista de comprobación de síntomas de Hopkins, y una variable binaria, que representaba si el encuestado había conseguido un empleo.

Aquí se describen las variables de esta demostración. Hay otras disponibles con las que también podría querer jugar.

• econ_hard: Nivel de dificultad económica antes del tratamiento con valores de 1 a 5.
• sexo: Variable indicadora del sexo. 1 = mujer
• edad: Edad en años.
• educ: Factor con cinco categorías para el nivel educativo.
• job_seek: Una escala continua que mide el nivel de autoeficacia en la búsqueda de empleo con valores de 1 a 5. La variable mediadora.
• depress2: Medida de los síntomas depresivos post-tratamiento. La variable de resultado.
• tratar: Variable indicadora de si el participante fue seleccionado al azar para el programa de formación JOBS II. 1 = asignación a la participación.

data(jobs, package = 'mediation')

Modelo

Dados estos datos, los modelos para el mediador y el resultado son los siguientes:

\N-

Por lo tanto, esperamos que el entrenamiento en habilidades laborales tenga un efecto negativo sobre la depresión (es decir, un aumento del bienestar), pero al menos parte de esto se debería a un efecto positivo en la búsqueda de empleo.

Como modelo gráfico, podríamos representarlo sucintamente como sigue.

Paquetes

Vamos a ver los siguientes paquetes para demostrar cómo se puede llevar a cabo el análisis de mediación en R:

• mediación
• lavaan
• psych
• brms

Aunque estos serán el foco de atención, también anotaré algunas otras alternativas, incluyendo Python y Stata.

mediación

Empezaremos con el paquete de mediación, ya que básicamente no requiere más capacidad de programación para llevar a cabo que la que ya se posee al ejecutar modelos de regresión estándar en R. El paquete proporciona el efecto de mediación causal promedio, definido de la siguiente manera a partir del archivo de ayuda y los artículos de Imai3:

El efecto de mediación causal promedio (ACME) representa la diferencia esperada en el resultado potencial cuando el mediador tomó el valor que se realizaría bajo la condición de tratamiento en oposición a la condición de control, mientras que el estado de tratamiento en sí se mantiene constante.

Nótese cómo esta definición se centra en los valores esperados o predichos condicionales al valor del tratamiento. Esta noción de contrafactualidad, o cómo sería la observación bajo el escenario opuesto, tiene una larga historia en la modelización en este punto. Piénsese que, si uno está en el grupo de tratamiento, tendría un valor específico para el mediador y, dado esto, tendría un valor esperado específico para el resultado. Sin embargo, podríamos plantear la misma observación como si estuviéramos en el grupo de control, y evaluar el efecto sobre el resultado a través del mediador de la misma manera. Podemos evaluar los posibles resultados manteniendo constante el tratamiento. Pensar en los cambios de los resultados teniendo en cuenta el valor del mediador no hace ninguna suposición sobre el tipo de modelo. Así es como el paquete de mediación es capaz de incorporar diferentes modelos para el mediador frente al resultado. Por ejemplo, el mediador podría ser binario, requiriendo un modelo de regresión logística, mientras que el modelo de resultado podría ser un modelo de supervivencia.

En nuestro ejemplo, nos ceñiremos a los modelos lineales estándar (normales). Observe también que, aunque nuestro tratamiento es una variable binaria, esto se generaliza al caso continuo, en el que consideramos el resultado de un movimiento de una unidad en el «tratamiento». Para que el paquete de mediación funcione, simplemente ejecutamos nuestros respectivos modelos para el mediador y el resultado, y luego utilizamos la función mediate para obtener el resultado 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)

Los resultados anteriores demuestran que el ACME no es estadísticamente distinto de cero, o de ninguna mediación. El efecto directo medio es negativo pero tampoco es estadísticamente notable, como tampoco lo es el efecto total (efecto indirecto + directo). También se proporciona la «proporción mediada» soi disant, que es la relación entre el efecto indirecto y el total. Sin embargo, no se trata de una proporción, e incluso puede ser negativa, por lo que en su mayoría es un número sin sentido.

Pros

• Modelos y sintaxis estándar de R
• Múltiples tipos de modelos tanto para el mediador como para el resultado
• Proporciona múltiples resultados simultáneamente
• Buena documentación y artículos asociados disponibles de forma gratuita
• Puede hacer mediación «moderada»

Limitaciones

.
• Uso del MASS4
• Modelos de efectos aleatorios simples
• La funcionalidad puede ser limitada con algunas complejidades del modelo
• Sin capacidades de variables latentes

lavaan

En el caso específico de que tanto los modelos de mediación como los de resultado sean modelos lineales estándar con una distribución normal para la variable objetivo, el efecto indirecto equivale al producto de las trayectorias a y b del diagrama anterior. El efecto directo es el camino c'. Una comparación del efecto directo independiente, que podríamos llamar c, frente a este efecto directo estimado en el modelo de mediación c', es tal que c - c' = a*b. Lo que se mencionó anteriormente podría ser ahora más claro, si a o b son casi cero, entonces el efecto indirecto sólo puede ser casi cero, por lo que es prudente investigar tales relaciones de antemano.

Este enfoque del producto de las trayectorias (o de la diferencia de coeficientes) es el que adoptaremos con el paquete lavaan y, de hecho, en el momento de escribir este artículo, es nuestra única forma de hacerlo. lavaan está orientado específicamente al modelado de ecuaciones estructurales, como el análisis de factores, los modelos de crecimiento y los modelos de mediación como el que estamos llevando a cabo aquí, y es muy recomendable para este tipo de modelos. Aunque se limita al caso del modelo lineal estándar para evaluar la mediación, es la única de nuestras herramientas que puede incorporar fácilmente variables latentes5. Por ejemplo, podríamos tener nuestro resultado de depresión como una variable latente subyacente a los ítems individuales del cuestionario. Además, también podríamos incorporar múltiples mediadores y múltiples resultados.

Para mantener las cosas como hemos estado discutiendo, etiquetaré los caminos a, b y c' en lavaan de acuerdo con la forma en que han sido representados previamente. Por lo demás, lavaan es muy fácil de usar y, en el caso de las variables observadas, utiliza la notación estándar de fórmulas de R para los modelos. Además, definimos los efectos de interés que queremos calcular con el operador :=. Especificamos el modelo en su totalidad como una simple cadena de caracteres, y luego usamos la función sem para hacer el análisis.

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

Vemos la misma salida anterior y podemos comparar nuestro parámetro indirect con el ACME que teníamos antes, el efecto direct se compara con el ADE, y el total se compara con el efecto total anterior. Los valores son esencialmente los mismos.

Note también que la salida muestra el valor \(R^2\) para ambos modelos. En el caso de job_seek, podemos ver que la razón por la que no estamos encontrando mucho en el camino de la mediación es porque las covariables involucradas no explican ninguna variación en el mediador para empezar. La investigación preliminar nos habría ahorrado el trabajo en este caso.

Pros

• Puede manejar múltiples mediadores
• Puede manejar múltiples «tratamientos»
• Puede manejar múltiples resultados
• Puede usar variables latentes
• Algún apoyo multinivel
• Puede hacer mediación moderada y mediación moderación (aunque no para variables latentes)

Limitaciones

• Requiere codificación adicional para estimar el efecto indirecto
• Efectos aleatorios simples
• Aunque los modelos podrían incorporar variables binarias u ordinales para el mediador/resultado, no hay una forma directa de calcular el efecto indirecto a la manera del paquete de mediación en esos escenarios.

piecewiseSEM

El paquete piecewiseSEM funciona de forma muy similar al paquete de mediación. Lo bueno de esto en relación con el paquete de mediación es que piecewiseSEM puede manejar tipos adicionales de modelos, así como proporcionar una salida adicional (por ejemplo, resultados estandarizados), opciones adicionales (por ejemplo, multigrupo, residuos correlacionados), y la visualización del modelo.

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

Podemos utilizar sus capacidades de trazado para crear una visualización rápida del modelo.

plot(mediation_result)

Desgraciadamente, no hay una manera automática de calcular los efectos indirectos en la actualidad, por lo que uno tendría que bootstrap los resultados a mano.

Pros

• Modelos y sintaxis estándar de R
• Múltiples tipos de modelos tanto para el mediador como para el resultado
• Algunos resultados de estilo SEM (por ejemplo ajuste, coeficientes estandarizados, AIC)
• Tráfico rápido de resultados
• Puede manejar múltiples mediadores, «tratamientos y resultados

Limitaciones

• No calcula automáticamente los efectos indirectos
• No tiene capacidad para variables latentes

psych

El paquete psych aprovecha el hecho de que en el caso del modelo lineal estándar, uno puede obtener los resultados a través de los modelos de regresión apropiados basados únicamente en las matrices de covarianza. Es muy similar a lavaan, aunque utiliza un enfoque de mínimos cuadrados ordinarios en lugar de máxima verosimilitud. Lo bueno aquí es una sintaxis que le permite centrarse sólo en el efecto de interés, o incluir todo, lo que es agradable si usted estaba interesado en los efectos indirectos para las dificultades económicas, la edad y el sexo también.

Para esta demostración utilizaremos la versión depurada usando el -, en lugar de +, para los efectos no relacionados con el tratamiento. Esto sólo significa que se incluyen en los modelos, pero no se muestran los resultados relativos a ellos. El mediador se identifica con (). Otra ventaja es un gráfico rápido de los resultados, que muestra la diferencia entre los efectos directos no ajustados y ajustados, y el intervalo apropiado de bootstrap.

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

Los mismos resultados, diferente empaquetado, pero posiblemente la ruta más fácil hasta ahora, ya que sólo requiere una llamada a la función. El paquete psicológico también maneja múltiples mediadores y resultados como un bono.

Pros

• La sintaxis más fácil, básicamente un modelo de una línea
• Planificación rápida de los resultados
• Puede manejar múltiples mediadores, ‘tratamientos’ y resultados
• Puede hacer mediación ‘moderada’

Limitaciones

• Limitado al modelo lineal estándar (lm)
• Uso de MASS

brms

Para nuestra siguiente demostración llegamos a lo que creo que es el paquete más potente, brms. El nombre significa Bayesian Regression Modeling with Stan, y Stan es un potente lenguaje de programación probabilístico para el análisis bayesiano. No voy a entrar en detalles sobre el análisis bayesiano, pero siéntase libre de ver mi documento que lo hace.

En general, hacemos lo mismo que antes, especificando el modelo mediador y el modelo de resultado. brms no hace nada especial para el análisis de mediación, pero su función de hipótesis puede permitirnos probar el enfoque del producto de los caminos. Además, el paquete sjstats proporcionará esencialmente los resultados de la misma manera que lo hace el paquete de mediación para nosotros, y para el caso, el paquete de mediación es básicamente un intento de solución bayesiana utilizando métodos frecuentistas de todos modos. Si tuviéramos diferentes distribuciones para el resultado y el mediador, nos resultaría relativamente fácil obtener estos valores medios de predicción y sus diferencias, ya que los enfoques bayesianos siempre piensan en distribuciones predictivas posteriores. En cualquier caso, aquí está el código.

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

En la salida, cualquier cosa con jobseek_* es un resultado para el modelo mediador, mientras que depress2_* es para el resultado. En este punto tenemos la misma historia de siempre, pero con el enfoque bayesiano tenemos cosas más divertidas que ver. Por ejemplo, podemos ver que no estamos capturando bien la asimetría del resultado de la depresión. Nuestros valores predichos frente a los observados no coinciden del todo. Estamos un poco mejor para el mediador, pero tal vez todavía un poco alto con algunas de nuestras predicciones basadas en el modelo.

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

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

Pros

• Sintaxis directa
• Extremadamente potente- Los modelos se limitan principalmente a la imaginación de uno
• Básicamente hace lo que el paquete de mediación aproxima
• Todas las ventajas de la inferencia bayesiana: diagnósticos, comprobaciones predictivas posteriores, comparación de modelos, etc.

Limitaciones

• Más fácil de estimar
• Cálculos «a mano» necesarios para ir más allá del modelo lineal estándar, pero esto ya es un enfoque común desde la perspectiva bayesiana
• Se requiere cierta comodidad con el enfoque bayesiano

Más complejidad

Algunos de los paquetes mencionados pueden manejar modelos más complejos o proporcionar enfoques adicionales para investigar los efectos indirectos.

Interacciones

Algunos modelos implican interacciones, ya sea para el modelo de mediación o para el resultado, y desafortunadamente esto se conoce a menudo como moderación mediada o mediación moderada. Personalmente no veo la ventaja de dar nombres ambiguos a lo que de otro modo podría ser un concepto sencillo (si bien no es un modelo tan directo), pero ese barco zarpó hace tiempo. No voy a entrar en detalles, pero la idea es que puedes tener un término de interacción en alguna parte del modelo, y la interacción puede implicar la variable de tratamiento, el mediador, o ambos.

Basta con decir que, dado que estamos utilizando herramientas de modelado estándar como lm y sus extensiones, la incorporación de interacciones es trivial para todos los paquetes anteriores, pero el tipo de enfoque del producto de las trayectorias no se mantiene (a*b != c').

Modelos lineales generalizados

En algunos casos nuestro mediador o resultado puede ser binario, de recuento, o algo donde asumir una distribución normal podría no ser la mejor idea. O puede que queramos investigar relaciones no lineales entre el tratamiento/mediador/resultado. O puede que tengamos datos con observaciones correlacionadas, como mediciones repetidas o similares. El paquete de mediación se enorgullece de esto en particular, pero brms puede hacer todo lo que puede hacer y más, aunque es posible que tenga que hacer un poco más de trabajo para calcular realmente el resultado. lavaan realmente puede hacer un conjunto limitado de modelos para variables binarias y ordinales, pero obtener la estimación indirecta adecuada requeriría un enfoque muy tedioso a mano.

Datos que faltan

A menudo, cuando se trata de tales datos, especialmente en las ciencias sociales, los datos suelen faltar en cualquiera de las covariables. A veces podemos descartarlas si no son demasiadas, pero en otros casos querremos hacer algo al respecto. Los paquetes lavaan, psych y brms proporcionan una o más formas de abordar la situación (por ejemplo, la imputación múltiple).

Alternativas

Hemos estado representando los modelos como redes de nodos, con arcos/aristas/caminos que los conectan. Nuestra discusión gira en torno a lo que se denomina grafos acíclicos dirigidos (DAG), en los que las flechas sólo pueden ir en una dirección, sin bucles de retroalimentación. El resultado de cualquier variable de resultado es una función de las flechas que la preceden, y condicionalmente independiente de otras. Algunos modelos teóricos pueden relajar esto, y otros pueden no tener flechas en absoluto, es decir, son no dirigidos, de tal manera que estamos interesados en sólo las conexiones (por ejemplo, con algunas redes sociales).

bnlearn

El paquete bnlearn permite la investigación de gráficos dirigidos, parcialmente dirigidos y no dirigidos. En términos de DAGs, podemos utilizarlo para duplicar esencialmente los modelos de mediación que hemos estado discutiendo. Sin embargo, lo bueno es que este paquete probará eficientemente los caminos para su inclusión en lugar de asumirlos, pero todavía podemos imponer restricciones teóricas según sea necesario. No sólo podemos entonces buscar los caminos de interés de una manera principista con las redes bayesianas y la teoría de los gráficos causales de Pearl como base, también tendremos herramientas para evitar aún más el exceso de ajuste a través de la validación cruzada.

Para el modelo inicial, nos aseguraremos de que los caminos existen entre el tratamiento – mediador, el tratamiento – resultado, y el mediador – resultado (la lista blanca). Desautorizaremos los caminos sin sentido como tener flechas hacia el tratamiento (que fue asignado al azar), el sexo, las dificultades económicas y la edad (la lista negra). Por lo demás, veremos lo que sugieren los datos.

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)

Vemos en el gráfico que las cosas han cambiado un poco. Por ejemplo, la edad ahora sólo se relaciona con la autoeficacia en la búsqueda de empleo, y el sexo sólo tiene un efecto sobre la depresión.

Si restringimos las trayectorias para que sólo sean lo que son en nuestros ejemplos anteriores, obtendríamos los mismos resultados.

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 

Lo principal a tener en cuenta es que los parámetros estimados son iguales a lo que obtuvimos con los paquetes anteriores. Es esencialmente equivalente a usar lavaan con el estimador de máxima verosimilitud por defecto.

Si usamos el tratamiento y el sexo como factores, bnlearn producirá modelos condicionales que son diferentes dependiendo del valor del factor que se tome. Es decir, se tendría un modelo distinto para cuando treatment == 'treatment' y otro para cuando treatment == control. En nuestro caso, esto sería idéntico a permitir que todo interactúe con el tratamiento, por ejemplo, lm( job_seek ~ treat * (econ_hard + sex + age)), y lo mismo para el modelo de depresión. Esto se extendería potencialmente a cualquier variable binaria (por ejemplo, incluyendo el sexo). Si el mediador es una variable binaria, esto es probablemente lo que querríamos hacer.

Python

El director del CSCAR, Kerby Shedden, ha impartido un taller de Python sobre modelos de mediación, así que aquí muestro la implementación de statsmodels. Sigue el enfoque de Imai, por lo que puede considerarse como la versión en Python del paquete de mediación. El resultado es esencialmente el mismo que el que se obtendría utilizando el tratamiento como variable factorial, donde se obtienen resultados separados para cada categoría de tratamiento. Esto es innecesario para nuestra demostración, por lo que sólo puede comparar los resultados ‘promedio’ con los resultados del paquete de mediación anterior.

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

Por último, proporciono una opción en Stata utilizando su comando sem. Stata facilita la obtención de los efectos indirectos en este ejemplo, pero lo hace para cada covariable, por lo que la salida es un poco verbosa para decir lo menos6. Para aquellos que trabajan con Stata, no necesitan un paquete SEM separado para obtener este tipo de resultados.

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)

Resumen

Los modelos con efectos indirectos requieren una cuidadosa consideración teórica para emplearlos en el análisis de datos. Sin embargo, si el modelo es apropiado para su situación de datos, es bastante fácil obtener resultados de una variedad de paquetes en R. Además, uno no necesita usar un paquete de modelado de ecuaciones estructurales para llevar a cabo un análisis con efectos indirectos, y de hecho, uno puede llegar lejos usando la sintaxis estándar de R. Para las variables estrictamente observadas, es decir, no latentes, no es necesaria ninguna herramienta de SEM, o incluso se recomienda.

Comparación de paquetes resumida

La siguiente tabla puede ayudar a uno a decidir qué paquete utilizar para sus necesidades dadas sus consideraciones teóricas.

1. Tengo un documento mucho más detallado sobre el SEM, incluyendo el análisis de mediación.︎

2. Por alguna razón no se ve mucho esto en la práctica, y uno se pregunta qué se hizo para que los datos fueran susceptibles de un modelo así si no estaba justificado.︎

3. Imai pone sus artículos a disposición en su página web.︎

4. MASS ha sido sustituido por otros durante más de una década en este punto, y en su mayoría sólo tiende a ensuciar su tidyverse y otros paquetes cuando se carga. Es un buen paquete (y fue genial en su día), pero si quieres usarlo en un paquete, sería bueno no cargarlo (u otros paquetes) en el entorno sólo para usar una o dos funciones. Principalmente sólo veo que se utiliza para mvrnorm (distribución normal multivariante) y glm.nb, pero hay otros paquetes con esa funcionalidad que proporcionarían beneficios adicionales, y no enmascaran las funciones de dplyr, que son de las más utilizadas en la comunidad R.︎

5. brms está trabajando en ello.︎

6. Las opciones en el código están ahí para suprimir/minimizar lo que puede ser.︎