This page edited by Carl F. Falk, Jan 2025

Overview

Although mediation is extremely difficult to implement with the desired inferences, there is some potential for using experience sampling to examine within-person mediation which may provide for more sound inferences than the usual cross-sectional (between-person) mediation. In this tutorial we work through a few examples and ways of implementing the within-person model. 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 workable examples. 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://davidakenny.net/cm/mediate.htm, http://afhayes.com/index.html.

Specifically for multilevel mediation analysis models, there are several foundational articles that may be of interest:

Some code also relies on that introduced by the following papers:

Outline

  1. Introduction to The Within-Person Mediation Model
  2. An Example using Multivariate Multilevel Model
    1. Data, Descriptives, Plots
    2. Estimating the Model Several Ways
    3. A Closer Look at Results
  3. Obtaining Confidence (or Credibility) Intervals for the Indirect Effect in 1-1-1 Mediation Models
    1. Bootstrapping
    2. Bayesian Estimation
  4. Conclusion

Prelim - Loading libraries used in this script.

library(psych)
library(ggplot2)
library(reshape2)
library(nlme)
library(lme4)
library(glmmTMB)
library(brms)

# install.packages("devtools")
# library(devtools)
# install_github("marklhc/bootmlm")
library(bootmlm) #to extract asymptotic covariance matrix of random effects
library(multilevelmediation) # from Falk, Vogel, Hammami, & Miočević (2024)
library(bmlm) # from Vuorre & Bolger (2018)

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_{ij} = d_{M_j} + a_{j}X_{ij} + \varepsilon_{M_{ij}}\]

