# 1 Quick review of moderation and mediation

## 1.2 Mediation

• Mediation tests a hypothetical causal chain where the effect of one variable (X) on another variable (Y) is mediated, or explained, by a third variable (M)
• X -> M -> Y
• For a review see Chapter 14: Mediation and Moderation

# 2 What is moderated mediation?

## 2.1 Conceptual definition

• Moderated mediation tests the influence of a fourth (or more) variable on the mediated relationship between X and Y
• The effect of the mediator is moderated by another variable
• X -> M -> Y (depending on Z)
• The moderation can occur on any and all paths in the mediation model (e.g., a path, b path, c path, or any combination of the three)

## 2.2 Practical definition and example

• The more time one spends in graduate school, the more job offers they have when they graduate
• This relationship is explained by increased publications (i.e., the more time spent in grad school, the more publications one has, and the more publications one has, the more job offers they get)
• However, this causal chain may only work for people who spend their time in graduate school wisely (i.e., spend time with Professor Demos)
• How does spending time with Professor Demos impact the causal chain between time spent in graduate school, publications, and job offers? Let’s find out…

# 3 Moderated mediation data example

## 3.1 Describe the dataset

We are going to simulate a dataset that measured the following:

• X = Time spent in graduate school (we will change the name to “time” when we create the data frame)
• Z = Time spent (hours per week) with Professor Demos in class or in office hours
• M = Number of publications in grad school
• Y = Number of job offers

## 3.2 Create the dataset

We are intentionally creating a moderated mediation effect here and we do so below by setting the relationships (the paths) between our causal chain variables and setting the relationships for our interaction terms

setwd("path of working directory here")

set.seed(42) #This makes sure that everyone gets the same numbers generated through rnorm function

a1 = -.59 #Set the path a1 strength (effect of X on M)

a2 = -.17 #Set path a2 strength (effect of Z on M)

a3 = .29 #Set path a3 strength (interaction between X and Z on M)

b = .59 #Set path b strength (effect of M on Y)

cdash1 = .27 #Set path c'1 strength (effect of X on Y)

cdash2 = .01 #Set path c'2 strength (effect of Z on Y)

cdash3 = -.01 #Set path c'3 strength (interaction betwee X and Z on Y)

Here we are creating the values of our variables for each subject

n <- 200 #Set sample size

X <- rnorm(n, 7, 1) #IV: Time spent in grad school (M = 7, SD = 1)

Z <- rnorm(n, 5, 1) #Moderator: Time spent (hours per week) with Professor Demos in class or in office hours (M = 5, SD = 1)

M <- a1*X + a2*Z + a3*X*Z + rnorm(n, 0, .1) #Mediator: Number of publications in grad school
#The mediator variable is created as a function of the IV, moderator, and their interaction with some random noise thrown in the mix

Y <- cdash1*X + cdash2*Z + cdash3*X*Z + b*M + rnorm(n, 0, .1) #DV: Number of job offers
#Similar to the mediator, the DV is a function of the IV, moderator, their interaction, and the mediator with some random noise thrown in the mix

Now we put it all together and make our data frame

Success.ModMed <- data.frame(jobs = Y, time = X, pubs = M, alex = Z) #Build our data frame and give it recognizable variable names

## 3.3 Examine the dataset and prepare for regression analyses

install.packages("psych") #install this package if not already installed

library(psych) #Helpful for common psych descriptive statistics

str(Success.ModMed) #Examine the structure of the dataset
## 'data.frame':    200 obs. of  4 variables:
##  $jobs: num 3.17 4.54 6.29 7.53 3.32 ... ##$ time: num  8.37 6.44 7.36 7.63 7.4 ...
##  $pubs: num 1.97 5.16 7.79 9.93 2.74 ... ##$ alex: num  3 5.33 6.17 7.06 3.62 ...
round(describe(Success.ModMed)[,c(2:5,8,9,13)], 2) #Put descriptive stats summary into table with only the columns of information that we care about
##        n mean   sd median  min   max   se
## jobs 200 4.60 1.21   4.45 1.83  9.17 0.09
## time 200 6.97 0.97   6.98 4.01  9.70 0.07
## pubs 200 5.14 1.97   4.82 0.01 12.24 0.14
## alex 200 5.01 0.95   4.97 2.30  7.46 0.07

Because we have interaction terms in our regression analyses, we need to mean center our IV and Moderator (Z)

