Package semMediation

Keon-Woong Moon

2016-08-07

Package semMediation is an extension of package lavaan. The main functions of package semMediation are as follows:

  1. makeEquation() : support function to make a mediation equations easily for lavaan.
  2. mediationPlot() : visualize the mediation effects

Installation

You can install semMediation package via github.

install.packages("devtools")
devtools::install_github("cardiomoon/semMediation")

makeEquation()

Consider a classical mediation setup with four variables: Y1 is the dependent variable, X1 is the predictor and M1 and M2 are mediators. For illustration, we create a toy dataset containing these variables, and fit a path analysis model that includes the direct effect of X1 on Y1 and the indirect effect of X1 on Y1 via M1 or M2 and second indirect effect via M1 and M2.

set.seed(1234)
X1 <- rnorm(100)
X2 <- c(rnorm(50),rnorm(50))
M1 <- 0.5*X1 + 0.5*X2 +rnorm(100)
M2 <-0.6*X1 +0.4*X2+rnorm(100)
Y1 <- 0.3*M1 + 0.4*M2 + rnorm(100)
Y2 <- 0.1*M1 + 0.5*M2 + rnorm(100)

data <- data.frame(X1 = X1, X2 = X2, Y1 = Y1, Y2 = Y2, M1 = M1, M2 = M2)
str(data)
'data.frame':   100 obs. of  6 variables:
 $ X1: num  -1.207 0.277 1.084 -2.346 0.429 ...
 $ X2: num  0.415 -0.475 0.066 -0.502 -0.826 ...
 $ Y1: num  -1.65548 -0.17509 0.00588 -1.35581 0.4317 ...
 $ Y2: num  0.424 -1.653 1.035 -0.481 1.769 ...
 $ M1: num  0.089 0.598 0.761 -0.723 0.113 ...
 $ M2: num  -1.1384 -0.9767 0.4976 -0.5986 -0.0493 ...

To make a mediatino equation, you can use makeEquation() function.

require(lavaan)
require(semPlot)
require(semMediation)
require(ggplot2)

model=makeEquation(X="X1",M=c("M1","M2"),Y="Y1")
cat(model)
# Mediation Effect
Y1~b1*M1+b2*M2+c1*X1
M1~a1*X1
M2~a2*X1+d1*M1
ind1:=a1*b1
ind2:=a2*b2
secondInd1:=a1*d1*b2
total1:=c1+a1*b1+a2*b2+a1*d1*b2
fit=sem(model,data=data)
summary(fit)
lavaan (0.5-22) converged normally after  19 iterations

  Number of observations                           100

  Estimator                                         ML
  Minimum Function Test Statistic                0.000
  Degrees of freedom                                 0
  Minimum Function Value               0.0000000000000

Parameter Estimates:

  Information                                 Expected
  Standard Errors                             Standard

Regressions:
                   Estimate  Std.Err  z-value  P(>|z|)
  Y1 ~                                                
    M1        (b1)    0.353    0.102    3.463    0.001
    M2        (b2)    0.307    0.098    3.126    0.002
    X1        (c1)   -0.040    0.136   -0.295    0.768
  M1 ~                                                
    X1        (a1)    0.565    0.112    5.025    0.000
  M2 ~                                                
    X1        (a2)    0.576    0.126    4.576    0.000
    M1        (d1)    0.286    0.100    2.865    0.004

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)
   .Y1                1.216    0.172    7.071    0.000
   .M1                1.263    0.179    7.071    0.000
   .M2                1.262    0.179    7.071    0.000

Defined Parameters:
                   Estimate  Std.Err  z-value  P(>|z|)
    ind1              0.200    0.070    2.851    0.004
    ind2              0.177    0.068    2.581    0.010
    secondInd1        0.050    0.026    1.947    0.051
    total1            0.386    0.126    3.068    0.002

To visualize model, you can use semPaths() function from semPlot package. But, semPaths() function is inappropriate for models with two or more mediators.

#par(family="NanumGothic")
semPaths(fit,rotation=2,whatLabels = "std")

mediationPlot() can visualize models with multiple mediators, multiple independent variables and/or multiple dependent variables.

mediationPlot(fit)

You can show standardized parameter estimates(default) or parameter estimates(est), or names(name) by settting the parameter “whatLabels”.

mediationPlot(fit,whatLabels = "name")

mediationPlot(fit,whatLabels = "est")

You can visualize indirect effects of this model.

mediationPlot(fit,regression=FALSE,residuals=FALSE,indirect=TRUE)

You can also visualize secondary indirect effect of this model.

mediationPlot(fit,regression=FALSE,residuals=FALSE,secondIndirect=TRUE)

More example

You can extend model with two predictors, two mediators and one response variable.