\[Y_{ij} = d_{Y_j} + b_{j}M_{ij} + c'_{j}X_{ij} + \varepsilon_{Y_{ij}}\] As is typical of multilevel modeling, here we use \(j\) to index level 2 units or clusters (e.g., people) and \(i\) to index level 1 units (e.g., repeated measures). This setup is often referred to as the 1-1-1 model because \(X\), \(X\), and \(Y\) are measured at level 1 (hence the presence of both \(i\) and \(j\) subscripts for each of these variables). Note also that the paths of the mediation analysis model, \(a_j\), \(b_j\) and \(c_j'\), also have \(j\) subscripts. This is because these path coefficients have random effects; that is, the effects may be different for each level 2 unit or cluster. As part of any competent modeling approach, the variances of \(a_j\), \(b_j\), and \(c_j'\) will be estimated (i.e., the random effects) as well as their covariances. The covariance between \(a_j\) and \(b_j\), or \(\sigma_{a_jb_j}\), becomes important later. But, researchers are often interested in interpreting fixed effects, that is, something akin to an average effect for each of these coefficients as their interpretation is of a similar flavor to regression coefficients. For instance Bolger and Laurenceau (2013) define \(\overline{a}_j = a\), \(\overline{b}_j = b\), \(\overline{c}_j' = c'\) as the paths of the mediation analysis model. Note, however, that that typical indirect effect when both \(a\) and \(b\) paths have random effects is not the simple product of these two. Instead, we have:

2. An Example using Multivariate Multilevel Model

We follow the example from Bolger & Laurenceau (2013) Chapter 9.

Data, descriptives, plots

Read in the data

Note: In a previous version of this tutorial, we loaded a .csv data file (B&Lmediation.csv) that already had centered data.

The data now exist as the BLch9 dataset in the R package bmlm:

data(BLch9)
data <- BLch9

If the variables of interest were not already centered appropriately, the isolate function from bmlm can accomplish this. For example, here we person-center three variables:

data <- isolate(data, by = "id", value = "fwkstrs")
data <- isolate(data, by = "id", value = "fwkdis")
data <- isolate(data, by = "id", value = "freldis")

And time also appears as centered in later analyses:

# some unnecessary attributes added unless wrapped in as.vector()
data$timec <- as.vector(scale(data[,"time"], center=TRUE, scale=FALSE))

The variables of interest in this example are fwkstrs_cw (the predictor: female work stress centered within),fwkdis_cw (the mediator: female work dissatisfaction centered within),freldis_cw (the outcome: female relationship dissatisfaction centered within). Note that these variables are state variables. We are not currently interested in between-person differences, so they have not been computed for this example, but isolate can do that as well (add an argument which = "both" or which = "between").

Though it turns out that 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. Take a look and compare a few values:

head(data)
##    id time fwkstrs   fwkdis  freldis          x          m          y
## 1 101    1       3 5.590119 3.034483  0.3333333  0.9828781 -1.4420726
## 2 101    2       3 5.535224 4.620690  0.3333333  0.9279833  0.1441343
## 3 101    3       3 3.888381 2.850575  0.3333333 -0.7188603 -1.6259807
## 4 101    4       4 5.352242 6.398467  1.3333333  0.7450007  1.9219121
## 5 101    5       1 4.483074 2.544061 -1.6666667 -0.1241668 -1.9324941
## 6 101    6       2 3.339433 5.164751 -0.6666667 -1.2678081  0.6881956
##   fwkstrs_cw  fwkdis_cw freldis_cw timec
## 1  0.3333333  0.9828781 -1.4420726   -10
## 2  0.3333333  0.9279833  0.1441343    -9
## 3  0.3333333 -0.7188603 -1.6259807    -8
## 4  1.3333333  0.7450007  1.9219121    -7
## 5 -1.6666667 -0.1241668 -1.9324941    -6
## 6 -0.6666667 -1.2678081  0.6881956    -5

Descriptives of Within-person Variables

Lets look at the descriptives.

#variables of interest
vars <- c("fwkstrs_cw","fwkdis_cw","freldis_cw","x","m","y")
#descriptives
describe(data[,vars])
##            vars    n mean   sd median trimmed  mad   min  max range skew
## fwkstrs_cw    1 2100    0 1.00   0.00       0 1.20 -2.90 2.95  5.86 0.05
## fwkdis_cw     2 2100    0 1.13  -0.01       0 1.15 -3.81 3.81  7.63 0.02
## freldis_cw    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
## fwkstrs_cw    -0.25 0.02
## fwkdis_cw     -0.11 0.02
## freldis_cw     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 plot each leg of the model for a few individuals.

x –> y fwkstr_cw –> freldis_cw

#within-person regressions, x --> y
ggplot(data=data[which(data$id <= 106),], aes(x=fwkstrs_cw, y=freldis_cw, group=id)) +
  geom_point(color="black") + 
  geom_smooth(method=lm, se=TRUE, fullrange=TRUE, color="black") + 
  xlab("x = Work Stress (within)") +
  ylab("y = Relationship Dissatisfaction (within)") +
  facet_wrap( ~ id)
## `geom_smooth()` using formula = 'y ~ x'

x –> m fwkstr_cw –> fwkdis_cw

#within-person regressions, x --> m
ggplot(data=data[which(data$id <= 106),], aes(x=fwkstrs_cw, y=fwkdis_cw, group=id)) +
  geom_point(color="blue") + 
  geom_smooth(method=lm, se=TRUE, fullrange=TRUE, color="blue") + 
  xlab("x = Work Stress (within)") +
  ylab("m = Work Dissatisfaction (within)") +
  facet_wrap( ~ id)
## `geom_smooth()` using formula = 'y ~ x'

m –> y fwkdis_cw –> freldis_cw

#within-person regressions, m --> y
ggplot(data=data[which(data$id <= 106),], aes(x=fwkdis_cw, y=freldis_cw, group=id)) +
  geom_point(color="red") + 
  geom_smooth(method=lm, se=TRUE, fullrange=TRUE, color="red") + 
  xlab("m = Work Dissatisfaction (within)") +
  ylab("y = Relationship Dissatisfaction (within)") +
  facet_wrap( ~ id)
## `geom_smooth()` using formula = 'y ~ x'

Restructuring the Data Set into Double-Record data.

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

Practically, fitting of the bivariate outcome model is facilitated through some data-restructuring, specifically “double-entry” data. More details can be found in the tutorial on the multivariate multilevel model or in Bauer et al (2006).

Two ways of restructuring are shown below.

Restructuring manually is done in just a few steps, and there are multiple different packages that could be used. Here, reshape2 is used for data restructuring. tidyverse packages could also be employed (e.g., tidyr and pivot_longer).

  1. Subset down to the variables of interest.
data1 <- data[,c("id","time","timec","fwkstrs_cw","fwkdis_cw", "freldis_cw","x","m","y")]
#look at  data set
head(data1)
##    id time timec fwkstrs_cw  fwkdis_cw freldis_cw          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
##            y
## 1 -1.4420726
## 2  0.1441343
## 3 -1.6259807
## 4  1.9219121
## 5 -1.9324941
## 6  0.6881956
tail(data1)
##       id time timec  fwkstrs_cw   fwkdis_cw freldis_cw           x           m
## 2095 200   16     5 -0.95238095 -1.11183723 -0.5148695 -0.95238095 -1.11183723
## 2096 200   17     6  0.04761905 -0.06883632  0.8491151  0.04761905 -0.06883632
## 2097 200   18     7 -0.95238095  0.01350586 -0.8060573 -0.95238095  0.01350586
## 2098 200   19     8 -1.95238095  1.68779680  0.2437511 -1.95238095  1.68779680
## 2099 200   20     9 -0.95238095 -1.90781162 -0.9516512 -0.95238095 -1.90781161
## 2100 200   21    10  1.04761905  0.70883980 -1.5187010  1.04761905  0.70883980
##               y
## 2095 -0.5148695
## 2096  0.8491151
## 2097 -0.8060573
## 2098  0.2437511
## 2099 -0.9516512
## 2100 -1.5187010
  1. Melt the two outcome variables (y and m) into a single z variable using the melt() function in the reshape2 library.
#melting data
names(data1)
## [1] "id"         "time"       "timec"      "fwkstrs_cw" "fwkdis_cw" 
## [6] "freldis_cw" "x"          "m"          "y"
datalong <- melt(data=data1,
               id.vars=c("id","time","timec","fwkstrs_cw","fwkdis_cw","freldis_cw","x"),
               na.rm=FALSE, variable.name="dv",
               value.name="z")
#look at updated data set
head(datalong)
##    id time timec fwkstrs_cw  fwkdis_cw freldis_cw          x dv          z
## 1 101    1   -10  0.3333333  0.9828781 -1.4420726  0.3333333  m  0.9828781
## 2 101    2    -9  0.3333333  0.9279833  0.1441343  0.3333333  m  0.9279833
## 3 101    3    -8  0.3333333 -0.7188603 -1.6259807  0.3333333  m -0.7188603
## 4 101    4    -7  1.3333333  0.7450007  1.9219121  1.3333333  m  0.7450007
## 5 101    5    -6 -1.6666667 -0.1241668 -1.9324941 -1.6666667  m -0.1241668
## 6 101    6    -5 -0.6666667 -1.2678081  0.6881956 -0.6666667  m -1.2678081
tail(datalong)
##       id time timec  fwkstrs_cw   fwkdis_cw freldis_cw           x dv
## 4195 200   16     5 -0.95238095 -1.11183723 -0.5148695 -0.95238095  y
## 4196 200   17     6  0.04761905 -0.06883632  0.8491151  0.04761905  y
## 4197 200   18     7 -0.95238095  0.01350586 -0.8060573 -0.95238095  y
## 4198 200   19     8 -1.95238095  1.68779680  0.2437511 -1.95238095  y
## 4199 200   20     9 -0.95238095 -1.90781162 -0.9516512 -0.95238095  y
## 4200 200   21    10  1.04761905  0.70883980 -1.5187010  1.04761905  y
##               z
## 4195 -0.5148695
## 4196  0.8491151
## 4197 -0.8060573
## 4198  0.2437511
## 4199 -0.9516512
## 4200 -1.5187010
  1. Adding dummy variables that indicate row entry (and reordering).
#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)

#reordering for convenience 
datalong <- datalong[order(datalong$id,datalong$timec,datalong$dm),]

#look at updated data set
head(datalong)
##       id time timec fwkstrs_cw  fwkdis_cw freldis_cw         x dv          z dy
## 2101 101    1   -10  0.3333333  0.9828781 -1.4420726 0.3333333  y -1.4420726  1
## 1    101    1   -10  0.3333333  0.9828781 -1.4420726 0.3333333  m  0.9828781  0
## 2102 101    2    -9  0.3333333  0.9279833  0.1441343 0.3333333  y  0.1441343  1
## 2    101    2    -9  0.3333333  0.9279833  0.1441343 0.3333333  m  0.9279833  0
## 2103 101    3    -8  0.3333333 -0.7188603 -1.6259807 0.3333333  y -1.6259807  1
## 3    101    3    -8  0.3333333 -0.7188603 -1.6259807 0.3333333  m -0.7188603  0
##      dm dvnum
## 2101  0     0
## 1     1     1
## 2102  0     0
## 2     1     1
## 2103  0     0
## 3     1     1
tail(datalong)
##       id time timec fwkstrs_cw  fwkdis_cw freldis_cw         x dv          z dy
## 4198 200   19     8  -1.952381  1.6877968  0.2437511 -1.952381  y  0.2437511  1
## 2098 200   19     8  -1.952381  1.6877968  0.2437511 -1.952381  m  1.6877968  0
## 4199 200   20     9  -0.952381 -1.9078116 -0.9516512 -0.952381  y -0.9516512  1
## 2099 200   20     9  -0.952381 -1.9078116 -0.9516512 -0.952381  m -1.9078116  0
## 4200 200   21    10   1.047619  0.7088398 -1.5187010  1.047619  y -1.5187010  1
## 2100 200   21    10   1.047619  0.7088398 -1.5187010  1.047619  m  0.7088398  0
##      dm dvnum
## 4198  0     0
## 2098  1     1
## 4199  0     0
## 2099  1     1
## 4200  0     0
## 2100  1     1

Some R packages to automate data restructuring in this way have begun to appear. For example, stack_bpg from the multilevelmediation package also accomplishes this restructuring, although column names are changed to work nicely with internal functions for the package.

datlong2 <- stack_bpg(data, L2ID = "id", X = "x", Y = "y", M = "m",
                      covars.m = c("time","timec"),
                      covars.y = c("time","timec"))
head(datlong2) # compare to datlong
## # A tibble: 6 × 12
##       X  L2id     Md  time timec Outcome      Z    Sy    Sm   SmX   SyX    SyM
##   <dbl> <int>  <dbl> <int> <dbl> <chr>    <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1 0.333   101  0.983     1   -10 Y       -1.44      1     0 0     0.333  0.983
## 2 0.333   101  0.983     1   -10 M        0.983     0     1 0.333 0      0    
## 3 0.333   101  0.928     2    -9 Y        0.144     1     0 0     0.333  0.928
## 4 0.333   101  0.928     2    -9 M        0.928     0     1 0.333 0      0    
## 5 0.333   101 -0.719     3    -8 Y       -1.63      1     0 0     0.333 -0.719
## 6 0.333   101 -0.719     3    -8 M       -0.719     0     1 0.333 0      0

Here, Outcome is a readable column that tells which outcome is being modeled for each row with Z containing the actual values of the outcome. Sy and Sm are indicators that will take the place of the intercepts in the same way as dy and dm in our manual restructuring. And SmX, SyX and SyM take the place of \(X\) when \(M\) is the outcome, and \(X\) and \(M\) when \(Y\) is the outcome, respectively. Note how when \(M\) is the outcome, SyX and SyM are already zero (since these two variables are only relevant when \(Y\) is the outcome), and when \(Y\) is the outcome, SmX has zero (since this variable is only relevant when \(M\) is the outcome). The same can be accomplished with the first data restructuring strategy if when entering predictors into the model, each predictor is multiplied by dy or dm as necessary.

Estimating the Model Several Ways

lme

Here we use lme from the nlme package because it can appropriately model different error variances for each outcome (i.e., heteroscedasticity). Those who are familiar with lmer from lme4 take note that lmer cannot do this.

Note also that dy and dm take the place of a single intercept, we can therefore force the function to not add some other intercept by specifying -1 or 0 in the model equation. Interestingly, on p. 185 of Bolger & Laurenceau (2013) there is mention of omitting dy and dm altogether as the intercepts will be zero because \(M\) and \(Y\) were person-mean centered; we include them here for now as it does not hurt and some such intercepts show up in the output on the next page (p. 186). Relatedly, although one would usually estimate random effects for the intercepts (Bauer, Preacher, and Gil , 2006)]), in this case the random effect variances will also be zero and are thus omitted.

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