Success.ModMed$time.c <- scale(Success.ModMed$time, center = TRUE, scale = FALSE)[,] #Scale returns a matrix so we have to make it a vector by indexing one column

Success.ModMed$alex.c <- scale(Success.ModMed$alex, center = TRUE, scale = FALSE)[,]

# 4 Moderated mediation analyses using “mediation” package

We will first create two regression models, one looking at the effect of our IVs (time spent in grad school, time spent with Alex, and their interaction) on our mediator (number of publications), and one looking at the effect of our IVs and mediator on our DV (number of job offers).

Next, we will examine the influence of our moderating variable (time spent with Alex) on the mediation effect of time spent in grad school on number of job offers, through number of publications. To do this, we will examine the mediation effect for those who spend a lot of time with Alex versus those who spend little time with Alex.

## 4.1 Create the necessary regression models

We need two regression models to use the mediation package

One model specifies the effect of our IV (time spent in grad school) on our Mediator (number of publications) [and in our case, our moderator (time spent with Alex) and the interaction]

The other model specifies the effect of the IV (time spent in grad school) and Mediator (number of publications) (and possibly moderator as well) on our DV (number of job offers)

install.packages("mediation") #install this first if not already installed

library(mediation)

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

Mod.Med.Model.1<-lm(pubs ~ time.c*alex.c, data = Success.ModMed) #This model predicts number of publications from time spent in grad school, time spent with alex, and the interaction between the two

summary(Mod.Med.Model.1)
##
## Call:
## lm(formula = pubs ~ time.c * alex.c, data = Success.ModMed)
##
## Residuals:
##      Min       1Q   Median       3Q      Max
## -0.24019 -0.06563 -0.00220  0.07089  0.32256
##
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)   5.162911   0.007308  706.49   <2e-16 ***
## time.c        0.857893   0.007518  114.12   <2e-16 ***
## alex.c        1.854224   0.007957  233.02   <2e-16 ***
## time.c:alex.c 0.309273   0.007948   38.91   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.103 on 196 degrees of freedom
## Multiple R-squared:  0.9973, Adjusted R-squared:  0.9973
## F-statistic: 2.418e+04 on 3 and 196 DF,  p-value: < 2.2e-16
• Significant main effect of time spent in grad school on number of publications
• Significant main effect of time spent with Alex on number of publications
• Significant interaction between time spent in grad school and time spent with Alex on number of publications
Mod.Med.Model.2<-lm(jobs ~ time.c*alex.c + pubs, data = Success.ModMed) #This model predicts number of job offers from time spent in grad school, time spent with alex, number of publications, and the interaction between time spent in grad school and time spent with alex

summary(Mod.Med.Model.2)
##
## Call:
## lm(formula = jobs ~ time.c * alex.c + pubs, data = Success.ModMed)
##
## Residuals:
##       Min        1Q    Median        3Q       Max
## -0.294627 -0.061596  0.003638  0.073073  0.226142
##
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)
## (Intercept)    1.601513   0.335390   4.775 3.52e-06 ***
## time.c         0.219265   0.056137   3.906 0.000129 ***
## alex.c        -0.048126   0.120647  -0.399 0.690404
## pubs           0.583970   0.064949   8.991  < 2e-16 ***
## time.c:alex.c -0.009457   0.021347  -0.443 0.658260
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.09367 on 195 degrees of freedom
## Multiple R-squared:  0.9941, Adjusted R-squared:  0.994
## F-statistic:  8191 on 4 and 195 DF,  p-value: < 2.2e-16
• Significant main effect of time spent in grad school on number of job offers
• No effect of time spent with Alex on number of job offers
• Significant main effect of number of publications on number of job offers
• No interaction between time spent in grad school and time spent with Alex on number of job offers

## 4.2 Examine the effect of our moderator on the mediation effect

In this mediation package we list the moderator as a covariate and set the levels to what we want

We can use the +/- 1SD from the mean (or another value that is theoretically important)

This allows us to view impact of the moderator on the direct and indirect effects

Lets look at grad students who spend little time with Alex first

#Moderator must be in both models for mediate to work.
low.alex<-mean(Success.ModMed$alex.c)-sd(Success.ModMed$alex.c) #Sets our level for 1 SD below mean of alex.c