model=makeEquation(X=c("X1","X2"),M=c("M1","M2"),Y="Y1")
cat(model)
# Mediation Effect
Y1~b1*M1+b2*M2+c1*X1+c2*X2
M1~a1*X1+a2*X2
M2~a3*X1+a4*X2+d1*M1
ind1:=a1*b1
ind2:=a2*b1
ind3:=a3*b2
ind4:=a4*b2
secondInd1:=a1*d1*b2+a2*d1*b2
total1:=c1+c2+a1*b1+a2*b1+a3*b2+a4*b2+a1*d1*b2+a2*d1*b2
fit=sem(model,data=data)
semPaths(fit)

mediationPlot(fit)

Example 3

You can extend model with two predictors, two mediators and two response variables.

model=makeEquation(X=c("X1","X2"),M=c("M1","M2"),Y=c("Y1","Y2"))
cat(model)
# Mediation Effect
Y1~b1*M1+b2*M2+c1*X1+c2*X2
Y2~b3*M1+b4*M2+c3*X1+c4*X2
M1~a1*X1+a2*X2
M2~a3*X1+a4*X2+d1*M1
ind1:=a1*b1
ind2:=a2*b1
ind3:=a3*b2
ind4:=a4*b2
ind5:=a1*b3
ind6:=a2*b3
ind7:=a3*b4
ind8:=a4*b4
secondInd1:=a1*d1*b2+a2*d1*b2
secondInd2:=a1*d1*b4+a2*d1*b4
total1:=c1+c2+a1*b1+a2*b1+a3*b2+a4*b2+a1*d1*b2+a2*d1*b2
total2:=c3+c4+a1*b3+a2*b3+a3*b4+a4*b4+a1*d1*b4+a2*d1*b4
fit=sem(model,data=data)
semPaths(fit)

mediationPlot(fit)

Example 4

model <- '
   # latent variables
     ind60 =~ x1 + x2 + x3
     dem60 =~ y1 + y2 + y3 + y4
     dem65 =~ y5 + y6 + y7 + y8
   # regressions
     dem60 ~ ind60
     dem65 ~ ind60 + dem60
'
fit <- sem(model,
           data=PoliticalDemocracy)
semPaths(fit,rotation=2,whatLabels = "std",groups="manlat",pastel=TRUE)

mediationPlot(fit)

Example 5

Next example is about the teacher’s intervention for students with ADHD. The ADHD data contained in semMediation package is a dataset contains measures about the teacher’s knowlege, empathy and intervention about attention-deficit hyperactivity disorder(ADHD).(Source:Effects of teacher’s knowledge and empathy on educational intervention for ADHD: Focused on the mediating effet of empathy. J Korean Acad Psychiatr Ment Health Nurs 2013:22;45-55.)

Step 1. Define latent variables

First of all, let’s start from defining the latent variables. You can define latent variables with =~ operator.

model='
knowledge =~ general+symptoms+treatmt
empathy =~ cognitiv+emotion+disposit+attitude
intervention =~ classrm+instruct
'

Step 2. Define the mediation effect

You can define mediation effects using the makeEquation() function. Here the knowledge is an independent variable(X) and intervention is a response variable(Y) and empathy is a mediator(M).

mediationModel=makeEquation(X="knowledge",M="empathy",Y="intervention")
cat(mediationModel)
# Mediation Effect
intervention~b1*empathy+c1*knowledge
empathy~a1*knowledge
ind1:=a1*b1
total1:=c1+a1*b1

You can add this mediation effects to the model. The final model is as follows.

model=paste0(model,mediationModel)
cat(model)

knowledge =~ general+symptoms+treatmt
empathy =~ cognitiv+emotion+disposit+attitude
intervention =~ classrm+instruct
# Mediation Effect
intervention~b1*empathy+c1*knowledge
empathy~a1*knowledge
ind1:=a1*b1
total1:=c1+a1*b1

Step 3. Fit the model

You can fit this model by using the lavaan::sem() function.

fit=sem(model,data=ADHD)
summary(fit,standardized= TRUE ,fit.measures= FALSE ,rsquare= FALSE ,modindices= FALSE )
lavaan (0.5-22) converged normally after  68 iterations

  Number of observations                           334

  Estimator                                         ML
  Minimum Function Test Statistic               31.185
  Degrees of freedom                                24
  P-value (Chi-square)                           0.149

Parameter Estimates:

  Information                                 Expected
  Standard Errors                             Standard

Latent Variables:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
  knowledge =~                                                          
    general           1.000                               1.850    0.748
    symptoms          0.582    0.068    8.508    0.000    1.077    0.635
    treatmt           0.746    0.086    8.631    0.000    1.380    0.692
  empathy =~                                                            
    cognitiv          1.000                               2.693    0.812
    emotion           0.950    0.058   16.371    0.000    2.558    0.821
    disposit          1.401    0.091   15.396    0.000    3.772    0.781
    attitude          0.964    0.060   16.155    0.000    2.595    0.812
  intervention =~                                                       
    classrm           1.000                               5.782    0.990
    instruct          0.853    0.046   18.628    0.000    4.934    0.886

