Overview

Although mediation is extremely difficult to implement in a way that supports the desired inferences, there is some potential for using experience sampling to examine within-person mediation. This approach may provide for more sound inferences than the usual cross-sectional (between-person) mediation. In this tutorial, we illustrate how the within-person (1-1-1) mediation model if fit in a multilevel modeling framework. Generally, we follow the example in Bolger and Laurenceau (2013) Chapter 9: Within-subject Mediation Analysis, but also make use of other resources.

Mediation is said to occur when the effect of one variable (x) on another (y) is transmitted through an intervening variable (m).

There are many resources on mediation, but not so many R-based examples and implementations of within-person mediation. We've done our best to compile these together into a workable example. A good set of general resources is here (note also their caution about one section of code at the beginning) https://stats.idre.ucla.edu/r/faq/how-can-i-perform-mediation-with-multilevel-data-method-2/ and also on the sites of researchers studying mediation, e.g., http://quantpsy.org/, http://afhayes.com/index.html, http://davidakenny.net/cm/mediate.htm, https://psychology.clas.asu.edu/research/labs/research-prevention-laboratory-mackinnon.

Outline

In this tutorial, we will cover...

  1. Introduction to the Within-Person Mediation Model
  2. An Example using a Multivariate Multilevel Model
  3. An Example using a Bayesian Multivariate Multilevel Model
  4. Moderated Within-Person Mediation Model
  5. Conclusion

Loading Libraries

Loading libraries used in this script.

library(ggplot2); theme_set(theme_bw()) #for data viz
library(lme4) #for multilevel models
library(nlme) #for multilevel models
library(psych) #for describing data
library(tidyr); library(dplyr) #for manipulating data 
library(bmlm) #for Bayesian 1-1-1 mediation
library(broom) #for converting results summary to df

1. Introduction to the Within-Person Mediation Model