low.alex #Check value of variable
## [1] -0.9470806
Mod.Med.LowAlex <- mediate(Mod.Med.Model.1, Mod.Med.Model.2,
covariates = list(alex.c = low.alex), boot = TRUE,
boot.ci.type = "bca", sims = 10, treat="time.c", mediator="pubs")
#The mediate function can handle different types of CI estimation. Here we are asking for bias-corrected and accelerated confidence intervals because this gives us more accurate confident interval estimates and corrects for deviation from normality
#We also have to specify our IV (treat) and Mediator(pubs)
#For demonstration I am only doing 10 simulations, but in reality you'd want to do at least 2,000

For a review on bootstrapping techniques, see Efron, 2003

• ACME: 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”)
summary(Mod.Med.LowAlex)
##
## Causal Mediation Analysis
##
## Nonparametric Bootstrap Confidence Intervals with the BCa Method
##
## (Inference Conditional on the Covariate Values Specified in covariates')
##
##                Estimate 95% CI Lower 95% CI Upper p-value
## ACME              0.330        0.263         0.39  <2e-16 ***
## ADE               0.228        0.166         0.30  <2e-16 ***
## Total Effect      0.558        0.532         0.57  <2e-16 ***
## Prop. Mediated    0.591        0.467         0.70  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sample Size Used: 200
##
##
## Simulations: 10
plot(Mod.Med.LowAlex, xlim = 0:1)

• Significant direct effect of time spent in grad school on job offers (for those who don’t spend a lot of time with Alex)
• Significant indirect effect of time spent in grad school on job offers through publications (for those who don’t spend a lot of time with Alex)

Now let’s look at grad students who spend a lot of time with Alex

high.alex<-mean(Success.ModMed$alex.c)+sd(Success.ModMed$alex.c)

high.alex
## [1] 0.9470806
Mod.Med.HighAlex <- mediate(Mod.Med.Model.1, Mod.Med.Model.2,
covariates = list(alex.c = high.alex), boot = TRUE,
boot.ci.type = "bca", sims = 10, treat="time.c", mediator="pubs")

summary(Mod.Med.HighAlex)   
##
## Causal Mediation Analysis
##
## Nonparametric Bootstrap Confidence Intervals with the BCa Method
##
## (Inference Conditional on the Covariate Values Specified in covariates')
##
##                Estimate 95% CI Lower 95% CI Upper p-value
## ACME              0.672        0.565         0.75  <2e-16 ***
## ADE               0.210        0.100         0.32  <2e-16 ***
## Total Effect      0.882        0.873         0.89  <2e-16 ***
## Prop. Mediated    0.762        0.641         0.89  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sample Size Used: 200
##
##
## Simulations: 10
plot(Mod.Med.HighAlex, xlim = 0:1)

• Significant direct effect of time spent in grad school on job offers (for those who spend a lot of time with Alex)
• Significant indirect effect of time spent in grad school on job offers through publications (for those who spend a lot of time with Alex)
• The indirect effect looks larger for those who spend a lot of time with Alex compared to those who don’t, but we can test this to make sure

The following code tests whether the difference between indirect effects at each level of the moderator is significantly different from zero

Mod.Med.TestAlex <- mediate(Mod.Med.Model.1, Mod.Med.Model.2, boot = TRUE,
boot.ci.type = "bca", sims = 10, treat="time.c", mediator="pubs")   #We don't specify anything about the moderator in this code yet

test.modmed(Mod.Med.TestAlex, covariates.1 = list(alex.c = low.alex),
covariates.2 = list(alex.c = high.alex), sims = 10) #Here we specify both levels of the moderator that we want to test
##
##  Test of ACME(covariates.1) - ACME(covariates.2) = 0
##
## data:  estimates from Mod.Med.TestAlex
## ACME(covariates.1) - ACME(covariates.2) = -0.3421, p-value <
## 2.2e-16
## alternative hypothesis: true ACME(covariates.1) - ACME(covariates.2) is not equal to 0
## 95 percent confidence interval:
##  -0.5540589 -0.2094565
##
##
##
## data:  estimates from Mod.Med.TestAlex
## alternative hypothesis: true ADE(covariates.1) - ADE(covariates.2) is not equal to 0
## 95 percent confidence interval:
##  -0.1190875  0.2212338
• We can see that the indirect effects are significantly different such that the effect of spending time in graduate school on getting job offers through publications is stronger for those students who spend a lot of time with Alex compared to those who do not
• There is no different in the size of the direct effects, however

