Mediation Analysis Examples Using R
Dr. Hairong Song
The code below demonstrates the uses of R packages for running basic mediation analysis and moderated mediation analysis.
Install packages
install.packages("mediation",repos = "http://cran.us.r-project.org")## package 'mediation' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
## C:\Users\song2170\AppData\Local\Temp\RtmpKAJjuw\downloaded_packagesinstall.packages("lavaan",repos = "http://cran.us.r-project.org")## package 'lavaan' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
## C:\Users\song2170\AppData\Local\Temp\RtmpKAJjuw\downloaded_packagesinstall.packages("psych",repos = "http://cran.us.r-project.org")## package 'psych' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
## C:\Users\song2170\AppData\Local\Temp\RtmpKAJjuw\downloaded_packagesinstall.packages("tidyverse",repos = "http://cran.us.r-project.org")## package 'tidyverse' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
## C:\Users\song2170\AppData\Local\Temp\RtmpKAJjuw\downloaded_packagesinstall.packages("knitr",repos = "http://cran.us.r-project.org")## package 'knitr' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
## C:\Users\song2170\AppData\Local\Temp\RtmpKAJjuw\downloaded_packagesLoad packages
library(mediation)
library(lavaan)
library(psych)
library(dplyr)
library(knitr)Read the data set framing in the mediation package
data("framing", package="mediation") # list all data sets in the package
framing = as.data.frame(framing) # read the data set "framing" in the package
framing %>%
headTail() %>%
kable()| cond | anx | age | educ | gender | income | emo | p_harm | tone | eth | treat | english | immigr | anti_info | cong_mesg | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 3 | a little anxious | 45 | high school | male | 13 | 7 | 6 | 0 | 1 | 0 | Oppose | 4 | 0 | 1 |
| 2 | 4 | somewhat anxious | 73 | bachelor’s degree or higher | male | 16 | 6 | 3 | 0 | 0 | 0 | Favor | 3 | 0 | 0 |
| 3 | 2 | a little anxious | 53 | some college | female | 3 | 8 | 7 | 1 | 0 | 0 | Strongly Oppose | 3 | 0 | 0 |
| 4 | 1 | not anxious at all | 45 | high school | male | 14 | 9 | 8 | 1 | 1 | 1 | Strongly Oppose | 4 | 0 | 1 |
| … | NA | NA | … | NA | NA | … | … | … | … | … | … | NA | … | … | … |
| 262 | 4 | somewhat anxious | 28 | bachelor’s degree or higher | female | 14 | 6 | 6 | 0 | 0 | 0 | Strongly Oppose | 1 | 0 | 1 |
| 263 | 2 | very anxious | 38 | bachelor’s degree or higher | female | 14 | 3 | 6 | 1 | 0 | 0 | Favor | 2 | 0 | 0 |
| 264 | 1 | somewhat anxious | 29 | bachelor’s degree or higher | male | 10 | 9 | 7 | 1 | 1 | 1 | Strongly Oppose | 4 | 0 | 1 |
| 265 | 4 | very anxious | 52 | some college | female | 1 | 3 | 2 | 0 | 0 | 0 | Oppose | 2 | 0 | 0 |
Mean center X and/or M and then compute the interaction terms
framing = framing %>% mutate(emo = emo-mean(emo)) # mean center continuous x and/or M;not to center Y
framing = framing %>% mutate(sex=recode(gender,"female"= 1,"male"= 0)) # recode gender to numerical
framing = framing %>% mutate(treat_emo = treat*emo,treat_sex = treat*sex,emo_sex = emo*sex) # create interaction terms1. Mediation Analysis Using mediation Package
Single-mediator model without covariates
set.seed (2024) # allow replication of the results
mediate = mediation::mediate #A mediate function is in both the "psych" and "mediation" packages. This allows us to use the correct mediate function from the "mediation" package
med.fit0 = lm(emo ~ treat, data = framing) # mediator model
out.fit0 = glm(cong_mesg ~ emo + treat, data = framing, family = binomial("probit")) # outcome model
med.out0 = mediate(med.fit0, out.fit0, treat = "treat", mediator = "emo",robustSE = TRUE, sims = 1000)
summary(med.out0)##
## Causal Mediation Analysis
##
## Quasi-Bayesian Confidence Intervals
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME (control) 0.0908 0.0433 0.15 <2e-16 ***
## ACME (treated) 0.0913 0.0433 0.15 <2e-16 ***
## ADE (control) 0.0115 -0.0985 0.13 0.85
## ADE (treated) 0.0120 -0.1098 0.14 0.85
## Total Effect 0.1028 -0.0274 0.23 0.12
## Prop. Mediated (control) 0.7994 -3.5839 6.00 0.12
## Prop. Mediated (treated) 0.8189 -3.3013 5.62 0.12
## ACME (average) 0.0910 0.0434 0.15 <2e-16 ***
## ADE (average) 0.0118 -0.1039 0.13 0.85
## Prop. Mediated (average) 0.8091 -3.4476 5.81 0.12
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sample Size Used: 265
##
##
## Simulations: 1000ACME: Average Causal Mediation Effect [total effect - direct effect]
ADE: Average Direct Effect [total effect - indirect effect]
Total Effect: Direct (ADE) + Indirect (ACME)
Prop. Mediated: Conceptually ACME / Total effect (This tells us how much of the total effect our indirect effect is “explaining”)
Single-mediator model with covariates
med.fit1 = lm(emo ~ treat + age + educ + sex + income, data = framing)
out.fit1 = glm(cong_mesg ~ emo + treat + age + educ + sex + income,
data = framing, family = binomial("probit"))
med.out1 = mediate(med.fit1, out.fit1, treat = "treat", mediator = "emo",
robustSE = TRUE, sims = 1000)
summary(med.out1)##
## Causal Mediation Analysis
##
## Quasi-Bayesian Confidence Intervals
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME (control) 0.0822 0.0303 0.14 <2e-16 ***
## ACME (treated) 0.0830 0.0317 0.14 <2e-16 ***
## ADE (control) 0.0163 -0.0982 0.13 0.80
## ADE (treated) 0.0171 -0.1067 0.14 0.80
## Total Effect 0.0993 -0.0279 0.24 0.14
## Prop. Mediated (control) 0.7317 -4.2983 7.98 0.14
## Prop. Mediated (treated) 0.7512 -3.8465 7.56 0.14
## ACME (average) 0.0826 0.0310 0.14 <2e-16 ***
## ADE (average) 0.0167 -0.1024 0.14 0.80
## Prop. Mediated (average) 0.7415 -4.1533 7.70 0.14
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sample Size Used: 265
##
##
## Simulations: 1000plot(med.out1)
Bootstrapping method for indirect effect
med.out11 = mediate(med.fit1, out.fit1, boot = TRUE, treat = "treat", mediator = "emo", sims = 1000) ## Running nonparametric bootstrapsummary(med.out11)##
## Causal Mediation Analysis
##
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME (control) 0.0856 0.0347 0.14 <2e-16 ***
## ACME (treated) 0.0867 0.0344 0.15 <2e-16 ***
## ADE (control) 0.0118 -0.1040 0.13 0.80
## ADE (treated) 0.0128 -0.1175 0.14 0.80
## Total Effect 0.0984 -0.0347 0.23 0.13
## Prop. Mediated (control) 0.8697 -3.4770 6.53 0.13
## Prop. Mediated (treated) 0.8805 -2.9485 6.19 0.13
## ACME (average) 0.0861 0.0349 0.14 <2e-16 ***
## ADE (average) 0.0123 -0.1101 0.13 0.80
## Prop. Mediated (average) 0.8751 -3.2136 6.33 0.13
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sample Size Used: 265
##
##
## Simulations: 1000Treatment and mediator interaction
It is possible that the mediation effect takes different values depending on the baseline treatment status. In such as situation, we can add an interaction term between the treatment and mediator to the outcome model.
med.fit2 = lm(emo ~ treat + age + educ + sex + income, data = framing)
# add the treatment and mediator interaction to the outcome model only
out.fit2 = glm(cong_mesg ~ treat + emo+ treat_emo + age + educ + sex + income,
data = framing, family = binomial("probit"))
med.out2 = mediate(med.fit2, out.fit2, treat = "treat", mediator = "emo",
robustSE = TRUE, sims = 1000)
summary(out.fit2)##
## Call:
## glm(formula = cong_mesg ~ treat + emo + treat_emo + age + educ +
## sex + income, family = binomial("probit"), data = framing)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.143865 0.501224 -2.282 0.02248 *
## treat -0.054659 0.229369 -0.238 0.81165
## emo 0.187288 0.041143 4.552 5.31e-06 ***
## treat_emo 0.078211 0.081841 0.956 0.33925
## age 0.005402 0.005540 0.975 0.32947
## educhigh school -0.182745 0.337420 -0.542 0.58810
## educsome college -0.581002 0.363739 -1.597 0.11020
## educbachelor's degree or higher -0.362888 0.370001 -0.981 0.32670
## sex -0.391459 0.175777 -2.227 0.02595 *
## income 0.081519 0.025700 3.172 0.00151 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 336.89 on 264 degrees of freedom
## Residual deviance: 278.39 on 255 degrees of freedom
## AIC: 298.39
##
## Number of Fisher Scoring iterations: 5summary(med.out2)##
## Causal Mediation Analysis
##
## Quasi-Bayesian Confidence Intervals
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME (control) 0.0759 0.0288 0.13 <2e-16 ***
## ACME (treated) 0.0741 0.0292 0.13 <2e-16 ***
## ADE (control) -0.0146 -0.1258 0.12 0.80
## ADE (treated) -0.0164 -0.1405 0.12 0.80
## Total Effect 0.0595 -0.0805 0.20 0.42
## Prop. Mediated (control) 0.8006 -7.9615 12.99 0.42
## Prop. Mediated (treated) 0.8190 -7.0620 12.15 0.42
## ACME (average) 0.0750 0.0287 0.13 <2e-16 ***
## ADE (average) -0.0155 -0.1335 0.12 0.80
## Prop. Mediated (average) 0.8098 -7.4063 12.65 0.42
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sample Size Used: 265
##
##
## Simulations: 1000plot(med.out2)
Moderated mediation
The following code estimates the mediation effect for each gender group.
# add interaction terms to both the mediator and the outcome models
med.fit3 = lm(emo ~ treat + sex + treat_sex + age + educ + income, data = framing)
out.fit3 = glm(cong_mesg ~ treat + emo + sex+ treat_sex + emo_sex + age + educ + income,
data = framing, family = binomial("probit"))
med.female = mediate(med.fit3, out.fit3, treat = "treat", mediator = "emo",
covariates = list(sex = 1), sims = 1000)
med.male = mediate(med.fit3, out.fit3, treat = "treat", mediator = "emo",
covariates = list(sex = 0), sims = 1000)
summary(med.fit3)##
## Call:
## lm(formula = emo ~ treat + sex + treat_sex + age + educ + income,
## data = framing)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.9090 -2.0655 -0.3542 1.7491 6.6581
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.655556 0.886806 1.867 0.06306 .
## treat 1.671838 0.535378 3.123 0.00200 **
## sex 0.226423 0.365054 0.620 0.53565
## treat_sex -0.608006 0.723071 -0.841 0.40121
## age 0.001259 0.010054 0.125 0.90042
## educhigh school -1.055228 0.634088 -1.664 0.09730 .
## educsome college -1.852526 0.666234 -2.781 0.00583 **
## educbachelor's degree or higher -3.062183 0.660378 -4.637 5.64e-06 ***
## income -0.033484 0.041748 -0.802 0.42326
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.533 on 256 degrees of freedom
## Multiple R-squared: 0.1918, Adjusted R-squared: 0.1666
## F-statistic: 7.596 on 8 and 256 DF, p-value: 4.12e-09summary(out.fit3)##
## Call:
## glm(formula = cong_mesg ~ treat + emo + sex + treat_sex + emo_sex +
## age + educ + income, family = binomial("probit"), data = framing)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.236617 0.507492 -2.437 0.014821 *
## treat 0.382214 0.296081 1.291 0.196735
## emo 0.205294 0.049761 4.126 3.7e-05 ***
## sex -0.220780 0.204660 -1.079 0.280693
## treat_sex -0.662937 0.411821 -1.610 0.107448
## emo_sex 0.001710 0.067793 0.025 0.979875
## age 0.004868 0.005605 0.869 0.385060
## educhigh school -0.205649 0.337480 -0.609 0.542280
## educsome college -0.639750 0.365862 -1.749 0.080358 .
## educbachelor's degree or higher -0.381445 0.369916 -1.031 0.302463
## income 0.086807 0.025738 3.373 0.000744 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 336.89 on 264 degrees of freedom
## Residual deviance: 276.55 on 254 degrees of freedom
## AIC: 298.55
##
## Number of Fisher Scoring iterations: 5summary(med.female)##
## Causal Mediation Analysis
##
## Quasi-Bayesian Confidence Intervals
##
## (Inference Conditional on the Covariate Values Specified in `covariates')
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME (control) 0.0957 0.0258 0.18 <2e-16 ***
## ACME (treated) 0.1029 0.0300 0.19 <2e-16 ***
## ADE (control) 0.1057 -0.0513 0.27 0.222
## ADE (treated) 0.1129 -0.0588 0.28 0.222
## Total Effect 0.2086 0.0324 0.37 0.012 *
## Prop. Mediated (control) 0.4446 0.1284 1.84 0.012 *
## Prop. Mediated (treated) 0.4967 0.1606 1.73 0.012 *
## ACME (average) 0.0993 0.0276 0.19 <2e-16 ***
## ADE (average) 0.1093 -0.0553 0.27 0.222
## Prop. Mediated (average) 0.4707 0.1445 1.76 0.012 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sample Size Used: 265
##
##
## Simulations: 1000summary(med.male)##
## Causal Mediation Analysis
##
## Quasi-Bayesian Confidence Intervals
##
## (Inference Conditional on the Covariate Values Specified in `covariates')
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME (control) 0.0999 0.0293 0.19 <2e-16 ***
## ACME (treated) 0.1048 0.0317 0.19 <2e-16 ***
## ADE (control) 0.1162 -0.0505 0.29 0.168
## ADE (treated) 0.1211 -0.0564 0.30 0.168
## Total Effect 0.2210 0.0429 0.38 0.026 *
## Prop. Mediated (control) 0.4351 0.1141 1.43 0.026 *
## Prop. Mediated (treated) 0.4649 0.1229 1.41 0.026 *
## ACME (average) 0.1023 0.0305 0.19 <2e-16 ***
## ADE (average) 0.1186 -0.0520 0.29 0.168
## Prop. Mediated (average) 0.4500 0.1186 1.42 0.026 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sample Size Used: 265
##
##
## Simulations: 1000The following step tests the statistical significance of the difference in the ACME and ADE between males and females.
med.init = mediate(med.fit3, out.fit3, treat = "treat", mediator = "emo", sims = 1000)
test.modmed(med.init, covariates.1= list(sex=1),covariates.2= list(sex=0))##
## Test of ACME(covariates.1) - ACME(covariates.2) = 0
##
## data: estimates from med.init
## ACME(covariates.1) - ACME(covariates.2) = -0.0047963, p-value = 0.918
## alternative hypothesis: true ACME(covariates.1) - ACME(covariates.2) is not equal to 0
## 95 percent confidence interval:
## -0.1118559 0.1048036
##
##
## Test of ADE(covariates.1) - ADE(covariates.2) = 0
##
## data: estimates from med.init
## ADE(covariates.1) - ADE(covariates.2) = -0.0067298, p-value = 0.966
## alternative hypothesis: true ADE(covariates.1) - ADE(covariates.2) is not equal to 0
## 95 percent confidence interval:
## -0.2318936 0.22000002. Mediation Analysis Using “lavaan” Package
Standard single-mediator model
model1 = 'cong_mesg ~ c*treat # direct effect
emo ~ a*treat # mediator
cong_mesg ~ b*emo
ab := a*b # indirect effect (a*b)
total := c + (a*b) # total effect
'
fit1 = sem(model1, data = framing)
summary(fit1, fit.measures=TRUE, standardized=TRUE,rsquare=TRUE) ## lavaan 0.6-19 ended normally after 1 iteration
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 5
##
## Number of observations 265
##
## Model Test User Model:
##
## Test statistic 0.000
## Degrees of freedom 0
##
## Model Test Baseline Model:
##
## Test statistic 54.185
## Degrees of freedom 3
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 1.000
## Tucker-Lewis Index (TLI) 1.000
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -795.365
## Loglikelihood unrestricted model (H1) -795.365
##
## Akaike (AIC) 1600.731
## Bayesian (BIC) 1618.629
## Sample-size adjusted Bayesian (SABIC) 1602.776
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.000
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.000
## P-value H_0: RMSEA <= 0.050 NA
## P-value H_0: RMSEA >= 0.080 NA
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.000
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## cong_mesg ~
## treat (c) 0.015 0.063 0.231 0.818 0.015 0.014
## emo ~
## treat (a) 1.480 0.379 3.906 0.000 1.480 0.233
## cong_mesg ~
## emo (b) 0.063 0.010 6.276 0.000 0.063 0.368
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .cong_mesg 0.191 0.017 11.511 0.000 0.191 0.862
## .emo 7.253 0.630 11.511 0.000 7.253 0.946
##
## R-Square:
## Estimate
## cong_mesg 0.138
## emo 0.054
##
## Defined Parameters:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## ab 0.093 0.028 3.316 0.001 0.093 0.086
## total 0.107 0.066 1.626 0.104 0.107 0.099Bootstrapping method for indirect effects
fit2 = sem(model1, data = framing,se ="bootstrap",bootstrap=5000)
parameterEstimates(fit2, ci=TRUE, level=0.95, boot.ci.type="perc", standardized=TRUE) ## lhs op rhs label est se z pvalue ci.lower ci.upper
## 1 cong_mesg ~ treat c 0.015 0.064 0.229 0.819 -0.111 0.138
## 2 emo ~ treat a 1.480 0.385 3.847 0.000 0.736 2.238
## 3 cong_mesg ~ emo b 0.063 0.010 6.491 0.000 0.044 0.081
## 4 cong_mesg ~~ cong_mesg 0.191 0.011 16.982 0.000 0.166 0.211
## 5 emo ~~ emo 7.253 0.480 15.121 0.000 6.280 8.141
## 6 treat ~~ treat 0.191 0.000 NA NA 0.191 0.191
## 7 ab := a*b ab 0.093 0.029 3.222 0.001 0.041 0.154
## 8 total := c+(a*b) total 0.107 0.068 1.571 0.116 -0.029 0.241
## std.lv std.all std.nox
## 1 0.015 0.014 0.031
## 2 1.480 0.233 0.534
## 3 0.063 0.368 0.368
## 4 0.191 0.862 0.862
## 5 7.253 0.946 0.946
## 6 0.191 1.000 0.191
## 7 0.093 0.086 0.197
## 8 0.107 0.099 0.228Treatment and mediator interaction
model_int = 'cong_mesg ~ c*treat # direct effect
emo ~ a*treat # mediator
cong_mesg ~ b*emo + d*emo:treat
ab := a*b # indirect effect (a*b)
total := c + (a*b) # total effect
'
fit_int = sem(model_int, data = framing)
summary(fit_int, fit.measures=TRUE, standardized=TRUE,rsquare=TRUE) ## lavaan 0.6-19 ended normally after 9 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 7
##
## Number of observations 265
##
## Model Test User Model:
##
## Test statistic 108.980
## Degrees of freedom 2
## P-value (Chi-square) 0.000
##
## Model Test Baseline Model:
##
## Test statistic 165.023
## Degrees of freedom 6
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.327
## Tucker-Lewis Index (TLI) -1.018
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -1267.155
## Loglikelihood unrestricted model (H1) -1212.665
##
## Akaike (AIC) 2548.310
## Bayesian (BIC) 2573.368
## Sample-size adjusted Bayesian (SABIC) 2551.174
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.449
## 90 Percent confidence interval - lower 0.380
## 90 Percent confidence interval - upper 0.523
## P-value H_0: RMSEA <= 0.050 0.000
## P-value H_0: RMSEA >= 0.080 1.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.209
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## cong_mesg ~
## treat (c) -0.008 0.063 -0.126 0.899 -0.008 -0.008
## emo ~
## treat (a) 1.480 0.379 3.906 0.000 1.480 0.233
## cong_mesg ~
## emo (b) 0.055 0.010 5.500 0.000 0.055 0.327
## emo:treat (d) 0.031 0.019 1.677 0.094 0.031 0.097
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .cong_mesg 0.190 0.016 11.511 0.000 0.190 0.885
## .emo 7.253 0.630 11.511 0.000 7.253 0.946
## emo:treat 2.075 0.180 11.511 0.000 2.075 1.000
##
## R-Square:
## Estimate
## cong_mesg 0.115
## emo 0.054
##
## Defined Parameters:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## ab 0.081 0.025 3.185 0.001 0.081 0.076
## total 0.073 0.065 1.127 0.260 0.073 0.069Moderated mediation analysis
We create arbitrary names for both gender groups and tell R that we are doing this through the ‘c’ function in R. For example, we label the ‘a’ parameter from our initial model as ‘c(ag1, ag2)’ – we simply added the ‘g1’ and ‘g2’ after ‘a’ to differentiate group 1 versus group 2.
model2 = 'cong_mesg ~ c(cg1,cg2)*treat # direct effect
emo ~ c(ag1,ag2)*treat # mediator
cong_mesg ~ c(bg1,bg2)*emo
ab_g1 := ag1*bg1 # indirect effect for group 1
ab_g2 := ag2*bg2 # indirect effect for group 2
total_g1 := cg1 + (ag1*bg1) # total effect
total_g2 := cg2 + (ag2*bg2)
'
fit2 = sem(model2, data = framing,group="sex")
summary(fit2, fit.measures=TRUE, standardized=TRUE,rsquare=TRUE) ## lavaan 0.6-19 ended normally after 1 iteration
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 14
##
## Number of observations per group:
## 0 126
## 1 139
##
## Model Test User Model:
##
## Test statistic 0.000
## Degrees of freedom 0
## Test statistic for each group:
## 0 0.000
## 1 0.000
##
## Model Test Baseline Model:
##
## Test statistic 58.344
## Degrees of freedom 6
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 1.000
## Tucker-Lewis Index (TLI) 1.000
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -790.855
## Loglikelihood unrestricted model (H1) -790.855
##
## Akaike (AIC) 1609.711
## Bayesian (BIC) 1659.827
## Sample-size adjusted Bayesian (SABIC) 1615.439
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.000
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.000
## P-value H_0: RMSEA <= 0.050 NA
## P-value H_0: RMSEA >= 0.080 NA
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.000
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
##
## Group 1 [0]:
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## cong_mesg ~
## treat (cg1) 0.108 0.096 1.123 0.261 0.108 0.094
## emo ~
## treat (ag1) 1.850 0.560 3.304 0.001 1.850 0.282
## cong_mesg ~
## emo (bg1) 0.068 0.015 4.607 0.000 0.068 0.387
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .cong_mesg 0.371 0.046 8.101 0.000 0.371 0.761
## .emo -0.557 0.273 -2.038 0.042 -0.557 -0.200
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .cong_mesg 0.195 0.025 7.937 0.000 0.195 0.820
## .emo 7.168 0.903 7.937 0.000 7.168 0.920
##
## R-Square:
## Estimate
## cong_mesg 0.180
## emo 0.080
##
##
## Group 2 [1]:
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## cong_mesg ~
## treat (cg2) -0.055 0.082 -0.672 0.502 -0.055 -0.054
## emo ~
## treat (ag2) 1.159 0.513 2.257 0.024 1.159 0.188
## cong_mesg ~
## emo (bg2) 0.058 0.013 4.369 0.000 0.058 0.354
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .cong_mesg 0.289 0.042 6.890 0.000 0.289 0.644
## .emo -0.211 0.268 -0.787 0.431 -0.211 -0.077
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .cong_mesg 0.177 0.021 8.337 0.000 0.177 0.879
## .emo 7.274 0.873 8.337 0.000 7.274 0.965
##
## R-Square:
## Estimate
## cong_mesg 0.121
## emo 0.035
##
## Defined Parameters:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## ab_g1 0.125 0.047 2.685 0.007 0.125 0.109
## ab_g2 0.067 0.033 2.006 0.045 0.067 0.067
## total_g1 0.233 0.100 2.337 0.019 0.233 0.204
## total_g2 0.012 0.085 0.143 0.886 0.012 0.012# test for overall group differences using Omnibus Wald-Test
all.constraints = 'ag1 == ag2
bg1 == bg2
cg1 == cg2'
lavTestWald(fit2, constraints = all.constraints) ## $stat
## [1] 3.186644
##
## $df
## [1] 3
##
## $p.value
## [1] 0.3637339
##
## $se
## [1] "standard"# test for specific group differences between paths using individual Wald-Test
lavTestWald(fit2, constraints = "ag1 == ag2") ## $stat
## [1] 0.8282337
##
## $df
## [1] 1
##
## $p.value
## [1] 0.3627838
##
## $se
## [1] "standard"Single-mediator model with covariates
model3 = 'cong_mesg ~ c*treat + age + sex + income # direct effect
emo ~ a*treat + age + sex + income # mediator
cong_mesg ~ b*emo
ab := a*b # indirect effect (a*b)
total := c + (a*b) # total effect
'
fit3 = sem(model3, data = framing)
summary(fit3, fit.measures=TRUE, standardized=TRUE,rsquare=TRUE) ## lavaan 0.6-19 ended normally after 1 iteration
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 11
##
## Number of observations 265
##
## Model Test User Model:
##
## Test statistic 0.000
## Degrees of freedom 0
##
## Model Test Baseline Model:
##
## Test statistic 76.092
## Degrees of freedom 9
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 1.000
## Tucker-Lewis Index (TLI) 1.000
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -784.411
## Loglikelihood unrestricted model (H1) -784.411
##
## Akaike (AIC) 1590.823
## Bayesian (BIC) 1630.200
## Sample-size adjusted Bayesian (SABIC) 1595.324
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.000
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.000
## P-value H_0: RMSEA <= 0.050 NA
## P-value H_0: RMSEA >= 0.080 NA
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.000
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## cong_mesg ~
## treat (c) 0.010 0.062 0.159 0.874 0.010 0.009
## age 0.002 0.002 1.225 0.221 0.002 0.069
## sex -0.122 0.053 -2.316 0.021 -0.122 -0.129
## income 0.020 0.007 2.936 0.003 0.020 0.166
## emo ~
## treat (a) 1.517 0.374 4.051 0.000 1.517 0.239
## age 0.009 0.010 0.859 0.390 0.009 0.051
## sex 0.099 0.329 0.300 0.764 0.099 0.018
## income -0.102 0.042 -2.422 0.015 -0.102 -0.144
## cong_mesg ~
## emo (b) 0.067 0.010 6.806 0.000 0.067 0.394
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .cong_mesg 0.181 0.016 11.511 0.000 0.181 0.816
## .emo 7.059 0.613 11.511 0.000 7.059 0.920
##
## R-Square:
## Estimate
## cong_mesg 0.184
## emo 0.080
##
## Defined Parameters:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## ab 0.102 0.029 3.481 0.000 0.102 0.094
## total 0.111 0.065 1.713 0.087 0.111 0.103