Regressions:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
  intervention ~                                                        
    empathy   (b1)    1.241    0.120   10.356    0.000    0.578    0.578
    knowledge (c1)    0.106    0.174    0.608    0.543    0.034    0.034
  empathy ~                                                             
    knowledge (a1)    0.304    0.100    3.033    0.002    0.209    0.209

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
   .general           2.691    0.405    6.649    0.000    2.691    0.440
   .symptoms          1.718    0.179    9.626    0.000    1.718    0.597
   .treatmt           2.076    0.252    8.244    0.000    2.076    0.522
   .cognitiv          3.735    0.386    9.683    0.000    3.735    0.340
   .emotion           3.162    0.334    9.466    0.000    3.162    0.326
   .disposit          9.110    0.882   10.330    0.000    9.110    0.390
   .attitude          3.480    0.359    9.692    0.000    3.480    0.341
   .classrm           0.652    1.504    0.433    0.665    0.652    0.019
   .instruct          6.638    1.209    5.488    0.000    6.638    0.214
    knowledge         3.423    0.549    6.240    0.000    1.000    1.000
   .empathy           6.934    0.812    8.542    0.000    0.956    0.956
   .intervention     21.955    2.353    9.331    0.000    0.657    0.657

Defined Parameters:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
    ind1              0.377    0.128    2.952    0.003    0.121    0.121
    total1            0.483    0.202    2.393    0.017    0.154    0.154

Step 4. Draw a plot

You can draw figure using semPlot::semPaths() function. Alternatively, you can draw a plot with semMediation::mediationPlot().

semPaths(fit,rotation=2,whatLabels = "std")

#mediationPlot(fit,whatLabels = "name",residuals=FALSE)
mediationPlot(fit,whatLabels = "std",residuals = FALSE,base_size=3)

Step 5. Draw a correlation plot

You can draw correlation plot by CorPlot function

require(mycor)
require(ggplot2)

corPlot(fit)

corTable(fit)

general

symptoms

treatmt

cognitiv

emotion

disposit

attitude

classrm

instruct

general

1

symptoms

0.48***

1

treatmt

0.52***

0.43***

1

cognitiv

0.14*

0.12*

0.20***

1

emotion

0.12*

0.09

0.11*

0.70***

1

disposit

0.12*

0.08

0.17**

0.62***

0.62***

1

attitude

0.06

0.07

0.13*

0.64***

0.68***

0.65***

1

classrm

0.11*

0.09

0.11*

0.47***

0.43***

0.52***

0.47***

1

instruct

0.08

0.06

0.13*

0.43***

0.38***

0.46***

0.42***

0.88***

1

p<0.001;p<0.01;p<0.05

Step 6. Make a model fit measures table

You can make a table summarizing model fit measures.

result=modelFitTable(fit)
df2Flextable(result,vanilla=TRUE)

χ2

df

p

CFI

GFI

AGFI

TLI

RMR

SRMR

RMSEA

AIC

BIC

31.19

24

0.15

1

0.98

0.96

0.99

0.43

0.03

0.03

14179.06

14259.09

Step 7. Make a estimates table

result=estimatesTable(fit,ci=TRUE)
df2Flextable(result,vanilla=TRUE)

Variables

Predictors

B

SE

z

p

β

knowledge

general

1(1-1)

0.00

0.75

knowledge

symptoms

0.58(0.45-0.72)

0.07

8.51

< 0.001

0.63

knowledge

treatmt

0.75(0.58-0.92)

0.09

8.63

< 0.001

0.69

empathy

cognitiv

1(1-1)

0.00

0.81

empathy

emotion

0.95(0.84-1.06)

0.06

16.37

< 0.001

0.82

empathy

disposit

1.4(1.22-1.58)

0.09

15.4

< 0.001

0.78

empathy

attitude

0.96(0.85-1.08)

0.06

16.16

< 0.001

0.81

intervention

classrm

1(1-1)

0.00

0.99

intervention

instruct

0.85(0.76-0.94)

0.05

18.63

< 0.001

0.89

intervention

empathy

1.24(1.01-1.48)

0.12

10.36

< 0.001

0.58

intervention

knowledge

0.11(-0.23-0.45)

0.17

0.61

0.543

0.03

empathy

knowledge

0.3(0.11-0.5)

0.10

3.03

0.002

0.21

Step 8. Make a mediation effects table

result=estimatesTable(fit,mediation=TRUE,ci=TRUE)
df2Flextable(result,vanilla=TRUE)

Variables

Predictors

label

B

SE

z

p

β

intervention

empathy

b1

1.24(1.01-1.48)

0.12

10.36

< 0.001

0.58

intervention

knowledge

c1

0.11(-0.23-0.45)

0.17

0.61

0.543

0.03

empathy

knowledge

a1

0.3(0.11-0.5)

0.10

3.03

0.002

0.21

indirect effect

a1*b1

ind1

0.38(0.13-0.63)

0.13

2.95

0.003

0.12

total effect

c1+a1*b1

total1

0.48(0.09-0.88)

0.20

2.39

0.017

0.15

You can draw direct and indirect effects.

mediationPlot(fit,whatLabels = "std",residuals = FALSE,indirect=TRUE,
              regression=TRUE,mediationOnly = TRUE)