## 4.3 Strengths and limitations of “mediation” package

• Code is fairly straightforward and makes intuitive sense in how to specify levels of moderators
• Compatible with many types of regression, including linear, glm, ordered, censored, quantile, GAM, and survival
• Limited in the types of moderated mediation models it can estimate
• Must include moderator in both models (meaning that you cannot model two of the most popular moderated mediation models, Hayes’ Model 7 and Model 14)
• Cannot handle highly complex mediational models with several causally dependent mediators and moderators
• However, structural equation model (SEM) programs can model more complex models, which we turn to next

# 5 Moderated mediation analyses using “lavaan” package

In “lavaan” we specify all regressions and relationships between our variables in one object

We can specify the effects we want to see in our output (e.g., direct, indirect, etc.)

We can also compute means and standard deviations for use in simple slopes analyses

After specifying all the necessary components, we fit the model using an SEM function

install.packages("lavaan") #install this first if not already installed

library(lavaan)

Mod.Med.Lavaan <- '
#Regressions
#These are the same regression equations from our previous example
#Except in this code we are naming the coefficients that are produced from the regression equations
#E.g., the regression coefficient for the effect of time on pubs is named "a1"
pubs ~ a1*time.c + a2*alex.c + a3*time.c:alex.c
jobs ~ cdash1*time.c + cdash2*alex.c + cdash3*time.c:alex.c + b1*pubs

#Mean of centered alex (for use in simple slopes)
#This is making a coefficient labeled "alex.c.mean" which equals the intercept because of the "1"
#(Y~1) gives you the intercept, which is the mean for our alex.c variable
alex.c ~ alex.c.mean*1

#Variance of centered alex (for use in simple slopes)
#This is making a coefficient labeled "alex.c.var" which equals the variance because of the "~~"
#Two tildes separating the same variable gives you the variance
alex.c ~~ alex.c.var*alex.c

#Indirect effects conditional on moderator (a1 + a3*ModValue)*b1
indirect.SDbelow := (a1 + a3*(alex.c.mean-sqrt(alex.c.var)))*b1
indirect.SDabove := (a1 + a3*(alex.c.mean+sqrt(alex.c.var)))*b1

#Direct effects conditional on moderator (cdash1 + cdash3*ModValue)
#We have to do it this way because you cannot call the mean and sd functions in lavaan package
direct.SDbelow := cdash1 + cdash3*(alex.c.mean-sqrt(alex.c.var))
direct.SDabove := cdash1 + cdash3*(alex.c.mean+sqrt(alex.c.var))

#Total effects conditional on moderator
total.SDbelow := direct.SDbelow + indirect.SDbelow
total.SDabove := direct.SDabove + indirect.SDabove

#Proportion mediated conditional on moderator
#To match the output of "mediate" package
prop.mediated.SDbelow := indirect.SDbelow / total.SDbelow
prop.mediated.SDabove := indirect.SDabove / total.SDabove

#Index of moderated mediation
#An alternative way of testing if conditional indirect effects are significantly different from each other
index.mod.med := a3*b1
'

Now we take the specified models and all of the effects we want to estimate and run them through the SEM function. The SEM function allows a completely user-defined model to be fit to the data, like our specifically defined moderated mediation model (the SEM function was designed to fit structural equation models, but can also fit “regular” regression models as well).