#lme mediation model
model.lme <- lme(fixed = z ~ 0 + dm + dy +  # intercepts
                             dm:fwkstrs_cw + dm:timec + # prediction of mediator 
                             dy:fwkdis_cw + dy:fwkstrs_cw + dy:timec, # prediction of outcome
                  random = ~ 0  +  dm:fwkstrs_cw + dy:fwkdis_cw + dy:fwkstrs_cw | id, # random effects for all three paths
                  weights = varIdent(form = ~ 1 | dvnum), #this invokes 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.62 12290.74 -6082.811
## 
## Random effects:
##  Formula: ~0 + dm:fwkstrs_cw + dy:fwkdis_cw + dy:fwkstrs_cw | id
##  Structure: General positive-definite, Log-Cholesky parametrization
##               StdDev     Corr         
## dm:fwkstrs_cw 0.25988626 dm:fw_ dy:fw_
## dy:fwkdis_cw  0.21720775 0.552        
## fwkstrs_cw:dy 0.08791889 0.448  0.939 
## Residual      0.92555835              
## 
## Variance function:
##  Structure: Different standard deviations per stratum
##  Formula: ~1 | dvnum 
##  Parameter estimates:
##        0        1 
## 1.000000 1.177733 
## Fixed effects:  z ~ 0 + dm + dy + dm:fwkstrs_cw + dm:timec + dy:fwkdis_cw + dy:fwkstrs_cw +      dy:timec 
##                     Value  Std.Error   DF   t-value p-value
## dm             0.00000000 0.02378708 4094  0.000000  1.0000
## dy             0.00000000 0.02019734 4094  0.000000  1.0000
## dm:fwkstrs_cw  0.19031195 0.03551311 4094  5.358921  0.0000
## dm:timec      -0.00583019 0.00397781 4094 -1.465677  0.1428
## dy:fwkdis_cw   0.14775645 0.02862813 4094  5.161234  0.0000
## dy:fwkstrs_cw  0.10532275 0.02266911 4094  4.646091  0.0000
## dy:timec      -0.00247490 0.00338825 4094 -0.730436  0.4652
##  Correlation: 
##               dm     dy     dm:fw_ dm:tmc dy:fwkd_ dy:fwks_
## dy             0.000                                       
## dm:fwkstrs_cw  0.000  0.000                                
## dm:timec       0.000  0.000 -0.012                         
## dy:fwkdis_cw   0.000  0.000  0.305  0.001                  
## dy:fwkstrs_cw  0.000  0.000  0.129  0.001  0.192           
## dy:timec       0.000  0.000  0.000  0.003  0.015   -0.014  
## 
## Standardized Within-Group Residuals:
##          Min           Q1          Med           Q3          Max 
## -3.547907404 -0.674772174 -0.006486992  0.669579417  3.349772750 
## 
## Number of Observations: 4200
## Number of Groups: 100

The results look good and closely match the results in the book! Yay!

The same can be done with the other version of the restructured data:

model.lme2 <- lme(fixed = Z ~ 0 + Sm + Sy + # intercepts
                             SmX + Sm:timec + # prediction of mediator 
                             SyM + SyX + Sy:timec, # prediction of outcome
                  random = ~ 0  + SmX + SyM + SyX | L2id, # random effects
                  weights = varIdent(form = ~ 1 | Sm), #this invokes separate sigma^{2}_{e} for each outcome
                  data=datlong2,
                  na.action=na.exclude,
                  control=lmeControl(opt="optim",maxIter = 200, msMaxIter = 200, niterEM = 50, msMaxEval = 400))
summary(model.lme2)
## Linear mixed-effects model fit by REML
##   Data: datlong2 
##        AIC      BIC    logLik
##   12195.62 12290.74 -6082.811
## 
## Random effects:
##  Formula: ~0 + SmX + SyM + SyX | L2id
##  Structure: General positive-definite, Log-Cholesky parametrization
##          StdDev     Corr       
## SmX      0.25988626 SmX   SyM  
## SyM      0.21720776 0.552      
## SyX      0.08791889 0.448 0.939
## Residual 0.92555835            
## 
## Variance function:
##  Structure: Different standard deviations per stratum
##  Formula: ~1 | Sm 
##  Parameter estimates:
##        0        1 
## 1.000000 1.177733 
## Fixed effects:  Z ~ 0 + Sm + Sy + SmX + Sm:timec + SyM + SyX + Sy:timec 
##                Value  Std.Error   DF   t-value p-value
## Sm        0.00000000 0.02378708 4094  0.000000  1.0000
## Sy        0.00000000 0.02019734 4094  0.000000  1.0000
## SmX       0.19031195 0.03551311 4094  5.358921  0.0000
## SyM       0.14775645 0.02862813 4094  5.161234  0.0000
## SyX       0.10532275 0.02266911 4094  4.646091  0.0000
## Sm:timec -0.00583019 0.00397781 4094 -1.465677  0.1428
## Sy:timec -0.00247490 0.00338825 4094 -0.730436  0.4652
##  Correlation: 
##          Sm     Sy     SmX    SyM    SyX    Sm:tmc
## Sy        0.000                                   
## SmX       0.000  0.000                            
## SyM       0.000  0.000  0.305                     
## SyX       0.000  0.000  0.129  0.192              
## Sm:timec  0.000  0.000 -0.012  0.001  0.001       
## Sy:timec  0.000  0.000  0.000  0.015 -0.014  0.003
## 
## Standardized Within-Group Residuals:
##          Min           Q1          Med           Q3          Max 
## -3.547907418 -0.674772171 -0.006486992  0.669579409  3.349772756 
## 
## Number of Observations: 4200
## Number of Groups: 100

The results here are identical, but just different labels for the coefficients. For example, compare the two Tables:

coefficient estimate
dm 0.0000000
dy 0.0000000
dm:fwkstrs_cw 0.1903119
dy:fwkdis_cw 0.1477565
dy:fwkstrs_cw 0.1053227
dm:timec -0.0058302
dy:timec -0.0024749
coefficient estimate
Sm 0.0000000
Sy 0.0000000
SmX 0.1903119
SyM 0.1477565
SyX 0.1053227
Sm:timec -0.0058302
Sy:timec -0.0024749

With the correspondence as follows:

  • dm or Sm - intercept when \(M\) is the outcome, \(d_{M_j}\)
  • dy or Sy - intercept when \(Y\) is the outcome, \(d_{Y_j}\)
  • dm:fwkstrs_cw or SmX - effect of \(X\) when \(M\) is the outcome (i.e., \(a\) path), \(a_j\)
  • dm:fwkdis_cw or SyM - effect of \(M\) when \(Y\) is the outcome (i.e., \(b\) path), \(b_j\)
  • dm:fwkstrs_cw or SyX - effect of \(X\) when \(Y\) is the outcome (i.e., direct effect), \(c'_j\)
  • dm:timec or Sm:timec - effect of covariate (timec) on \(M\)
  • dy:timec or Sy:timec - effect of covariate (timec) on \(Y\)

lmer

Note that in a previous version of this tutorial, lmer from the lme package was also demonstrated based in part on this tutorial: https://stats.idre.ucla.edu/r/faq/how-can-i-perform-mediation-with-multilevel-data-method-2/ But, generally lmer cannot work with heterogeneous error variances and is now omitted.

glmmTMB

If one wants something more modern than nlme, but more flexible (i.e., supports heteroscedasticity), then glmmTMB may be an option. Let’s see if we can estimate the same model.

fitglmmtmb <- glmmTMB(Z ~ 0 + Sm + Sy +
                        SmX + Sm:timec +
                        SyX + SyM + Sy:timec + 
                     (0 + SmX + SyM + SyX | L2id),
                   dispformula =  ~ 0 + Sm + Sy,
                   family=gaussian,
                   data=datlong2,
                   na.action=na.exclude,
                   control = glmmTMBControl(optimizer = optim, optArgs=list(method="BFGS")),
                   REML = TRUE)

It works! (Apparently after fiddling with the optimization options a bit). And results are again pretty close to the results obtained previously:

summary(fitglmmtmb)
##  Family: gaussian  ( identity )
## Formula:          
## Z ~ 0 + Sm + Sy + SmX + Sm:timec + SyX + SyM + Sy:timec + (0 +  
##     SmX + SyM + SyX | L2id)
## Dispersion:         ~0 + Sm + Sy
## Data: datlong2
## 
##      AIC      BIC   logLik deviance df.resid 
##  12195.6  12290.8  -6082.8  12165.6     4185 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name Variance Std.Dev. Corr      
##  L2id     SmX  0.067960 0.26069            
##           SyM  0.047187 0.21723  0.55      
##           SyX  0.008022 0.08957  0.38 0.92 
##  Residual            NA      NA            
## Number of obs: 4200, groups:  L2id, 100
## 
## Conditional model:
##            Estimate Std. Error z value Pr(>|z|)    
## Sm        2.625e-12  2.378e-02   0.000    1.000    
## Sy        6.612e-13  2.020e-02   0.000    1.000    
## SmX       1.903e-01  3.559e-02   5.348 8.88e-08 ***
## SyX       1.048e-01  2.337e-02   4.484 7.31e-06 ***
## SyM       1.488e-01  2.911e-02   5.110 3.22e-07 ***
## Sm:timec -5.837e-03  3.981e-03  -1.466    0.143    
## Sy:timec -2.448e-03  3.392e-03  -0.722    0.470    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Dispersion model:
##    Estimate Std. Error z value Pr(>|z|)    
## Sm  0.08600    0.01582   5.435 5.48e-08 ***
## Sy -0.07731    0.01584  -4.880 1.06e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

It would probably help if future work ironed out how well glmmTMB performed versus nlme here. Note that diagnosing and troubleshooting model estimation issues is an ongoing problem: https://cran.r-project.org/web/packages/glmmTMB/vignettes/troubleshooting.html

multilevelmediation

The multilevelmediation package contains a utility function that combines both restructuring and model fitting in one step. It also recently offers support for model fitting using either nlme (default) or glmmTMB. Let’s see if it also can yield the same results. If so, it can be used for bootstrapping for inferences and confidence intervals. One difference, at least at the time of writing, involves many arguments to modmed.mlm instead of formula-like syntax.

In brief,

  • L2ID allows defining the level 2 ID variable in the dataset.
  • X, Y, M, covars.m and covars.y arguments define the variables to include in the analysis. Covariates (here, time centered) are optional.
  • random.a, random.b, random.cprime, random.int.m, and random.int.y are boolean arguments and define whether random effects should be present for the a-path, b-path, c’-path and the intercepts.
  • control are options passed to lme or whatever function is used to estimate the data.
fitmm <- modmed.mlm(data1, L2ID = "id", X = "x", Y = "y", M = "m",
           covars.m = "timec", covars.y = "timec",
           random.a = TRUE, random.b = TRUE, random.cprime = TRUE,
           random.int.m = FALSE, random.int.y = FALSE,
           control=lmeControl(opt="optim",maxIter = 200, msMaxIter = 200, niterEM = 50, msMaxEval = 400))

summary(fitmm$model)
## Linear mixed-effects model fit by REML
##   Data: tmp 
##        AIC      BIC    logLik
##   12195.62 12290.74 -6082.811
## 
## Random effects:
##  Formula: ~0 + SmX + SyM + SyX | L2id
##  Structure: General positive-definite, Log-Cholesky parametrization
##          StdDev     Corr       
## SmX      0.25988626 SmX   SyM  
## SyM      0.21720775 0.552      
## SyX      0.08791888 0.448 0.939
## Residual 0.92555835            
## 
## Variance function:
##  Structure: Different standard deviations per stratum
##  Formula: ~1 | Sm 
##  Parameter estimates:
##        0        1 
## 1.000000 1.177733 
## Fixed effects:  as.formula(fixed.formula) 
##                Value  Std.Error   DF   t-value p-value
## Sm        0.00000000 0.02378708 4094  0.000000  1.0000
## Sy        0.00000000 0.02019734 4094  0.000000  1.0000
## SmX       0.19031195 0.03551311 4094  5.358921  0.0000
## SyX       0.10532275 0.02266911 4094  4.646091  0.0000
## SyM       0.14775645 0.02862813 4094  5.161234  0.0000
## Sm:timec -0.00583019 0.00397781 4094 -1.465677  0.1428
## Sy:timec -0.00247490 0.00338825 4094 -0.730436  0.4652
##  Correlation: 
##          Sm     Sy     SmX    SyX    SyM    Sm:tmc
## Sy        0.000                                   
## SmX       0.000  0.000                            
## SyX       0.000  0.000  0.129                     
## SyM       0.000  0.000  0.305  0.192              
## Sm:timec  0.000  0.000 -0.012  0.001  0.001       
## Sy:timec  0.000  0.000  0.000 -0.014  0.015  0.003
## 
## Standardized Within-Group Residuals:
##          Min           Q1          Med           Q3          Max 
## -3.547907410 -0.674772173 -0.006486992  0.669579411  3.349772752 
## 
## Number of Observations: 4200
## Number of Groups: 100

Results look similar to before! In particular, the results are essentially the same as that in model.lme, but using different notation and with the effects listed in the output in a slightly different order. The notation is the same as that in model.lme2.

Let’s look more closely at the results.

A Closer Look at Results

The fixed effects …

We are working from model.lme, but recall we could also look at model.lme2, fitglmmtmb, or fitmm.

  • dm or Sm
  • dy or Sy
  • dm:fwkstrs_cw or SmX
  • dm:fwkdis_cw or SyM
  • dm:fwkstrs_cw or SyX
  • dm:timec or Sm:timec
  • dy:timec or Sy:timec
#fixed effects
round(FEmatrix <- fixef(model.lme), 4)
##            dm            dy dm:fwkstrs_cw      dm:timec  dy:fwkdis_cw 
##        0.0000        0.0000        0.1903       -0.0058        0.1478 
## dy:fwkstrs_cw      dy:timec 
##        0.1053       -0.0025

Lets interpret …
0 = intercept in the m as outcome model (We expect this to be zero because we person-centered all the data) 0 = intercept in the y as outcome model (We expect this to be zero because we person-centered all the data) 0.1903 = effect of x –> m (a: slope of work stressors predicting work dissatisfaction)
-0.0058 = effect of time –> m (time trend in work dissatisfaction)
0.1478= effect of m –> y (b: slope of work dissatisfaction predicting relationship dissatisfaction)
0.1053 = effect of x –> y (c’: slope of work stressors predicting relationship dissatisfaction, after adjusting for work dissatisfaction)
-0.0025 = effect of time –> y (time trend in relationship dissatisfaction)

Lets put the parameters into named objects, as these will be useful later when showing how calculations for the indirect effect work out.

#making parameter objects
a <- FEmatrix[3]
b <- FEmatrix[5]
cprime <- FEmatrix[6]

The random effects …

The random effects can be output this way though are also visible when using summary() on any given model object. Note that depending on the function used and how extracted, sometimes the standard deviation of random effects are given, sometimes the variance is given.

(REmatrix <- getVarCov(model.lme))
## Random effects variance covariance matrix
##               dm:fwkstrs_cw dy:fwkdis_cw fwkstrs_cw:dy
## dm:fwkstrs_cw      0.067541     0.031139     0.0102360
## dy:fwkdis_cw       0.031139     0.047179     0.0179360
## fwkstrs_cw:dy      0.010236     0.017936     0.0077297
##   Standard Deviations: 0.25989 0.21721 0.087919

They are a little easier to interpret and work with in this form …

The variance of the a paths, \(\sigma^2_{a}\),is 0.0675409

sig2_a <- REmatrix[1,1]

The variance of the b paths, \(\sigma^2_{b}\),is 0.0471792

sig2_b <- REmatrix[2,2]

The variance of the c’ paths, \(\sigma^2_{c'}\),is 0.0077297

sig2_cprime <- REmatrix[3,3]

The residual variance of the outcome and mediator variable are less obvious and the result varies a bit depending on whether nlme or glmmTMB is used, as well as how the command specifying heteroscedasticity is defined.

Consider output (look way back) from summary(model.lme) when the nlme package was used to estimate the model. In the random effects section, there is a residual SD listed as: 0.9255583. That’s the residual SD for one of the outcomes; to know which, look at the variance function from that output and how we defined the residual structure: Formula: ~1 | dvnum. If you look at dvnum in the data, it has a “1” when M is the outcome, and 0 when Y is the outcome. The coding works similarly to dummy coding; the residual SD listed here is when dvnum is equal to 0. So, it turns out that 0.9255583 is the standard deviation when dvnum is 0, or for Y.

Therefore residual variance of the outcome variable, \(\sigma^2_{ey}\), is given by “Residual” = 0.8566583;

Now look further down under the variance function of the output. There are some parameter estimates:

model.lme$modelStruct$varStruct
## Variance function structure of class varIdent representing
##        0        1 
## 1.000000 1.177733

Note the second value, 1.1777333. If we do, 0.9255583 \(\times\) 1.1777333, we end up with the residual SD for M, 1.0900609.

Therefore, the residual variance of the mediator variable, \(\sigma^2_{em}\), is 1.1882327.

In R code, it is possible to automate computing of both of the above…

model.lme$sigma^2 # residual variance for Y
## [1] 0.8566583
(model.lme$sigma*1/attr(model.lme$modelStruct$varStruct,"weights")[2])^2 # residual variance for M
##        1 
## 1.188233

Both of these values are fairly close to that provided in the book. Note that a similar logic can be applied to obtain the residual variances under model.lme2; in fact, there the exact same values are provided.

Finally… the covariance between the \(a_{j}\) and \(b_{j}\) paths, \(\sigma_{ajbj}\) is 0.0311386

covajbj <- REmatrix[1,2]

Output from glmmTMB is slightly easier to work with for obtaining residual variances ; documentation appears to suggest that the Dispersion model output contains the residual variances, but these SDs are represented on a log scale. Converting them to variances is actually not that hard. We can extract them in the following way:

resids <- fixef(fitglmmtmb)$disp
resids
##          Sm          Sy 
##  0.08599913 -0.07731211

Then attempt to convert them to variances:

(exp(resids))^2
##        Sm        Sy 
## 1.1876758 0.8567371

Hey, those are also really close to the output given in the book.

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

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

ie = ab + \(\sigma_{ajbj}\)
= 0.1903119*0.1477565 + 0.0311386 = 0.059

indirecteffect <- a*b + covajbj

Calculating total effect c

c = c’ + ab + \(\sigma_{ajbj}\) = 0.1053227 + 0.1903119*0.1477565 + 0.0311386 = 0

totaleffect <- cprime + a*b + covajbj

Calculating percent of total effect accounted for by mediated effect

=(ab + $\sigma_{ajbj}$) / (c' + ab + $\sigma_{ajbj}$)    
= ie / c    
= 0.0592584 / 0.1645812 
  = 0.3600559
percentmediated <- 100*(indirecteffect/totaleffect)

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

= $\sigma_{ajbj}$ / (ab + $\sigma_{ajbj}$)    
= $\sigma_{ajbj}$ / ie    
= 0.0311386 / 0.0592584
= 0.5254713
percentcovariance <- 100*(covajbj/indirecteffect)

All these calculations match the results in B&L’s Chapter 9 pretty well.

  1. Obtaining confidence (or credibility) intervals for the indirect effect

Obtaining Confidence (or Credibility) Intervals for the Indirect Effect in 1-1-1 Mediation Models

We also often want to obtain confidence intervals for the indirect effect, ab + \(\sigma_{ajbj}\) = 0.059. This directly provides us with some information regarding the uncertainty in the point estimate and is often used for making inferences. For example, any 95% confidence interval that does not include 0 is often taken as rejection of the null hypothesis that the indirect effet is nil at the \(\alpha=.05\) level.

There are several methods for testing Mediation effects for single level Mediation models. Many of these have been implemented in R packages (e.g., mediation). As computational speed has increased, confidence intervals for indirect effects are now often obtained using Monte Carlo, bootstrapping, or Bayesian approachs. Some of these approaches are being incorporated into recent R packages (e.g., multilevelmediation). In the current tutorial, we present bootstrapping and Bayesian estimation. A forthcoming revision shall include the Monte Carlo method. Note that in a previous version of this tutorial, lmer from the lme4 package was used with the Monte Carlo method. However, this approach did not properly model heteroscedasticity and has been removed.

Bootstrapping

One way to obtain CIs is to conduct nonparametric bootstrapping. This approach has been popular for quite some time with single-level models, but has only recently been studied for multilevel models (Falk, Vogel, Hammami, & Miočević, 2024).

Note that a previous version of this tutorial mentioned a convenience function indirectMLM.R, which would do bootstrapping. Falk et al note (footnote 3, p. 758) that this tool did not perform well in some preliminary simulations. It is therefore not recommended for use.

The multilevelmediation package instead can do some bootstrapping methods. Included are:

  • Case resampling at level 2
    • Only resample clusters or level 2 units
  • Case resampling at level 1
    • Only resample level 1 units (observations within clusters)
  • Case resampling of both level 2 and level 1
    • Resample level 2 units, then level 1 observations within each
  • Residual-based bootstrap
    • Resamples residuals at level 2 and level 1

Which approach works well depends a bit on the data and presence of certain random effects. Case resampling only at level 1 fails quite miserably if certain random effects are present, and probably is not a safe bet. Case resampling at level 2 or both level 2 and 1 worked fairly well in simulations, but are best if cluster sizes are about equal. If cluster sample size (i.e., number of level 1 units) is not equal, then the residual-based bootstrap is probably a good choice.

To do bootstrapping, boot.modmed.mlm.custom from the multilevelmediation package can be used. The syntax is almost the same as the modmed.mlm function:

boot.result<-boot.modmed.mlm.custom(data1, L2ID = "id", X = "x", Y = "y", M = "m",
           covars.m = "timec", covars.y = "timec",
           random.a = TRUE, random.b = TRUE, random.cprime = TRUE,
           random.int.m = FALSE, random.int.y = FALSE,
           boot.type="caseboth", nrep = 1000,
           parallel.type="parallel",ncores=8,seed=2299,
           control=lmeControl(opt="optim",maxIter = 200, msMaxIter = 200, niterEM = 50, msMaxEval = 400)
           )

Additional options here include:

  • boot.type - "caseboth" resamples at both levels; other options are "case2", "case1", or "resid" for the other kinds of bootstrap
  • nrep - number of bootstrap replications. Set to something large like 1000 or more for more stable intervals, though larger number will mean it takes longer
  • parallel.type and ncores allow one to do parallel processing; here, we use 8 processing cores and the parallel package.

One issue with bootstrapping multilevel models is that these are often quite slow to estimate, so please be patient. Once finished, we can extract confidence intervals for the indirect effect and some of the other paths in the model:

extract.boot.modmed.mlm(boot.result, type="indirect", ci.conf=.95)
## $CI
##       2.5%      97.5% 
## 0.02751116 0.09207604 
## 
## $est
## [1] 0.05925842
extract.boot.modmed.mlm(boot.result, type="a", ci.conf=.95)
## $CI
##      2.5%     97.5% 
## 0.1157208 0.2589108 
## 
## $est
## [1] 0.1903119
extract.boot.modmed.mlm(boot.result, type="b", ci.conf=.95)
## $CI
##       2.5%      97.5% 
## 0.08688845 0.20638547 
## 
## $est
## [1] 0.1477565
extract.boot.modmed.mlm(boot.result, type="cprime", ci.conf=.95)
## $CI
##      2.5%     97.5% 
## 0.0462796 0.1607349 
## 
## $est
## [1] 0.1053227

Bayesian Estimation

Bayesian estimation of the model can also be accomplished in a number of ways. A full introduction of Bayesian analysis is somewhat outside of the scope of this tutorial, but we provide a few references for interested readers: van de Schoot et al (2014); van de School and Depaoli (2014), and Yuan and MacKinnon (2009). The Bayesian approach represents a fundamentally different way of thinking about the truth of parameters and estimation; it allows the possibility of incorporating the analysts prior beliefs regarding the values of some parameters. In cases where there is prior information (e.g., based on past research), this approach then clearly has advantages. Use of uninformative priors (often the default) often results in point estimates and credibility intervals (similar to confidence intervals, but technically different in interpretation) that are very similar to those from a frequentist analysis. For our purposes, a Bayesian analysis allows obtaining a joint (posterior) distribution among all parameter estimates, and it is then relatively easy to obtain credibility intervals for any function of model parameters, such as the indirect effect.

bmlm

The bmlm package can do Bayesian estimation of the 1-1-1 mediation analysis model and the underlying code that executes the analysis is the popular statistical computing language Stan. The package itself is also co-authored by Bolger, one of the authors of the book we are referencing for this example. However, one current limitation of the package is that it does not yet support covariates, making inclusion of the timec variable not possible in the analysis. Regardless, let’s take a look at the typical syntax and output without this covariate (here, refresh=0 just suppresses iteration output and can be removed):

bmlm.result <- mlm(data1, id="id", x="x", m="m", y="y", refresh=0)
## Estimating model, please wait.
## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess
mlm_summary(bmlm.result)
##   Parameter Mean   SE Median 2.5% 97.5% n_eff Rhat
## 1         a 0.19 0.04   0.19 0.12  0.26  2957    1
## 2         b 0.15 0.03   0.15 0.09  0.21  2739    1
## 3        cp 0.10 0.02   0.10 0.06  0.15  5044    1
## 4        me 0.06 0.01   0.06 0.03  0.09  1539    1
## 5         c 0.16 0.03   0.16 0.11  0.21  4138    1
## 6       pme 0.36 0.08   0.35 0.21  0.54  1720    1
  • a - a path
  • b - b path
  • cp - direct effect
  • me - mediated effect or indirect effect
  • c - total effect
  • pme - percent mediated effect

For more information regarding the package and its capabilities, we refer the reader to the related publication: Vuorre & Bolger (2018) and the documentation for bmlm. Note that the package can incorporate some prior information regarding model parameters and includes additional visualizations.

multilevelmediation

On the other hand, multilevelmediation recently also supports estimation of the model using Bayesian methods and can handle covariates. Under the hood, the package uses the brms package in R Bürkner (2017), which also incidentally conducts the analysis using Stan. To conduct the analysis below, we enable parallel processing–by default 4 chains are used, so we use 4 processing cores.

Note that most of the syntax for modmed.mlm.brms is the same as that used before by modmed.mlm. Thus, no further explanation is given here.

options(mc.cores=4) # set number of processing cores

# run analysis
brms.result<-modmed.mlm.brms(data1, L2ID = "id", X = "x", Y = "y", M = "m",
           covars.m = "timec", covars.y = "timec",
           random.a = TRUE, random.b = TRUE, random.cprime = TRUE,
           random.int.m = FALSE, random.int.y = FALSE
           )

Once the model is finished estimating, the following should allow access to the brms object that is estimated:

summary(brms.result$model)
##  Family: gaussian 
##   Links: mu = identity; sigma = log 
## Formula: Z ~ 0 + Sm + Sy + SmX + SyX + SyM + Sm:timec + Sy:timec + (0 + SmX + SyM + SyX | L2id) 
##          sigma ~ 0 + Sm + Sy
##    Data: tmp (Number of observations: 4200) 
##   Draws: 4 chains, each with iter = 7000; warmup = 3500; thin = 1;
##          total post-warmup draws = 14000
## 
## Multilevel Hyperparameters:
## ~L2id (Number of levels: 100) 
##              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(SmX)          0.26      0.04     0.19     0.34 1.00     7577     9829
## sd(SyM)          0.22      0.03     0.17     0.28 1.00     8184    10493
## sd(SyX)          0.08      0.03     0.02     0.15 1.00     6155     5992
## cor(SmX,SyM)     0.51      0.15     0.19     0.79 1.00     4223     6360
## cor(SmX,SyX)     0.31      0.31    -0.36     0.85 1.00    10035     9617
## cor(SyM,SyX)     0.70      0.23     0.13     0.97 1.00    10762    10096
## 
## Regression Coefficients:
##          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Sm           0.00      0.02    -0.05     0.05 1.00    37974     9374
## Sy          -0.00      0.02    -0.04     0.04 1.00    34785     9724
## SmX          0.19      0.04     0.12     0.26 1.00    11403    10852
## SyX          0.10      0.02     0.06     0.15 1.00    22342    12090
## SyM          0.15      0.03     0.09     0.21 1.00    12704    12065
## Sm:timec    -0.01      0.00    -0.01     0.00 1.00    35095    10431
## Sy:timec    -0.00      0.00    -0.01     0.00 1.00    32061     9941
## sigma_Sm     0.09      0.02     0.06     0.12 1.00    28180     9990
## sigma_Sy    -0.08      0.02    -0.11    -0.04 1.00    27327    10310
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

When conducting Bayesian data analysis, there are also some diagnostics that should be checked. We provide a very brief overview here. Note that the bottom of this output should list some information that is useful in diagnosing any estimation problems. For example, we want Rhat values that are close to 1 (higher would indicate convergence problems, e.g., such as 1.05 or higher) and Bulk_ESS and Tail_ESS values that are also nontrivial (some recommendations put this at 100 per chain, but we probably want higher with larger models).

In addition, it is common to examine traceplots in order to visually assess convergence:

plot(brms.result$model) # print all traceplots

On the left-hand side is a histogram of the posterior for each parameter. On the right-hand side are traceplots. These are plotting draws from the posterior distribution for each chain of the analysis and plotted are the values of the parameter estimates (y-axis) across iterations (the x-axis). In brief, we usually expect unimodal posteriors and want a jumbled mess where the chains criss-cross randomly.

So long as estimation looks ok (it does here), the column Estimate from the above output provides a point estimate (mediation of posterior) for each parameter, Est.Error is similar to a standard error (posterior standard deviation), and two columns list the lower (l-95% CI) and upper (u-95% CI) bounds of a credibility interval.

multilevelmediation then also automates computation of credibility intervals for the indirect effect:

extract.modmed.mlm.brms(brms.result, "indirect")
## $CI
## # A tibble: 1 × 10
##   variable   mean median     sd    mad   q2.5  q97.5  rhat ess_bulk ess_tail
##   <chr>     <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl> <dbl>    <dbl>    <dbl>
## 1 indirect 0.0588 0.0579 0.0141 0.0138 0.0339 0.0892  1.00    6469.    9329.
## 
## $draws
## # A draws_matrix: 3500 iterations, 4 chains, and 1 variables
##     variable
## draw indirect
##   1     0.075
##   2     0.040
##   3     0.058
##   4     0.043
##   5     0.040
##   6     0.044
##   7     0.076
##   8     0.072
##   9     0.061
##   10    0.043
## # ... with 13990 more draws

We may report the median of the posterior as the point estimate for the indirect effect, 0.058, and we have a 95% credibility interval indicated by the q2.5 and q97.5 parts of the output: [round(res.indirect$CI$q2.5,3), round(res.indirect$CI$q97.5,3)].

brms

Finally, for didactic purposes, it may also be useful to see how one could custom define a brms model. Suppose we use the same data that was restructured using stack_bpg, i.e., datlong2, earlier in this tutorial.

The equations to estimate the model are encapsulated in bf() and we have split them over three lines:

  • The following lists Z as the outcome variable and all fixed effects: Z ~ 0 + Sm + Sy + SmX + SyX + SyM + Sm:timec + Sy:timec
    • Refer back to the example earlier in this tutorial using lme for an explanation of all of the terms here.
  • The next line contains the random effects: (0 + SmX + SyM + SyX | L2id ).
    • Don’t forget the + symbol between the fixed and random effects
  • This line, separated by a comma from the random effects, ensures that the errors are heteroscedastic: sigma ~ 0 + Sm + Sy
brms.result2 <- brm(formula = bf(Z ~ 0 + Sm + Sy + SmX + SyX + SyM + Sm:timec + Sy:timec  + 
                                   (0  + SmX + SyM + SyX | L2id ),
                                 sigma ~ 0 + Sm + Sy),
                         data = datlong2)

4. Conclusion

Ok - we have attempted to push into the within-person (1-1-1) mediation models. The software to make estimation of these models easy in R has advanced in the last several years since initially publishing this tutorial. In particular, there are several approaches for obtaining CIs for the indirect effect that were not easily available before. Additional work could remain in making the output more user friendly, visualizing raw data, providing effect sizes, and so on.

Finally, while Bauer, Preacher, and Gil (2006) mention moderated mediation in multilevel models, additional work is needed on how to do this in R. multilevelmediation supports some moderated mediation models, but work remains on evaluating the approach in simulations and illustrating how to conduct such an analysis.

Good Luck!