In basic form, mediation occurs when the effect of an independent variable (x) on a dependent variable (y) is transmitted via a mediator variable (m). This "mediation effect" is often referred to as the indirect effect of X on Y through M. Mediation models allow researchers to test simple hypotheses about "causal processes." Mediation models often involve parsing the total effect (c) of X on Y into a direct effect (c') and an indirect effect (a × b). These coefficients can be derived by fitting a set of simultaneous equations to sample data using linear regression or path analysis. Specifically, in the within-person world,

\[M_{it} = d_{Mi} + a_{i}X_{it} + e_{Mit}\] \[Y_{it} = d_{Yi} + b_{i}M_{it} + c'_{i}X_{it} + e_{Yit}\] where \(d_{Mi}\) and \(d_{Yi}\) are person-specific intercept terms.

2. An Example using a Multivariate Multilevel Model

We follow the example from Bolger & Laurenceau (2013) Chapter 9. In Session R we examined the extent to which daily work stressors were associated with daily relationship dissatisfaction. Here, we examine the extent to which work dissatisfaction can explain that association. For parsimony of presentation, we use data from the female partners only.

The data are daily repeated measures of work stress, work dissatisfaction, and relationship satisfaction for females within a dyad.

Loading Data

#set filepath for data file
filepath <- "https://quantdev.ssri.psu.edu/sites/qdev/files/B&Lmediation.csv"

#read in the .csv file using the url() function
data <- read.csv(file = url(filepath), header = TRUE)

# Examine first few rows of the data set
head(data, 10)
##     id time timec  freldis   fwkdis fwkstr fwkstrc    fwkdisc    freldisc
## 1  101    1   -10 3.034483 5.590119      3    0.03  0.3601189 -1.60551724
## 2  101    2    -9 4.620690 5.535224      3    0.03  0.3052242 -0.01931034
## 3  101    3    -8 2.850575 3.888381      3    0.03 -1.3416194 -1.78942529
## 4  101    4    -7 6.398467 5.352242      4    1.03  0.1222415  1.75846743
## 5  101    5    -6 2.544061 4.483074      1   -1.97 -0.7469259 -2.09593870
## 6  101    6    -5 5.164751 3.339433      2   -0.97 -1.8905672  0.52475096
## 7  101    7    -4 2.704981 4.135407      3    0.03 -1.0945929 -1.93501916
## 8  101    8    -3 5.003831 5.800549      4    1.03  0.5705489  0.36383142
## 9  101    9    -2 4.099617 5.434584      3    0.03  0.2045837 -0.54038314
## 10 101   10    -1 5.471264 4.830741      2   -0.97 -0.3992589  0.83126437
##      fwkstrcb   fwkdiscb  freldiscb   fwkstrcw   fwkdiscw  freldiscw          x
## 1  -0.3033333 -0.6227591 -0.1634446  0.3333333  0.9828781 -1.4420726  0.3333333
## 2  -0.3033333 -0.6227591 -0.1634446  0.3333333  0.9279833  0.1441343  0.3333333
## 3  -0.3033333 -0.6227591 -0.1634446  0.3333333 -0.7188603 -1.6259807  0.3333333
## 4  -0.3033333 -0.6227591 -0.1634446  1.3333333  0.7450007  1.9219121  1.3333333
## 5  -0.3033333 -0.6227591 -0.1634446 -1.6666667 -0.1241668 -1.9324941 -1.6666667
## 6  -0.3033333 -0.6227591 -0.1634446 -0.6666667 -1.2678081  0.6881956 -0.6666667
## 7  -0.3033333 -0.6227591 -0.1634446  0.3333333 -0.4718337 -1.7715745  0.3333333
## 8  -0.3033333 -0.6227591 -0.1634446  1.3333333  1.1933081  0.5272760  1.3333333
## 9  -0.3033333 -0.6227591 -0.1634446  0.3333333  0.8273428 -0.3769385  0.3333333
## 10 -0.3033333 -0.6227591 -0.1634446 -0.6666667  0.2235002  0.9947090 -0.6666667
##             m          y
## 1   0.9828781 -1.4420726
## 2   0.9279833  0.1441343
## 3  -0.7188603 -1.6259807
## 4   0.7450007  1.9219121
## 5  -0.1241668 -1.9324941
## 6  -1.2678081  0.6881956
## 7  -0.4718337 -1.7715745
## 8   1.1933081  0.5272760
## 9   0.8273428 -0.3769385
## 10  0.2235002  0.9947090

The variables of interest in this example are fwkstrcw (the predictor: female work stress centered within person),fwkdiscw (the mediator: female work dissatisfaction centered within person),freldiscw (the outcome: female relationship dissatisfaction centered within person). Note that these variables are state variables. They have already been person-mean centered. We are not currently interested in between-person differences, so they have been separated and set aside. In sum, these data are already person-centered. It may take some steps to get your data to this form.

Copies of the three variables of interest have been made and conveniently labeled as x,m, and y for easy conceptualization and operationalization of the mediation model.

Descriptives of within-person variables

Describe the data.

#variables of interest
vars <- c("fwkstrcw", "fwkdiscw", "freldiscw", "x", "m", "y")
#descriptives
describe(data[, vars])
##           vars    n mean   sd median trimmed  mad   min  max range skew
## fwkstrcw     1 2100    0 1.00   0.00       0 1.20 -2.90 2.95  5.86 0.05
## fwkdiscw     2 2100    0 1.13  -0.01       0 1.15 -3.81 3.81  7.63 0.02
## freldiscw    3 2100    0 0.99   0.00       0 0.98 -3.52 3.22  6.74 0.02
## x            4 2100    0 1.00   0.00       0 1.20 -2.90 2.95  5.86 0.05
## m            5 2100    0 1.13  -0.01       0 1.15 -3.81 3.81  7.63 0.02
## y            6 2100    0 0.99   0.00       0 0.98 -3.52 3.22  6.74 0.02
##           kurtosis   se
## fwkstrcw     -0.25 0.02
## fwkdiscw     -0.11 0.02
## freldiscw     0.07 0.02
## x            -0.25 0.02
## m            -0.11 0.02
## y             0.07 0.02

As expected, the variables all have mean = 0, as they have been person-centered. Standard deviations near one are not purposeful, just coincidence.

Plots of within-person associations

We can generate scatterplots showing the associations between variables for each person for each link in our proposed mediation effect.

x --> y: The c path.
fwkstrcw --> freldiscw

#within-person regressions, x --> y
ggplot(data = data[which(data$id <= 106), ], 
       aes(x = fwkstrcw, y = freldiscw, group = id)) +
  geom_point(color = "black", alpha = .7) + 
  xlab("x = Work Stress (within)") +
  ylab("y = Relationship Dissatisfaction (within)") +
  theme_bw() +
  facet_wrap( ~ id)

x --> m: The a path. fwkstrcw --> fwkdiscw

#within-person regressions, x --> m
ggplot(data = data[which(data$id <= 106),], 
       aes(x = fwkstrcw, y = fwkdiscw, group = id)) +
  geom_point(color = "blue",  alpha = .7) + 
  xlab("x = Work Stress (within)") +
  ylab("m = Work Dissatisfaction (within)") +
  theme_bw()+
  facet_wrap( ~ id)

m --> y: The b path.
fwkdiscw --> freldiscw

#within-person regressions, m --> y
ggplot(data = data[which(data$id <= 106),], 
       aes(x = fwkdiscw, y = freldiscw, group = id)) +
  geom_point(color = "red",  alpha = .7) + 
  xlab("m = Work Dissatisfaction (within)") +
  ylab("y = Relationship Dissatisfaction (within)") +
  theme_bw()+
  facet_wrap( ~ id)

Restructuring the Data Set into Double-entry (stacked) data.

As noted in the introduction, the mediation model can be conceived as a multivariate multilevel model where both y and m are outcome variables.

Practically, fitting the bivariate outcome model is facilitated through some data-restructuring, specifically "double-entry" data.

Restructuring is done in just a few steps.

  1. Subset down to the variables of interest.
data1 <- dplyr::select(data, id, time, timec, fwkstrcw, fwkdiscw, freldiscw, x, m, y)
#look at data set
head(data1, 10)
##     id time timec   fwkstrcw   fwkdiscw  freldiscw          x          m
## 1  101    1   -10  0.3333333  0.9828781 -1.4420726  0.3333333  0.9828781
## 2  101    2    -9  0.3333333  0.9279833  0.1441343  0.3333333  0.9279833
## 3  101    3    -8  0.3333333 -0.7188603 -1.6259807  0.3333333 -0.7188603
## 4  101    4    -7  1.3333333  0.7450007  1.9219121  1.3333333  0.7450007
## 5  101    5    -6 -1.6666667 -0.1241668 -1.9324941 -1.6666667 -0.1241668
## 6  101    6    -5 -0.6666667 -1.2678081  0.6881956 -0.6666667 -1.2678081
## 7  101    7    -4  0.3333333 -0.4718337 -1.7715745  0.3333333 -0.4718337
## 8  101    8    -3  1.3333333  1.1933081  0.5272760  1.3333333  1.1933081
## 9  101    9    -2  0.3333333  0.8273428 -0.3769385  0.3333333  0.8273428
## 10 101   10    -1 -0.6666667  0.2235002  0.9947090 -0.6666667  0.2235002
##             y
## 1  -1.4420726
## 2   0.1441343
## 3  -1.6259807
## 4   1.9219121
## 5  -1.9324941
## 6   0.6881956
## 7  -1.7715745
## 8   0.5272760
## 9  -0.3769385
## 10  0.9947090
  1. Combine the two outcome variables (y and m) into a single z variable, resulting in a "longer" dataset. We will use the pivot_longer() function in the tidyr library.
# Making data longer by putting m and y into a single column, z
datalong <- pivot_longer(data = data1,
                         cols = c("m", "y"), # variables we want to combine
                         names_to = "dv", # column that has the variable names
                         values_to = "z") # column with m and y values

# Reorder rows for convenience 
datalong <- arrange(datalong, id, time, dv)

#look at updated data set
head(datalong, 10)
## # A tibble: 10 x 9
##       id  time timec fwkstrcw fwkdiscw freldiscw      x dv         z
##    <int> <int> <int>    <dbl>    <dbl>     <dbl>  <dbl> <chr>  <dbl>
##  1   101     1   -10    0.333    0.983    -1.44   0.333 m      0.983
##  2   101     1   -10    0.333    0.983    -1.44   0.333 y     -1.44 
##  3   101     2    -9    0.333    0.928     0.144  0.333 m      0.928
##  4   101     2    -9    0.333    0.928     0.144  0.333 y      0.144
##  5   101     3    -8    0.333   -0.719    -1.63   0.333 m     -0.719
##  6   101     3    -8    0.333   -0.719    -1.63   0.333 y     -1.63 
##  7   101     4    -7    1.33     0.745     1.92   1.33  m      0.745
##  8   101     4    -7    1.33     0.745     1.92   1.33  y      1.92 
##  9   101     5    -6   -1.67    -0.124    -1.93  -1.67  m     -0.124
## 10   101     5    -6   -1.67    -0.124    -1.93  -1.67  y     -1.93
  1. Add dummy variables that indicate row entry.
#adding the double indicators 
datalong$dy <- ifelse(datalong$dv == "y", 1, 0)
datalong$dm <- ifelse(datalong$dv == "m", 1, 0)
datalong$dvnum <- ifelse(datalong$dv == "m", 1, 0)

#look at updated data set
head(datalong, 10)
## # A tibble: 10 x 12
##       id  time timec fwkstrcw fwkdiscw freldiscw      x dv         z    dy    dm
##    <int> <int> <int>    <dbl>    <dbl>     <dbl>  <dbl> <chr>  <dbl> <dbl> <dbl>
##  1   101     1   -10    0.333    0.983    -1.44   0.333 m      0.983     0     1
##  2   101     1   -10    0.333    0.983    -1.44   0.333 y     -1.44      1     0
##  3   101     2    -9    0.333    0.928     0.144  0.333 m      0.928     0     1
##  4   101     2    -9    0.333    0.928     0.144  0.333 y      0.144     1     0
##  5   101     3    -8    0.333   -0.719    -1.63   0.333 m     -0.719     0     1
##  6   101     3    -8    0.333   -0.719    -1.63   0.333 y     -1.63      1     0
##  7   101     4    -7    1.33     0.745     1.92   1.33  m      0.745     0     1
##  8   101     4    -7    1.33     0.745     1.92   1.33  y      1.92      1     0
##  9   101     5    -6   -1.67    -0.124    -1.93  -1.67  m     -0.124     0     1
## 10   101     5    -6   -1.67    -0.124    -1.93  -1.67  y     -1.93      1     0
## # … with 1 more variable: dvnum <dbl>

The data set datalong is now ready for analysis.

Running the mediation model -

We are now ready to start running the mediation model as a multivariate multilevel model with two outcomes.

We'll construct a model for the mediator variable as outcome with fwkstrcw (x) and timec as predictors; and a model for the y variable as outcome with fwkstrcw (x) and fwkdiscw (m) and timec as predictors. All the variables have been person centered so that there is no need for intercept terms.

We use the two dummy variables (dm and dy) to turn on and off the parameters for each of the outcomes. The parameters invoked with \(dm\) are for the mediator as outcome model to get the a path, and parameters invoked with \(dy\) are for the y as outcome model to get the b and c' paths.

Setting up and fitting the 1-1-1 mediation model using lme().

#lme mediation model
model_lme <- lme(fixed = z ~ -1 + 
                             dm + dm:fwkstrcw + dm:timec + #m as outcome
                             dy + dy:fwkdiscw + dy:fwkstrcw + dy:timec, #y as outcome
                  random = ~ -1  +  dm:fwkstrcw + dy:fwkdiscw + dy:fwkstrcw | id, 
                  weights = varIdent(form = ~ 1 | dvnum), #separate sigma^{2}_{e} for each outcome
                  data = datalong,
                  na.action = na.exclude,
                  control = lmeControl(opt = "optim", maxIter = 200, msMaxIter = 200, niterEM = 50, msMaxEval = 400))

summary(model_lme)
## Linear mixed-effects model fit by REML
##  Data: datalong 
##        AIC      BIC    logLik
##   12195.67 12290.79 -6082.837
## 
## Random effects:
##  Formula: ~-1 + dm:fwkstrcw + dy:fwkdiscw + dy:fwkstrcw | id
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev     Corr         
## dm:fwkstrcw 0.26103020 dm:fwk dy:fwk
## dy:fwkdiscw 0.21769909 0.542        
## fwkstrcw:dy 0.08798665 0.441  0.932 
## Residual    1.08982601              
## 
## Variance function:
##  Structure: Different standard deviations per stratum
##  Formula: ~1 | dvnum 
##  Parameter estimates:
##         1         0 
## 1.0000000 0.8492324 
## Fixed effects: z ~ -1 + dm + dm:fwkstrcw + dm:timec + dy + dy:fwkdiscw + dy:fwkstrcw +      dy:timec 
##                   Value  Std.Error   DF   t-value p-value
## dm          -0.00000623 0.02378243 4094 -0.000262  0.9998
## dy          -0.00011824 0.02019650 4094 -0.005855  0.9953
## dm:fwkstrcw  0.19030912 0.03559543 4094  5.346448  0.0000
## dm:timec    -0.00583706 0.00397737 4094 -1.467570  0.1423
## dy:fwkdiscw  0.14780833 0.02866743 4094  5.155966  0.0000
## fwkstrcw:dy  0.10548277 0.02267412 4094  4.652121  0.0000
## timec:dy    -0.00247302 0.00338844 4094 -0.729839  0.4655
##  Correlation: 
##             dm     dy     dm:fwk dm:tmc dy:fwk fwkst:
## dy           0.000                                   
## dm:fwkstrcw -0.001  0.000                            
## dm:timec     0.000  0.000 -0.012                     
## dy:fwkdiscw  0.000  0.000  0.301  0.001              
## fwkstrcw:dy  0.000 -0.001  0.127  0.001  0.190       
## timec:dy     0.000  0.000  0.000  0.003  0.015 -0.014
## 
## Standardized Within-Group Residuals:
##          Min           Q1          Med           Q3          Max 
## -3.547996904 -0.675619227 -0.007149205  0.670459292  3.349086939 
## 
## Number of Observations: 4200
## Number of Groups: 100

The results look good and closely match the results in the Bolger & Laurenceau (2013) book! Soon, we will fill more interpretation in here of the lme() model output (and use this model to obtain confidence intervals).

Fixed Effects

#pulling out fixed effects info
FE <- fixef(model_lme)
FE
##            dm            dy   dm:fwkstrcw      dm:timec   dy:fwkdiscw 
## -6.234341e-06 -1.182449e-04  1.903091e-01 -5.837061e-03  1.478083e-01 
##   fwkstrcw:dy      timec:dy 
##  1.054828e-01 -2.473016e-03

Let's interpret ...

Intercepts:

  • 0.0000 = intercept in the m as outcome model (\(d_{Mi}\): We expect this to be zero because we person-centered all the data)
  • 0.0001 = intercept in the y as outcome model (\(d_{Yi}\): We expect this to be zero (or extremely close) because we person-centered all the data)

Effects on M (fwkdiscw):

  • 0.19 = effect of x --> m (a: work stressors predicting work dissatisfaction)

  • -0.006 = effect of time --> m (time trend in work dissatisfaction)

Effects on Y (freldiscw):

  • 0.148 = effect of m --> y (b: work dissatisfaction predicting relationship dissatisfaction)

  • 0.105 = effect of x --> y (c': work stressors predicting relationship dissatisfaction, after adjusting for work dissatisfaction)

  • -0.002 = effect of time --> y (time trend in relationship dissatisfaction)

Let's put the parameters into named objects, as these will be useful later.

#making parameter objects
a <- as.numeric(FE[3])
a
## [1] 0.1903091
b <- as.numeric(FE[5])
b
## [1] 0.1478083
cprime <- as.numeric(FE[6])
cprime
## [1] 0.1054828

Random Effects

#pulling out random effects info
VarCorr(model_lme)
## id = pdLogChol(-1 + dm:fwkstrcw + dy:fwkdiscw + dy:fwkstrcw) 
##             Variance    StdDev     Corr         
## dm:fwkstrcw 0.068136768 0.26103020 dm:fwk dy:fwk
## dy:fwkdiscw 0.047392894 0.21769909 0.542        
## fwkstrcw:dy 0.007741651 0.08798665 0.441  0.932 
## Residual    1.187720737 1.08982601

They are a little easier to interpret and work with when placed in objects ...

The variance of the a paths,

#The *variance* of the *a* paths
sig2_a <- as.numeric(VarCorr(model_lme)["dm:fwkstrcw", "Variance"])
sig2_a
## [1] 0.06813677

\(\sigma^2_{a}\),= 0.0681368

The variance of the b paths,

sig2_b <- as.numeric(VarCorr(model_lme)["dy:fwkdiscw", "Variance"])
sig2_b
## [1] 0.04739289

\(\sigma^2_{b}\),is 0.0473929

The variance of the c' paths,

sig2_cprime <- as.numeric(VarCorr(model_lme)["fwkstrcw:dy", "Variance"])
sig2_cprime
## [1] 0.007741651

\(\sigma^2_{c'}\),is 0.0077417

The residual variance of the mediator variable,

sig2_em <- (1.00*as.numeric(VarCorr(model_lme)["Residual", "StdDev"]))^2
sig2_em
## [1] 1.187721

\(\sigma^2_{em}\), calculated as \((0.8492324 + 1.187721)^{2}\) = 1.1877207.

The residual variance of the outcome variable,

sig2_ey <- (0.8492324*as.numeric(VarCorr(model_lme)["Residual", "StdDev"]))^2
sig2_ey
## [1] 0.856579

\(\sigma^2_{ey}\), is calculated as \((0.8492324 + 1.187721)^{2}\) = 0.856579.

These residual variances differ from the values in the book. I am not sure why.

The covariance between the \(a_{j}\) and \(b_{j}\) paths We need to convert the correlation to a covariance, which we can by using a cortocov function. We can just feed in the correlation coefficient and the two variances, and it will return a covariance.

cortocov <- function (r, var1, var2) {
  cov=r*((var1*var2)^0.5)
  return(cov)
}

covajbj <- cortocov(r = as.numeric(VarCorr(model_lme)["dy:fwkdiscw","Corr"]),
         var1 = as.numeric(VarCorr(model_lme)["dm:fwkstrcw","Variance"]),
         var2 = as.numeric(VarCorr(model_lme)["dy:fwkdiscw","Variance"]))

\(\sigma_{ajbj}\) is calculated as \(\sigma_{aj}*r_{ajbj}*\sigma_{bj}\) = 0.0307997

Yay! We can now use all these to calculate the indirect and total effects.

Calculating mediated effect (e.g., indirect effect)

IndirectEffects = ab + \(\sigma_{ajbj}\)
= (0.1903*0.1478) + 0.0308
= 0.0589

indirecteffect <- a*b + covajbj   
indirecteffect
## [1] 0.05892899

Calculating total effect c

c = c' + ab + \(\sigma_{ajbj}\)
= 0.1054 + (0.1903*0.1478) + 0.0308
= 0.1644

totaleffect <- cprime + a*b + covajbj 
totaleffect
## [1] 0.1644118

Calculating percent of total effect accounted for by mediated effect

PercentMediated = (ab + \(\sigma_{ajbj}\)) / (c' + ab + \(\sigma_{ajbj}\))
= IndirectEffect / c
= 0.0589 / 0.1644
= 0.3584

percentmediated <- 100*(indirecteffect/totaleffect)  
percentmediated
## [1] 35.84232

Calculating percent of mediated effect accounted for by ab covariance (\(\sigma_{ajbj}\))

PercentCovariance = \(\sigma_{ajbj}\) / (ab + \(\sigma_{ajbj}\))
= \(\sigma_{ajbj}\) / ie
= 0.0308 / 0.0589
= 0.5226

percentcovariance <- 100*(covajbj/indirecteffect) 
percentcovariance
## [1] 52.26581

All these calculations match the results in Bolger & Laurenceau (2013) Chapter 9 pretty well.

3. An Example using a Bayesian Multivariate Multilevel Model

As we have seen, performing a multilevel mediation analysis required several data processing steps. An additional limitation is that there is currently not an easy way to run tests on the indirect, such as generating its 95% confidence interval. This issue is particularly complicated with indirect effects in general, as the product of two coefficients (a and b in this case) often results in a skewed distribution.

Presently, Bayesian statistics offer a relatively more straightforward way of estimating the uncertainty of an indirect effect. A comprehensive introduction to Bayesian statistics is beyond the scope of this tutorial, but in brief, Bayesian statistics involve generating a distribution of possible parameter values, which allows the analyst to make probability statements about the evidence in favor of (or against) a research hypothesis, given the dataset. Although a feature of Bayesian statistics is the ability to incorporate one's prior beliefs about parameters values through the specification of a prior distribution, it is also feasible to set what is called a "noninformative" prior. This simply means that our prior belief is that parameter values can range from -\(\infty\) to +\(\infty\). With this kind of noninformative prior, we end up with results that are comparable to those that we would have obtained using maximum likelihood estimation under the classical frequentist approach. We will leverage this feature of Bayesian statistics to help us estimate a 1-1-1 mediation effect and its corresponding uncertainty, in the form of a 95% credibility interval.

3a. Example Using the bmlm package in R

We will do this in two ways. First, we will use the bmlm package for R, which stands for "Bayesian multilevel mediation." This package offers wrapper functions that fit a Bayesian model using Stan (Bayesian software). This makes things easier for us because we don't need to know how to program in Stan to use it. You can read more about this package and Bayesian 1-1-1 mediation in Vuorre & Bolger (2018).

A detailed tutorial showing how bmlm can be used for the Chapter 9 data is available here: https://mvuorre.github.io/bmlm/articles/bmlm-blch9/bmlm-blch9.html. We will reproduce and review some of the core steps from this tutorial in this document (full credit to bmlm developer Matti Vuorre!).

First, we can use the mlm function to simply specify the dataset, id variable, x variable, mediator, and outcome. It will take these inputs and automatically fit a stacked multilevel model and compute the indirect effect using more or less the same steps we went through above. Note that we do not need to create the stacked dataset. We only need to have centered our variables.

Note that a limitation of bmlm is that it does not (currently) accommodate the inclusion of covariates, so our results will be a little different from those generated earlier since we are no longer controlling for time. In our earlier model, though, we found that time had essentially no relationship to either our mediator or y variables.

library(bmlm)
# This can take a few minutes to run
fit <- mlm(d = data, 
           id = "id",
           x = "fwkstrcw",
           m = "fwkdiscw",
           y = "freldiscw",
           iter = 4000, 
           cores = 4)
mlm_summary(fit, digits = 3)
##   Parameter  Mean    SE Median  2.5% 97.5% n_eff Rhat
## 1         a 0.189 0.036  0.190 0.118 0.260 10544    1
## 2         b 0.150 0.030  0.151 0.091 0.209  9246    1
## 3        cp 0.104 0.023  0.104 0.058 0.150 11303    1
## 4        me 0.057 0.014  0.056 0.032 0.086  5230    1
## 5         c 0.161 0.026  0.160 0.110 0.214 10801    1
## 6       pme 0.356 0.081  0.350 0.211 0.525  5490    1

Our indirect effect estimate and its corresponding 95% interval is given by the term me (mediated effect). We can see that it is extremely close to the indirect effect we computed in our earlier step (.06), but an advantage is that this time we can gauge the uncertainty around the effect, 95% CI [.03, .09]. We can say that the indirect effect for the typical person is .06, but it could be as small as .03 or as large as .09.

We can also use helpful plotting functions in bmlm to visualize our results.

We can generate a classic mediation path diagram.

mlm_path_plot(fit, level = .95, text = T,
              xlab = "Work\nstressors",
              mlab = "Work\ndissatisfaction",
              ylab = "Relationship\ndissatisfaction", digits = 2)

We can also visualize the between-subject heterogeneity in the indirect effect. This is one of the most interesting features of 1-1-1 mediation: Because every subject has their own slope for the a-path and b-path, they each have their own indirect effect as well. Here, the subjects are displayed in order of the size of their indirect effect, and the indirect effect in on the y-axis. We can easily see that although some participants show something similar to the average mediated effect (shown in \(\color{red}{\text{red}}\)), others show no mediated effect and others show a mediated effect several times larger than what we found for the prototypical person.

# pull out exact fixed effect for indirect effect
fit_tidy <- tidy(fit)
me_fixed <- subset(fit_tidy, term == "me")$estimate
mlm_pars_plot(fit, pars = c("me", "u_me"), type = "coef", level = .95) + 
  geom_hline(yintercept = me_fixed, lty = "dashed", color = "red", alpha = .5)

3b. Example Using the brms package in R

The second option for fitting a Bayesian version of the 1-1-1 (and generating an uncertainty estimate around the indirect effect) is to use the brms package for R. brms is an acronym that stands for Bayesian Regression Modeling using Stan. It is an extremely powerful and flexible program that can be used for regression models, multilevel models, and generalized linear models.

Data setup

brms can accommodate "true" multivariate models, so we don't need to stack our dataset to fit a 1-1-1 mediation model.

Instead, we will define two "formulas" to fit simultaneously in a single model: one formula in which M is the outcome, and a second in which Y is the outcome. As you will see from the code below, the code is very similar to lme4 syntax. There is one addition we need to make, however, which is specific to multivariate models in brms. We'll add in |p| as part of the random effects statements. The selection of p is arbitrary. The important thing here is to choose the same letter to use in both formulas, as this tells the model that we want to model correlations among random effects from each of the two sub-models that will make up our overall mediation model. This is essential for 1-1-1 mediation so we can get an estimate of the correlation between the a and b paths, which we will later convert into a covariance.

library (brms)

xm <- bf(fwkdiscw ~ -1 + timec + fwkstrcw + (-1 + fwkstrcw |p| id))

my <- bf(freldiscw ~  -1 + timec + fwkstrcw + fwkdiscw + (-1 + fwkstrcw + fwkdiscw |p| id))

Fit model

We can use these formulas and plug them into the brm function, which is a general modeling function from brms that uses syntax similar to the lme4 package. We add in set_rescor(FALSE) to let the model know that we don't want to model a residual correlation between M and Y.

## Note: This will probably take awhile to run.
# You could un-comment the "iter = " statement and fill in a smaller number of iterations to test out the model. You will likely need to increase the number of iterations (or use defaults) to achieve model convergence.


fit_brm <- brm(xm + my + set_rescor(FALSE),
               data = data,
               # iter = [fill in number here to test],
               seed = 111
               )

As you can see, the output looks a bit different compared to what we saw with nlme. Instead of having dm:[name of predictor], our predictors automatically come labeled with the suffix indicating the DV to which they correspond. For example, fwkdiscw_timec is the effect of time on female partner work dissatisfaction (fwkdiscw), equivalent to dm:timec in the nlme version of our model.

We also get separate standard deviation of error terms for each outcome automatically under the Family Specific Parameters section. sigma_fwkdiscw is the residual for M (fwkdiscw), and freldiscw is the residual for Y (freldiscw).

print(fit_brm, digits = 4)
##  Family: MV(gaussian, gaussian) 
##   Links: mu = identity; sigma = identity
##          mu = identity; sigma = identity 
## Formula: fwkdiscw ~ -1 + timec + fwkstrcw + (-1 + fwkstrcw | p | id) 
##          freldiscw ~ -1 + timec + fwkstrcw + fwkdiscw + (-1 + fwkstrcw + fwkdiscw | p | id) 
##    Data: data (Number of observations: 2100) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Group-Level Effects: 
## ~id (Number of levels: 100) 
##                                            Estimate Est.Error l-95% CI u-95% CI
## sd(fwkdiscw_fwkstrcw)                        0.2638    0.0356   0.1967   0.3375
## sd(freldiscw_fwkstrcw)                       0.0845    0.0317   0.0170   0.1472
## sd(freldiscw_fwkdiscw)                       0.2226    0.0269   0.1728   0.2777
## cor(fwkdiscw_fwkstrcw,freldiscw_fwkstrcw)    0.3097    0.3097  -0.3433   0.8520
## cor(fwkdiscw_fwkstrcw,freldiscw_fwkdiscw)    0.5132    0.1558   0.1880   0.7868
## cor(freldiscw_fwkstrcw,freldiscw_fwkdiscw)   0.6971    0.2331   0.1407   0.9751
##                                              Rhat Bulk_ESS Tail_ESS
## sd(fwkdiscw_fwkstrcw)                      1.0005     2350     3037
## sd(freldiscw_fwkstrcw)                     1.0013     1356      859
## sd(freldiscw_fwkdiscw)                     1.0006     2200     2856
## cor(fwkdiscw_fwkstrcw,freldiscw_fwkstrcw)  1.0003     2424     2176
## cor(fwkdiscw_fwkstrcw,freldiscw_fwkdiscw)  1.0026     1352     2453
## cor(freldiscw_fwkstrcw,freldiscw_fwkdiscw) 1.0030      548      750
## 
## Population-Level Effects: 
##                    Estimate Est.Error l-95% CI u-95% CI   Rhat Bulk_ESS
## fwkdiscw_timec      -0.0059    0.0040  -0.0135   0.0020 1.0004     6185
## fwkdiscw_fwkstrcw    0.1902    0.0366   0.1184   0.2630 1.0001     3707
## freldiscw_timec     -0.0024    0.0034  -0.0092   0.0043 1.0020     5509
## freldiscw_fwkstrcw   0.1042    0.0233   0.0583   0.1506 1.0007     6734
## freldiscw_fwkdiscw   0.1492    0.0296   0.0911   0.2068 1.0005     4625
##                    Tail_ESS
## fwkdiscw_timec         2561
## fwkdiscw_fwkstrcw      3463
## freldiscw_timec        2920
## freldiscw_fwkstrcw     2872
## freldiscw_fwkdiscw     3246
## 
## Family Specific Parameters: 
##                 Estimate Est.Error l-95% CI u-95% CI   Rhat Bulk_ESS Tail_ESS
## sigma_fwkdiscw    1.0909    0.0176   1.0574   1.1266 1.0030     7452     2792
## sigma_freldiscw   0.9264    0.0147   0.8984   0.9551 1.0021     6952     2893
## 
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, and Rhat is the potential 
## scale reduction factor on split chains (at convergence, Rhat = 1).

Compute indirect effect

To compute the indirect effect, we first need to convert the correlation of a and b paths to a covariance.

We can do that with this function:

cortocov <- function (r, var1, var2) {
  cov=r*((var1*var2)^0.5)
  return(cov)
}

In order to apply this function, we need to pull out the relevant model parameters: the random effects for a and b (as standard deviations) and their correlation.

We can accomplish this by accessing the posterior distributions from our model. As touched on above, Bayesian estimation involves an iterative process that generates a bunch of possible effect sizes for that parameter. The estimates returned from the model are a summary statistic of all of these effect sizes (usually a mean or median).

## use the posterior_samples() function to pull out the posterior distributions for all model parameters
## these are all possible effect sizes for each parameter
med_post <- posterior_samples(fit_brm)

Now, we have a new dataframe that contains all the possible effect sizes for each model parameter. We can extract the columns for the parameters we need to get from an a-b correlation to a covariance.

We see that after performing this conversion, we get many possible values for cov(a, b). But, if we take the mean (a summary of these possible values), we get .03, which matches what we found using nlme above.

# plug in SD and correlation corresponding to a and b paths for cortocov function

med_post$covab <- cortocov(
  # vector of posterior samples corresponding to correlation of a and b paths
  med_post$cor_id__fwkdiscw_fwkstrcw__freldiscw_fwkdiscw,
  
  # vector of posterior samples corresponding to SD of a path --> convert to variance
  med_post$sd_id__fwkdiscw_fwkstrcw^2,
  
  # vector of posterior samples corresponding to SD of b path --> convert to variance
  med_post$sd_id__freldiscw_fwkdiscw^2)

round(mean(med_post$covab), digits = 4) # rounds to .03
## [1] 0.0302

Now that we have cov(a,b), we can add that in to get our indirect effect: a*b + cov(a,b). We can again use the posterior distributions for a and b and multiply them together. Note that for cov(a,b), we should use the entire vector of possible values rather than the summary statistic only. This will allow us to generate a 95% credibility interval for the indirect effect!

indirect_effect <- 
  med_post$b_fwkdiscw_fwkstrcw *   # a path
   med_post$b_freldiscw_fwkdiscw +   # b path
  med_post$covab   # cov(a, b) as calculated above

We've defined the indirect effect, and this gives us a whole bunch of values--aka, a posterior distribution for the indirect effect, or a vector with all possible effect sizes for the indirect effect.

We can summarize this vector to get a single estimate for the indirect effect, and, importantly, put a credibility interval around this estimate to gauge the (un)certainty. We can use the quantile function to do this. If we ask for .025 and .975 probabilities, this will give us the 95% interval (with 2.5% on each end of the distribution). Asking for the .50 probability will give us the median estimate, although we could also compute the mean (in this particular example, the estimates for the mean and median are essentially the same).

We can see that the 95% credibility interval excludes 0, so we can be quite confident that there is a non-zero indirect effect. In the frequentist framework, we would consider this indirect effect to be statistically significant.

round(quantile(indirect_effect, probs = c(.025, .5, .975)), digits = 3) 
##  2.5%   50% 97.5% 
## 0.034 0.058 0.089
round(mean(indirect_effect), digits = 3)
## [1] 0.059

4. Moderated Within-Person Mediation Model

Using a 1-1-1 mediation analysis, we found that participants varied in the size of their indirect effects. Where does this heterogeneity in the indirect effect come from? We will next investigate whether a between-person variable might moderate the pattern of within-person mediation we found in our example and help explain the heterogeneity we observed.

Imagine that roughly half the participants in this study completed a mindfulness intervention (treat = +.5) vs. a control intervention (treat = -.5). Does having completed a mindfulness intervention weaken the link between work stressors and relationship dissatisfaction, by way of work dissatisfaction?

Simulate Data

First, we need to simulate a variable that contains values for our intervention condition. In 'real' datasets, you would, of course, already have this variable. So this simulation is only for illustrative purposes to help us walk through this fake example.

# simulate a new variable, treat. +.5 = treatment group, -.5 = control group
re <- data.frame(ranef(model_lme))
re$id <- row.names(re)
re <- re %>% arrange(dm.fwkstrcw)
set.seed(111)
re$treat <- c(sample(c(.5, -.5), nrow(re)/2, replace = T, prob = c(.6, .4)), sample(c(.5, -.5), nrow(re)/2, replace = T, prob = c(.4, .6)))

# add this new variable to the original data, resulting in a new dataframe, datalong2
datalong2 <- merge(datalong, dplyr::select(re, id, treat), by = "id")

Run Moderated Mediation Model

#lme moderated mediation model
model_lme_mod <- lme(fixed = z ~ -1 + 
                       # m as outcome
                       dm + dm:fwkstrcw + dm:timec + 
                       # add in treat as main effect and moderator
                       dm:treat + dm:fwkstrcw:treat + 
                       
                       # y as outcome
                       dy + dy:fwkdiscw + dy:fwkstrcw + dy:timec +
                       #  add in treat as main effect and moderator
                       dy:treat + dy:fwkdiscw:treat + dy:fwkstrcw:treat,
                  random = ~ -1  +  dm:fwkstrcw + dy:fwkdiscw + dy:fwkstrcw | id, 
                  weights = varIdent(form = ~ 1 | dvnum), #separate sigma^{2}_{e} for each outcome
                  data = datalong2,
                  na.action = na.exclude,
                  control = lmeControl(opt = "optim", maxIter = 200, msMaxIter = 200, niterEM = 50, msMaxEval = 400))

summary(model_lme_mod)
## Linear mixed-effects model fit by REML
##  Data: datalong2 
##        AIC      BIC    logLik
##   12216.11 12342.91 -6088.057
## 
## Random effects:
##  Formula: ~-1 + dm:fwkstrcw + dy:fwkdiscw + dy:fwkstrcw | id
##  Structure: General positive-definite, Log-Cholesky parametrization
##             StdDev     Corr         
## dm:fwkstrcw 0.25245235 dm:fwk dy:fwk
## dy:fwkdiscw 0.20332688 0.494        
## fwkstrcw:dy 0.09569396 0.438  0.936 
## Residual    1.09019951              
## 
## Variance function:
##  Structure: Different standard deviations per stratum
##  Formula: ~1 | dvnum 
##  Parameter estimates:
##         1         0 
## 1.0000000 0.8489391 
## Fixed effects: z ~ -1 + dm + dm:fwkstrcw + dm:timec + dm:treat + dm:fwkstrcw:treat +      dy + dy:fwkdiscw + dy:fwkstrcw + dy:timec + dy:treat + dy:fwkdiscw:treat +      dy:fwkstrcw:treat 
##                         Value  Std.Error   DF   t-value p-value
## dm                 0.00000644 0.02383347 4089  0.000270  0.9998
## dy                -0.00012377 0.02023291 4089 -0.006117  0.9951
## dm:fwkstrcw        0.19484965 0.03507723 4089  5.554875  0.0000
## dm:timec          -0.00598793 0.00398054 4089 -1.504300  0.1326
## dm:treat           0.00001287 0.04766693 4089  0.000270  0.9998
## dy:fwkdiscw        0.15292577 0.02766610 4089  5.527550  0.0000
## fwkstrcw:dy        0.10543314 0.02312540 4089  4.559193  0.0000
## timec:dy          -0.00223900 0.00339372 4089 -0.659748  0.5095
## treat:dy          -0.00024755 0.04046582 4089 -0.006117  0.9951
## dm:fwkstrcw:treat -0.13915606 0.07018825 4089 -1.982612  0.0475
## treat:dy:fwkdiscw -0.16672114 0.05533202 4089 -3.013105  0.0026
## fwkstrcw:treat:dy -0.01146192 0.04630526 4089 -0.247529  0.8045
##  Correlation: 
##                   dm     dy     dm:fwk dm:tmc dm:trt dy:fwk fwkst: tmc:dy
## dy                 0.000                                                 
## dm:fwkstrcw       -0.001  0.000                                          
## dm:timec           0.000  0.000 -0.015                                   
## dm:treat          -0.060  0.000 -0.001  0.000                            
## dy:fwkdiscw        0.000  0.000  0.262  0.001  0.000                     
## fwkstrcw:dy        0.000 -0.001  0.133  0.000  0.000  0.192              
## timec:dy           0.000  0.000  0.000  0.003  0.000  0.020 -0.021       
## treat:dy           0.000 -0.060  0.000  0.000  0.000  0.000 -0.001  0.000
## dm:fwkstrcw:treat -0.001  0.000 -0.077  0.034 -0.001 -0.016 -0.010 -0.001
## treat:dy:fwkdiscw  0.000  0.000 -0.016  0.000  0.000 -0.063  0.032 -0.020
## fwkstrcw:treat:dy  0.000 -0.001 -0.010  0.001  0.000  0.033 -0.104  0.053
##                   trt:dy dm:fw: trt:d:
## dy                                    
## dm:fwkstrcw                           
## dm:timec                              
## dm:treat                              
## dy:fwkdiscw                           
## fwkstrcw:dy                           
## timec:dy                              
## treat:dy                              
## dm:fwkstrcw:treat  0.000              
## treat:dy:fwkdiscw  0.000  0.261       
## fwkstrcw:treat:dy -0.001  0.133  0.191
## 
## Standardized Within-Group Residuals:
##          Min           Q1          Med           Q3          Max 
## -3.544284304 -0.673535395 -0.009187647  0.668227972  3.327378684 
## 
## Number of Observations: 4200
## Number of Groups: 100

Fixed Effects

Because treat is effect coded, the values we obtain for main effects are estimates averaging over treatment condition.

#pulling out fixed effects info
FE_mod <- lme4::fixef(model_lme_mod)
FE_mod
##                dm                dy       dm:fwkstrcw          dm:timec 
##      6.437258e-06     -1.237743e-04      1.948496e-01     -5.987934e-03 
##          dm:treat       dy:fwkdiscw       fwkstrcw:dy          timec:dy 
##      1.287452e-05      1.529258e-01      1.054331e-01     -2.239001e-03 
##          treat:dy dm:fwkstrcw:treat treat:dy:fwkdiscw fwkstrcw:treat:dy 
##     -2.475485e-04     -1.391561e-01     -1.667211e-01     -1.146192e-02

Let's interpret ...

Intercepts:

  • 0 = intercept in the m as outcome model (\(d_{Mi}\): We expect this to be zero because we person-centered all the data)

  • 0 = intercept in the y as outcome model (\(d_{Yi}\): We expect this to be zero because we person-centered all the data)

Effects on M (fwkdiscw):

  • 0.195 = effect of x --> m (a: work stressors predicting work dissatisfaction across averaging conditions)

  • -0.006 = effect of time --> m (time trend in work dissatisfaction averaging across conditions)

  • 0 = effect of treatment --> m

  • -0.139 = interaction of x and treatment --> m (moderation of a: differences in effect of work stressors predicting work dissatisfaction as a function of treatment )

Effects on Y (freldiscw):

  • 0.153 = effect of m --> y (b: work dissatisfaction predicting relationship dissatisfaction averaging across conditions)

  • 0.105 = effect of x --> y (c': work stressors predicting relationship dissatisfaction, after adjusting for work dissatisfaction, averaging across conditions)

  • -0.002 = effect of time --> y (time trend in relationship dissatisfaction, averaging across conditions)

  • 0 = effect of treatment --> y

  • -0.011 = interaction of x and treatment --> y (moderation of c: differences in effect of work stressors predicting relationship dissatisfaction as a function of treatment)

  • -0.167 = interaction of m and treatment --> y (moderation of b: differences in effect of work dissatisfaction predicting relationship dissatisfaction as a function of treatment)

#making parameter objects

# mediation paths
a_mod <- as.numeric(lme4::fixef(model_lme_mod)["dm:fwkstrcw"])
a_mod
## [1] 0.1948496
b_mod <- as.numeric(lme4::fixef(model_lme_mod)["dy:fwkdiscw"])
b_mod
## [1] 0.1529258
cprime_mod <- as.numeric(lme4::fixef(model_lme_mod)["fwkstrcw:dy"])
cprime_mod
## [1] 0.1054331
# moderation of mediation paths
aXtreat_mod <- as.numeric(lme4::fixef(model_lme_mod)["dm:fwkstrcw:treat"])
aXtreat_mod
## [1] -0.1391561
bXtreat_mod <- as.numeric(lme4::fixef(model_lme_mod)["treat:dy:fwkdiscw"])
bXtreat_mod
## [1] -0.1667211
cprimeXtreat_mod <- as.numeric(lme4::fixef(model_lme_mod)["fwkstrcw:treat:dy"])
cprimeXtreat_mod
## [1] -0.01146192

Random Effects

Random effects can be obtained using similar code above. Of particular interest, we will need the covariance of a and b that is left over after taking our moderator, treat, into account (i.e., their residual covariance).

As before, we need to convert the correlation of these random effects to a covariance.

covajbj_mod <- cortocov(r = as.numeric(VarCorr(model_lme_mod)["dy:fwkdiscw","Corr"]),
         var1 = as.numeric(VarCorr(model_lme_mod)["dm:fwkstrcw","Variance"]),
         var2 = as.numeric(VarCorr(model_lme_mod)["dy:fwkdiscw","Variance"]))
covajbj_mod
## [1] 0.02535719

Calculating mediated effect (e.g., indirect effect)

The indirect effect is again a combination of a*b + cov(ab), but this time we have moderation on both the a and b paths to account for.

Treatment Group (treat = +.5)

First, we compute the indirect effect for the average person in the treatment group.

treat_weight <- .5 #reflects effect coding of treat = +.5

indirecteffect_treat_mod <- (a_mod + treat_weight*aXtreat_mod)*
                            (b_mod + treat_weight*bXtreat_mod) +
                              covajbj_mod   
indirecteffect_treat_mod
## [1] 0.03407174

Control Group (treat = -.5)

Next, we compute the indirect effect for the average person in the control group.

control_weight <- -.5 #reflects effect coding of treat = -.5

indirecteffect_control_mod <- (a_mod + control_weight*aXtreat_mod)*
                              (b_mod + control_weight*bXtreat_mod) +
                              covajbj_mod   
indirecteffect_control_mod
## [1] 0.08783784

Person-specific Mediation Estimates

# pull out person-specific coefficients
ind_ranefs <- data.frame(coef(model_lme_mod))
ind_ranefs$id <- row.names(ind_ranefs)

# combine with info about treatment condition
ind_ranefs_treat <- merge(ind_ranefs, dplyr::select(re, id, treat), by = "id")


# compute person specific indirect effect
ind_ranefs_treat$indirect <-
  # moderated a path
  (ind_ranefs_treat$dm.fwkstrcw + ind_ranefs_treat$treat * ind_ranefs_treat$dm.fwkstrcw.treat) *
  # moderated b path
  (ind_ranefs_treat$dy.fwkdiscw + ind_ranefs_treat$treat * ind_ranefs_treat$treat.dy.fwkdiscw)

Visualizing Moderated 1-1-1 Mediation Effect

# create an empty y-variable strictly for plotting purposes
ind_ranefs_treat$y <- NA


ggplot(ind_ranefs_treat, aes(x = indirect, y = y, color = as.factor(treat))) +
  # add in data points showing subject-specific mediation effect, colored by treat
  geom_point(position = position_jitter(h = .03), alpha = .4, size = 2) + 
  
  # add in large points to indicate the mediation effect for 
  # the average person in each group
  geom_point(aes(x = indirecteffect_control_mod, y = y), color = "red", size = 6) + 
  geom_point(aes(x = indirecteffect_treat_mod, y = y), color = "navyblue", size = 6) +
  
  scale_color_manual(values = c("red", "blue"), labels = c("Control", "Treatment")) +
  theme_bw()+
  xlab("Indirect Effect") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank()) +
  labs(color = "Intervention Condition") +
  theme(legend.position="bottom")

As this figure shows, the 1-1-1 mediation effect differs between conditions. The average person in the treatment condition showed an attenuated indirect effect of daily work stress on relationship dissatisfaction via work dissatisfaction, compared to those in the control condition. Note that we still see some unexplained between-person heterogeneity, as indicated by the smaller dots in the figure above.

To obtain 95% intervals around the indirect effect at varying levels of the moderator, a version of the code provided above using the brms package could be adapted.

For further information and discussion of moderated mediation in multilevel models, you may wish to check out this paper by Bauer, Preacher, and Gil: http://quantpsy.org/pubs/bauer_preacher_gil_2006.pdf. There are some interesting extensions that can be used to obtain bootstrapped confidence intervals for the mediated effects.

5. Conclusion

We have briefly reviewed within-person (1-1-1) mediation models, as well as how to include a between-person moderator of the X-M-Y links.

As always, please be careful with the mediational inferences.

Thanks for playing!

sessionInfo()
## R version 3.6.1 (2019-07-05)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Mojave 10.14.6
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] brms_2.10.0   broom_0.5.2   bmlm_1.3.11   Rcpp_1.0.4.6  dplyr_0.8.3  
##  [6] tidyr_1.0.0   psych_1.8.12  nlme_3.1-140  lme4_1.1-16   Matrix_1.2-17
## [11] ggplot2_3.3.0
## 
## loaded via a namespace (and not attached):
##   [1] minqa_1.2.4           colorspace_1.4-1      rjson_0.2.20         
##   [4] ggridges_0.5.1        rsconnect_0.8.15      htmlTable_1.13.1     
##   [7] markdown_1.1          corpcor_1.6.9         base64enc_0.1-3      
##  [10] rstudioapi_0.10       lavaan_0.6-5          rstan_2.19.2         
##  [13] DT_0.8                fansi_0.4.0           bridgesampling_0.7-2 
##  [16] codetools_0.2-16      splines_3.6.1         mnormt_1.5-5         
##  [19] knitr_1.28            shinythemes_1.1.2     glasso_1.10          
##  [22] zeallot_0.1.0         bayesplot_1.7.0       Formula_1.2-3        
##  [25] nloptr_1.2.1          cluster_2.1.0         png_0.1-7            
##  [28] shiny_1.3.2           compiler_3.6.1        backports_1.1.5      
##  [31] assertthat_0.2.1      cli_1.1.0             later_0.8.0          
##  [34] acepack_1.4.1         htmltools_0.4.0       prettyunits_1.0.2    
##  [37] tools_3.6.1           igraph_1.2.4.1        coda_0.19-3          
##  [40] gtable_0.3.0          glue_1.4.1            reshape2_1.4.3       
##  [43] vctrs_0.2.0           crosstalk_1.0.0       xfun_0.14            
##  [46] stringr_1.4.0         ps_1.3.0              miniUI_0.1.1.1       
##  [49] mime_0.9              ggm_2.3               lifecycle_0.1.0      
##  [52] gtools_3.8.1          MASS_7.3-51.4         zoo_1.8-6            
##  [55] scales_1.0.0          BDgraph_2.61          colourpicker_1.0     
##  [58] promises_1.0.1        Brobdingnag_1.2-6     parallel_3.6.1       
##  [61] inline_0.3.15         huge_1.3.2            shinystan_2.5.0      
##  [64] RColorBrewer_1.1-2    yaml_2.2.1            pbapply_1.4-2        
##  [67] gridExtra_2.3         loo_2.1.0             StanHeaders_2.18.1-10
##  [70] rpart_4.1-15          latticeExtra_0.6-28   stringi_1.4.6        
##  [73] dygraphs_1.1.1.6      checkmate_1.9.4       pkgbuild_1.0.5       
##  [76] rlang_0.4.6           pkgconfig_2.0.3       d3Network_0.5.2.1    
##  [79] matrixStats_0.54.0    evaluate_0.14         lattice_0.20-38      
##  [82] purrr_0.3.3           rstantools_1.5.1      htmlwidgets_1.3      
##  [85] labeling_0.3          processx_3.4.1        tidyselect_0.2.5     
##  [88] plyr_1.8.4            magrittr_1.5          R6_2.4.0             
##  [91] generics_0.0.2        Hmisc_4.2-0           pillar_1.4.2         
##  [94] whisker_0.4           foreign_0.8-71        withr_2.1.2          
##  [97] xts_0.11-2            survival_2.44-1.1     abind_1.4-5          
## [100] nnet_7.3-12           tibble_2.1.3          crayon_1.3.4         
## [103] fdrtool_1.2.15        utf8_1.1.4            rmarkdown_2.3        
## [106] jpeg_0.1-8            grid_3.6.1            qgraph_1.6.3         
## [109] data.table_1.14.0     pbivnorm_0.6.0        callr_3.3.1          
## [112] threejs_0.3.1         digest_0.6.25         xtable_1.8-4         
## [115] httpuv_1.5.1          stats4_3.6.1          munsell_0.5.0        
## [118] shinyjs_1.0