#Fit model
Mod.Med.SEM <- sem(model = Mod.Med.Lavaan,
data = Success.ModMed,
se = "bootstrap",
bootstrap = 10)
#Fit measures
summary(Mod.Med.SEM,
fit.measures = FALSE,
standardize = TRUE,
rsquare = TRUE)
• The first chunk of the output show fit indices related to SEM (not really applicable for our purposes)
• The second part of the output shows our regression formulas
• The end of the output shows the specified direct, indirect, total, proportion mediated effects
## lavaan (0.5-23.1097) converged normally after 107 iterations
##
##   Number of observations                           200
##
##   Estimator                                         ML
##   Minimum Function Test Statistic               12.599
##   Degrees of freedom                                 2
##   P-value (Chi-square)                           0.002
##
## Parameter Estimates:
##
##   Information                                 Observed
##   Standard Errors                            Bootstrap
##   Number of requested bootstrap draws               10
##   Number of successful bootstrap draws              10
##
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   pubs ~
##     time.c    (a1)    0.858    0.007  123.187    0.000    0.858    0.425
##     alex.c    (a2)    1.854    0.007  283.790    0.000    1.854    0.892
##     tm.c:l.   (a3)    0.309    0.009   35.101    0.000    0.309    0.149
##   jobs ~
##     time.c  (cds1)    0.219    0.070    3.154    0.002    0.219    0.175
##     alex.c  (cds2)   -0.048    0.151   -0.319    0.750   -0.048   -0.037
##     tm.c:l. (cds3)   -0.009    0.024   -0.389    0.697   -0.009   -0.007
##     pubs      (b1)    0.584    0.081    7.217    0.000    0.584    0.942
##
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   time.c ~~
##     time.c:alex.c    -0.008    0.117   -0.067    0.947   -0.008   -0.009
##
## Intercepts:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##     alex.c  (al..)   -0.000    0.087   -0.000    1.000   -0.000   -0.000
##    .pubs              5.163    0.007  780.718    0.000    5.163    2.629
##    .jobs              1.602    0.419    3.825    0.000    1.602    1.316
##     time.c           -0.000    0.057   -0.000    1.000   -0.000   -0.000
##     tm.c:l.          -0.074    0.062   -1.195    0.232   -0.074   -0.078
##
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##     alex.c  (al..)    0.892    0.041   21.967    0.000    0.892    1.000
##    .pubs              0.010    0.001    9.366    0.000    0.010    0.003
##    .jobs              0.009    0.001    8.741    0.000    0.009    0.006
##     time.c            0.945    0.106    8.919    0.000    0.945    1.000
##     tm.c:l.           0.889    0.141    6.309    0.000    0.889    1.000
##
## R-Square:
##                    Estimate
##     pubs              0.997
##     jobs              0.994
##
## Defined Parameters:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##     indirect.SDblw    0.330    0.060    5.526    0.000    0.330    0.260
##     indirect.SDabv    0.672    0.106    6.355    0.000    0.672    0.540
##     direct.SDbelow    0.228    0.050    4.529    0.000    0.228    0.182
##     direct.SDabove    0.210    0.096    2.183    0.029    0.210    0.168
##     total.SDbelow     0.559    0.015   36.414    0.000    0.559    0.443
##     total.SDabove     0.882    0.015   59.084    0.000    0.882    0.708
##     prp.mdtd.SDblw    0.591    0.095    6.217    0.000    0.591    0.588
##     prop.mdtd.SDbv    0.762    0.113    6.741    0.000    0.762    0.763
##     index.mod.med     0.181    0.028    6.381    0.000    0.181    0.140

We can also call for bootstrapped confidence interval parameter estimates of all of our effects

#Bootstraps
parameterEstimates(Mod.Med.SEM,
boot.ci.type = "bca.simple",
level = .95, ci = TRUE,
standardized = FALSE)[c(19:27),c(4:10)] #We index the matrix to only display columns we are interested in
##                    label   est    se      z pvalue ci.lower ci.upper
## 19      indirect.SDbelow 0.330 0.060  5.526  0.000    0.229    0.439
## 20      indirect.SDabove 0.672 0.106  6.355  0.000    0.471    0.837
## 21        direct.SDbelow 0.228 0.050  4.529  0.000    0.141    0.322
## 22        direct.SDabove 0.210 0.096  2.183  0.029    0.046    0.386
## 23         total.SDbelow 0.559 0.015 36.414  0.000    0.546    0.575
## 24         total.SDabove 0.882 0.015 59.084  0.000    0.852    0.892
## 25 prop.mediated.SDbelow 0.591 0.095  6.217  0.000    0.415    0.757
## 26 prop.mediated.SDabove 0.762 0.113  6.741  0.000    0.549    0.948
## 27         index.mod.med 0.181 0.028  6.381  0.000    0.131    0.221
• Our estimates and confidence intervals are almost identical to the “mediation” package estimates
• The difference is most likely a result of bootstrap estimation differences (e.g., lavaan uses bias-corrected but not accelerated bootstrapping for their confidence intervals)

## 5.1 Strengths and limitations of “lavaan” package

• Extremely customizable
• Can also model latent variables if your measurement model requires it
• Tedious! It took me several hours to figure out how the naming conventions worked
• A lot of up front coding required meaning you kind of need to know exactly what you’re looking for in your model