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 wihtin-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://afhayes.com/index.html, http://davidakenny.net/cm/mediate.htm, https://psychology.clas.asu.edu/research/labs/research-prevention-laboratory-mackinnon.

Outline

  1. Introduction to The Within-Person Mediation Model
  2. An Example using Multivariate Multilevel Model
  3. An Example using Multi-Step Estimation and Bootstrapping
  4. Conclusion

Prelim - Loading libraries used in this script.

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

#install.packages("devtools")
#library(devtools)
#install_github("marklhc/bootmlm")  m 
library(bootmlm) #to extract asymptotic covariance matrix of random effects

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_{jt} = a0_{j} + a_{j}X_{jt} + e_{Mjt}\] \[Y_{jt} = b0_{j} + b_{j}M_{jt} + c'_{i}X_{jt} + e_{Yjt}\]

2. An Example using Multivariate Multilevel Model

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

Preliminaries -

Read in the data

#Daily 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
## 1  -0.3033333 -0.6227591 -0.1634446  0.3333333  0.9828781 -1.4420726
## 2  -0.3033333 -0.6227591 -0.1634446  0.3333333  0.9279833  0.1441343
## 3  -0.3033333 -0.6227591 -0.1634446  0.3333333 -0.7188603 -1.6259807
## 4  -0.3033333 -0.6227591 -0.1634446  1.3333333  0.7450007  1.9219121
## 5  -0.3033333 -0.6227591 -0.1634446 -1.6666667 -0.1241668 -1.9324941
## 6  -0.3033333 -0.6227591 -0.1634446 -0.6666667 -1.2678081  0.6881956
## 7  -0.3033333 -0.6227591 -0.1634446  0.3333333 -0.4718337 -1.7715745
## 8  -0.3033333 -0.6227591 -0.1634446  1.3333333  1.1933081  0.5272760
## 9  -0.3033333 -0.6227591 -0.1634446  0.3333333  0.8273428 -0.3769385
## 10 -0.3033333 -0.6227591 -0.1634446 -0.6666667  0.2235002  0.9947090
##             x          m          y
## 1   0.3333333  0.9828781 -1.4420726
## 2   0.3333333  0.9279833  0.1441343
## 3   0.3333333 -0.7188603 -1.6259807
## 4   1.3333333  0.7450007  1.9219121
## 5  -1.6666667 -0.1241668 -1.9324941
## 6  -0.6666667 -1.2678081  0.6881956
## 7   0.3333333 -0.4718337 -1.7715745
## 8   1.3333333  1.1933081  0.5272760
## 9   0.3333333  0.8273428 -0.3769385
## 10 -0.6666667  0.2235002  0.9947090

The variables of interest in this example are fwkstrcw (the predictor: female work stress centered within),fwkdiscw (the mediator: female work dissatisfaction centered within),freldiscw (the outcome: female relationship dissatisfaction centered within). 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. (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

Lets look at the descriptives.

#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 plot each leg of the model for a few individuals.

x –> y 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") + 
  geom_smooth(method=lm, se=TRUE, fullrange=TRUE, color="black") + 
  xlab("x = Work Stress (within)") +
  ylab("y = Relationship Dissatisfaction (within)") +
  facet_wrap( ~ id)

x –> m 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") + 
  geom_smooth(method=lm, se=TRUE, fullrange=TRUE, color="blue") + 
  xlab("x = Work Stress (within)") +
  ylab("m = Work Dissatisfaction (within)") +
  facet_wrap( ~ id)

m –> y 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") + 
  geom_smooth(method=lm, se=TRUE, fullrange=TRUE, color="red") + 
  xlab("m = Work Dissatisfaction (within)") +
  ylab("y = Relationship Dissatisfaction (within)") +
  facet_wrap( ~ id)

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.

Restructuring is done in just a few steps. 1. Subset down to the variables of interest.

data1 <- data[,c("id","time","timec","fwkstrcw","fwkdiscw", "freldiscw","x","m","y")]
#look at  data set
head(data1)
##    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
##            y
## 1 -1.4420726
## 2  0.1441343
## 3 -1.6259807
## 4  1.9219121
## 5 -1.9324941
## 6  0.6881956
tail(data1)
##       id time timec    fwkstrcw    fwkdiscw  freldiscw           x
## 2095 200   16     5 -0.95238095 -1.11183723 -0.5148695 -0.95238095
## 2096 200   17     6  0.04761905 -0.06883632  0.8491151  0.04761905
## 2097 200   18     7 -0.95238095  0.01350586 -0.8060573 -0.95238095
## 2098 200   19     8 -1.95238095  1.68779680  0.2437511 -1.95238095
## 2099 200   20     9 -0.95238095 -1.90781162 -0.9516512 -0.95238095
## 2100 200   21    10  1.04761905  0.70883980 -1.5187010  1.04761905
##                m          y
## 2095 -1.11183723 -0.5148695
## 2096 -0.06883632  0.8491151
## 2097  0.01350586 -0.8060573
## 2098  1.68779680  0.2437511
## 2099 -1.90781162 -0.9516512
## 2100  0.70883980 -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"     "fwkstrcw"  "fwkdiscw"  "freldiscw"
## [7] "x"         "m"         "y"
datalong <- melt(data=data1,
               id.vars=c("id","time","timec","fwkstrcw","fwkdiscw","freldiscw","x"),
               na.rm=FALSE, variable.name="dv",
               value.name="z")
#look at updated data set
head(datalong)
##    id time timec   fwkstrcw   fwkdiscw  freldiscw          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    fwkstrcw    fwkdiscw  freldiscw           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  fwkstrcw   fwkdiscw  freldiscw         x dv
## 2101 101    1   -10 0.3333333  0.9828781 -1.4420726 0.3333333  y
## 1    101    1   -10 0.3333333  0.9828781 -1.4420726 0.3333333  m
## 2102 101    2    -9 0.3333333  0.9279833  0.1441343 0.3333333  y
## 2    101    2    -9 0.3333333  0.9279833  0.1441343 0.3333333  m
## 2103 101    3    -8 0.3333333 -0.7188603 -1.6259807 0.3333333  y
## 3    101    3    -8 0.3333333 -0.7188603 -1.6259807 0.3333333  m
##               z dy dm dvnum
## 2101 -1.4420726  1  0     0
## 1     0.9828781  0  1     1
## 2102  0.1441343  1  0     0
## 2     0.9279833  0  1     1
## 2103 -1.6259807  1  0     0
## 3    -0.7188603  0  1     1
tail(datalong)
##       id time timec  fwkstrcw   fwkdiscw  freldiscw         x dv
## 4198 200   19     8 -1.952381  1.6877968  0.2437511 -1.952381  y
## 2098 200   19     8 -1.952381  1.6877968  0.2437511 -1.952381  m
## 4199 200   20     9 -0.952381 -1.9078116 -0.9516512 -0.952381  y
## 2099 200   20     9 -0.952381 -1.9078116 -0.9516512 -0.952381  m
## 4200 200   21    10  1.047619  0.7088398 -1.5187010  1.047619  y
## 2100 200   21    10  1.047619  0.7088398 -1.5187010  1.047619  m
##               z dy dm dvnum
## 4198  0.2437511  1  0     0
## 2098  1.6877968  0  1     1
## 4199 -0.9516512  1  0     0
## 2099 -1.9078116  0  1     1
## 4200 -1.5187010  1  0     0
## 2100  0.7088398  0  1     1

Running Mediation Model -

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

#lme mediation model
model.lme <- lme(fixed = z ~ -1 +
                             dm:fwkstrcw + dm:timec + 
                             dy:fwkdiscw + dy:fwkstrcw + dy:timec,
                  random = ~ -1  +  dm:fwkstrcw + dy:fwkdiscw + dy:fwkstrcw | id, 
                  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)

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

Note: Previously, we had made a coding error (now corrected above) that motivated us to try fitting the model another way. Soon, we will fill more interpretataion in here of the lme() model output (and use this model to obtain confidence intervals).

We also fit the model using the lmer() function in the lme package. (Note, however, that generally lmer() cannot work with heterogeneous variances. We follow the implementation described here (https://stats.idre.ucla.edu/r/faq/how-can-i-perform-mediation-with-multilevel-data-method-2/).

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

#lmer mediation model
mediationmodel.lmer <- lmer(formula = z ~ 0 + 
                                dm + dm:fwkstrcw + dm:timec + 
                                dy + dy:fwkdiscw + dy:fwkstrcw + dy:timec + 
                                (0 + dm:fwkstrcw + dy:fwkdiscw + dy:fwkstrcw | id) + 
                                (0 + dm | timec), 
                  data=datalong) 
#viewing and putting model summary into an object

(medmodelsummary <- summary(mediationmodel.lmer))
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## z ~ 0 + dm + dm:fwkstrcw + dm:timec + dy + dy:fwkdiscw + dy:fwkstrcw +  
##     dy:timec + (0 + dm:fwkstrcw + dy:fwkdiscw + dy:fwkstrcw |  
##     id) + (0 + dm | timec)
##    Data: datalong
## 
## REML criterion at convergence: 12218.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5990 -0.6597 -0.0078  0.6676  3.0724 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr     
##  id       dm:fwkstrcw 0.076191 0.27603           
##           dy:fwkdiscw 0.042141 0.20528  0.55     
##           fwkstrcw:dy 0.006514 0.08071  0.39 0.98
##  timec    dm          0.001014 0.03184           
##  Residual             1.021865 1.01087           
## Number of obs: 4200, groups:  id, 100; timec, 21
## 
## Fixed effects:
##               Estimate Std. Error t value
## dm           2.331e-05  2.313e-02   0.001
## dy          -1.123e-04  2.206e-02  -0.005
## dm:fwkstrcw  1.911e-01  3.562e-02   5.366
## dm:timec    -5.997e-03  3.870e-03  -1.550
## dy:fwkdiscw  1.490e-01  2.885e-02   5.165
## fwkstrcw:dy  1.042e-01  2.412e-02   4.320
## timec:dy    -2.474e-03  3.692e-03  -0.670
## 
## Correlation of Fixed Effects:
##             dm     dy     dm:fwk dm:tmc dy:fwk fwkst:
## dy           0.000                                   
## dm:fwkstrcw  0.000  0.000                            
## dm:timec     0.000  0.000 -0.011                     
## dy:fwkdiscw  0.000  0.000  0.303  0.001              
## fwkstrcw:dy  0.000 -0.001  0.101  0.001  0.139       
## timec:dy     0.000  0.000  0.000  0.002  0.016 -0.014

It converges!

Let’s look more closely at the results.

The fixed effects …

#fixed effects
(FEmatrix <- coef(medmodelsummary))
##                  Estimate  Std. Error      t value
## dm           2.330974e-05 0.023128061  0.001007856
## dy          -1.122918e-04 0.022059146 -0.005090489
## dm:fwkstrcw  1.911433e-01 0.035622537  5.365798326
## dm:timec    -5.996781e-03 0.003869876 -1.549605506
## dy:fwkdiscw  1.490231e-01 0.028850334  5.165385548
## fwkstrcw:dy  1.042102e-01 0.024124453  4.319689812
## timec:dy    -2.474453e-03 0.003692396 -0.670148218

Lets interpret …
2.330973e-05 = intercept in the m as outcome model (a0: We expect this to be zero because we person-centered all the data)
-1.122919e-04 = intercept in the y as outcome model (b0: We expect this to be zero because we person-centered all the data)
1.911433e-01 = effect of x –> m (a: slope of work stressors predicting work dissatisfaction)
-5.996781e-03 = effect of time –> m (time trend in work dissatisfaction)
1.490231e-01 = effect of m –> y (b: slope of work dissatisfaction predicting relationship dissatisfaction)
1.042101e-01 = effect of x –> y (c’: slope of work stressors predicting relationship dissatisfaction, after adjusting for work dissatisfaction)
-2.474453e-03 = effect of time –> y (time trend in relationship dissatisfaction)

Lets put the parameters into named objects, as these will be useful later when calculating confidence intervals.

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

The random effects …

The random effects can be output this way

#random effects
VarCorr(mediationmodel.lmer)
##  Groups   Name        Std.Dev. Corr       
##  id       dm:fwkstrcw 0.276028            
##           dy:fwkdiscw 0.205283 0.551      
##           fwkstrcw:dy 0.080710 0.386 0.982
##  timec    dm          0.031842            
##  Residual             1.010873

But, we found a bit better way to extract …

(REmatrix <- as.data.frame(VarCorr(mediationmodel.lmer)))
##        grp        var1        var2        vcov      sdcor
## 1       id dm:fwkstrcw        <NA> 0.076191327 0.27602776
## 2       id dy:fwkdiscw        <NA> 0.042141076 0.20528292
## 3       id fwkstrcw:dy        <NA> 0.006514040 0.08070960
## 4       id dm:fwkstrcw dy:fwkdiscw 0.031234112 0.55121825
## 5       id dm:fwkstrcw fwkstrcw:dy 0.008598060 0.38594241
## 6       id dy:fwkdiscw fwkstrcw:dy 0.016277623 0.98245570
## 7    timec          dm        <NA> 0.001013944 0.03184249
## 8 Residual        <NA>        <NA> 1.021864597 1.01087319

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

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

sig2_a <- REmatrix[1,4]

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

sig2_b <- REmatrix[2,4]

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

sig2_cprime <- REmatrix[3,4]

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

The residual variance of the mediator variable, \(\sigma^2_{em}\), is calculated as \((1.0219 + 0.0010)\) = 1.0228785.

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, \(\sigma_{ajbj}\) is 0.0312341

covajbj <- REmatrix[4,4]

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

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

ie = ab + \(\sigma_{ajbj}\)
= (1.911433e-01*1.490231e-01) + 0.031234118
= 0.05971889

indirecteffect <- a*b + covajbj

Calculating total effect c

c = c’ + ab + \(\sigma_{ajbj}\)
= 1.042101e-01 + (1.911433e-01*1.490231e-01) + 0.031234118
= 0.163929

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.05971889 / 0.163929    
  = 0.3642973    
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.031234118 / 0.05971889    
= 0.5230191  
percentcovariance <- 100*(covajbj/indirecteffect)

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

Obtaining Confidence Intervals for the Indirect Effect in 1-1-1 Mediation Models

We also want to obtain confidence intervals for the indirect effects, ie = ab + \(\sigma_{ajbj}\) = 0.05971889

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 or bootstrapping approach. For 1-1-1 Mediation models, Bauer et al. (2006) adapted the approach used by MacKinnon, Lockwood, and Williams (2004), and coded this into an online calculator (which produces and runs R code).

As described at http://quantpsy.org/medmc/medmc111.htm … “This approach can be used as long as seven pieces of information are available from the results of a multilevel Mediation model. These estimates can be found using most multilevel modeling software. The estimates will include fixed and random effects as well as estimates from the asymptotic covariance matrix (ACM), or the covariance matrix of the model parameters.”

The information is used to simulate repeated sampling of the indirect effects, the many estimates of which are used to compute confidence intervals for the observed indirect effect.

We have extracted many of the needed pieces of information already. We now obtain the additional ones, and then run the R code that does the simulations (and would be produced by the on-line calculator). (Of note, the notation used on the website is a bit different than what is used here.)

The piece we are missing is the asymptotic covariance matrix (ACM). For our lmer() model, this can be obtained using vcov(object) (see http://quantpsy.org/interact/acov.htm)

#Getting asymptotic covariance matrix for the fixed effects
(ACM_FE <- vcov(mediationmodel.lmer))
## 7 x 7 Matrix of class "dpoMatrix"
##                        dm            dy   dm:fwkstrcw      dm:timec
## dm           5.349072e-04  1.127835e-09 -4.043144e-07  1.464154e-09
## dy           1.127835e-09  4.866059e-04  3.171200e-09 -5.387093e-11
## dm:fwkstrcw -4.043144e-07  3.171200e-09  1.268965e-03 -1.497342e-06
## dm:timec     1.464154e-09 -5.387093e-11 -1.497342e-06  1.497594e-05
## dy:fwkdiscw -1.753144e-08  3.280616e-08  3.114940e-04  1.411046e-07
## fwkstrcw:dy -1.655596e-08 -5.059067e-07  8.686131e-05  8.319960e-08
## timec:dy     4.751544e-10  1.833909e-09  1.999805e-08  2.697543e-08
##               dy:fwkdiscw   fwkstrcw:dy      timec:dy
## dm          -1.753144e-08 -1.655596e-08  4.751544e-10
## dy           3.280616e-08 -5.059067e-07  1.833909e-09
## dm:fwkstrcw  3.114940e-04  8.686131e-05  1.999805e-08
## dm:timec     1.411046e-07  8.319960e-08  2.697543e-08
## dy:fwkdiscw  8.323417e-04  9.704167e-05  1.745953e-06
## fwkstrcw:dy  9.704167e-05  5.819892e-04 -1.279899e-06
## timec:dy     1.745953e-06 -1.279899e-06  1.363379e-05
#Pulling the asymptotic covariances for the fixed effects of interest
vara <- ACM_FE[3,3]
varb <- ACM_FE[5,5]
covab <- ACM_FE[3,5]

#Pulling the asymptotic covariances for the random effects
#Using the vcov_vc() in the bootmlm library
(ACM_RE <- vcov_vc(mediationmodel.lmer, sd_cor = FALSE, print_names = TRUE))
##                                var_dm:fwkstrcw|id
## var_dm:fwkstrcw|id                   1.964042e-04
## cov_dy:fwkdiscw.dm:fwkstrcw|id      -3.210168e-06
## cov_fwkstrcw:dy.dm:fwkstrcw|id      -3.376113e-06
## var_dy:fwkdiscw|id                  -1.507729e-05
## cov_fwkstrcw:dy.dy:fwkdiscw|id      -5.547562e-06
## var_fwkstrcw:dy|id                  -2.404047e-06
## var_dm|timec                        -8.677451e-07
## sigma2                              -1.861085e-05
##                                cov_dy:fwkdiscw.dm:fwkstrcw|id
## var_dm:fwkstrcw|id                              -3.210168e-06
## cov_dy:fwkdiscw.dm:fwkstrcw|id                   1.285140e-06
## cov_fwkstrcw:dy.dm:fwkstrcw|id                  -2.199432e-06
## var_dy:fwkdiscw|id                              -1.109432e-05
## cov_fwkstrcw:dy.dy:fwkdiscw|id                  -4.824872e-06
## var_fwkstrcw:dy|id                              -1.855785e-06
## var_dm|timec                                     4.680061e-08
## sigma2                                           2.474099e-06
##                                cov_fwkstrcw:dy.dm:fwkstrcw|id
## var_dm:fwkstrcw|id                              -3.376113e-06
## cov_dy:fwkdiscw.dm:fwkstrcw|id                  -2.199432e-06
## cov_fwkstrcw:dy.dm:fwkstrcw|id                   1.660580e-06
## var_dy:fwkdiscw|id                              -3.770488e-06
## cov_fwkstrcw:dy.dy:fwkdiscw|id                  -1.290609e-06
## var_fwkstrcw:dy|id                              -6.403165e-07
## var_dm|timec                                     4.148771e-08
## sigma2                                           8.095173e-07
##                                var_dy:fwkdiscw|id
## var_dm:fwkstrcw|id                  -1.507729e-05
## cov_dy:fwkdiscw.dm:fwkstrcw|id      -1.109432e-05
## cov_fwkstrcw:dy.dm:fwkstrcw|id      -3.770488e-06
## var_dy:fwkdiscw|id                   2.576054e-06
## cov_fwkstrcw:dy.dy:fwkdiscw|id      -1.463768e-06
## var_fwkstrcw:dy|id                  -1.216005e-06
## var_dm|timec                        -1.026542e-07
## sigma2                              -2.091520e-07
##                                cov_fwkstrcw:dy.dy:fwkdiscw|id
## var_dm:fwkstrcw|id                              -5.547562e-06
## cov_dy:fwkdiscw.dm:fwkstrcw|id                  -4.824872e-06
## cov_fwkstrcw:dy.dm:fwkstrcw|id                  -1.290609e-06
## var_dy:fwkdiscw|id                              -1.463768e-06
## cov_fwkstrcw:dy.dy:fwkdiscw|id                   8.023921e-09
## var_fwkstrcw:dy|id                              -2.710824e-07
## var_dm|timec                                    -4.270987e-08
## sigma2                                           1.161193e-07
##                                var_fwkstrcw:dy|id  var_dm|timec
## var_dm:fwkstrcw|id                  -2.404047e-06 -8.677451e-07
## cov_dy:fwkdiscw.dm:fwkstrcw|id      -1.855785e-06  4.680061e-08
## cov_fwkstrcw:dy.dm:fwkstrcw|id      -6.403165e-07  4.148771e-08
## var_dy:fwkdiscw|id                  -1.216005e-06 -1.026542e-07
## cov_fwkstrcw:dy.dy:fwkdiscw|id      -2.710824e-07 -4.270987e-08
## var_fwkstrcw:dy|id                   8.839147e-08 -1.440627e-08
## var_dm|timec                        -1.440627e-08  1.406784e-05
## sigma2                               1.589347e-07 -5.128078e-06
##                                       sigma2
## var_dm:fwkstrcw|id             -1.861085e-05
## cov_dy:fwkdiscw.dm:fwkstrcw|id  2.474099e-06
## cov_fwkstrcw:dy.dm:fwkstrcw|id  8.095173e-07
## var_dy:fwkdiscw|id             -2.091520e-07
## cov_fwkstrcw:dy.dy:fwkdiscw|id  1.161193e-07
## var_fwkstrcw:dy|id              1.589347e-07
## var_dm|timec                   -5.128078e-06
## sigma2                          5.213382e-04
#Pulling the asymptotic covariances for the fixed effects of interest
varcovajbj=ACM_RE[2,2]

Took a bit if searching to find a function that extracted the ACM for the random effects, but it works! See https://rdrr.io/github/marklhc/bootmlm/man/vcov_vc.html (Thanks to Chang Liu for locating the solution!)

Once we have everything we need (and named correctly), we use Preacher et al’s online calculator code to obtain the confidence intervals. Preacher, K. J., & Selig, J. P. (2010, July). Monte Carlo method for assessing multilevel Mediation: An interactive tool for creating confidence intervals for indirect effects in 1-1-1 multilevel models [Computer software]. Available from http://quantpsy.org/.

################################################
# This code was generated by the online        #
# calculator. No need to generate again.       #
################################################
a=a #0.1911433
b=b #0.1490231
covajbj=covajbj #0.03123412
vara=vara #0.001268965
varb=varb #0.0008323414
covab=covab #0.0003114941
varcovajbj=varcovajbj #0.0001
rep=20000
conf=95
dvec=rnorm(rep)
avec=dvec*sqrt(vara)+a
bvec=dvec*covab/sqrt(vara)+sqrt(varb)*rnorm(rep,sd=sqrt(1-(covab^2)/(vara*varb)))+b
cvec=rnorm(rep)*sqrt(varcovajbj)+covajbj
ab=avec*bvec+cvec
low=(1-conf/100)/2
upp=((1-conf/100)/2)+(conf/100)
LL=quantile(ab,low)
UL=quantile(ab,upp)
LL4=format(LL,digits=4)
UL4=format(UL,digits=4)
################################################
# The number of columns in the histogram can   #
# be changed by replacing 'FD' below with      #
# an integer value.                            #
################################################
hist(ab,breaks='FD',col='skyblue',xlab=paste(conf,'% Confidence Interval ','LL',LL4,'  UL',UL4),
main='Distribution of Indirect Effect')

So, we now have the following. The indirect effect is 0.0597189, with 95% CI of [0.0445,0.07962]

Yay!

3. An Example using Multi-Step Estimation and Bootstrapping

There is another way that people are estimating the multilevel mediation model. In this approach, the two models are run separately,

\[M_{jt} = a0_{j} + a_{j}X_{jt} + eM_{jt}\] and \[Y_{jt} = b0_{j} + b_{j}M_{jt} + c'_{i}X_{jt} + eY_{jt}\] and run many times, bootstrapping. The parameters are summarized to get a point estimate and confidence intervals.

A convenience function indirectMLM.R was written by Elizabeth Page-Gould. It can be downloaded along with a worked through example from her website here, http://www.page-gould.com/r/indirectmlm/

We use that code to run the example from Bolger & Laurenceau (2013) Chapter 9 (and ceck consistency with above appraoch).

This approach makes use of the boot() function in the boot library, and the source files with the indirectMLM.R functions.

#### Load Packages ####
library(boot)
## 
## Attaching package: 'boot'
## The following object is masked from 'package:psych':
## 
##     logit
source("indirectMLM.R")

This approach makes use of the original single-entry data. We already have these data in place.

#### Read Data ####
data1 <- data1

head(data1)
##    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
##            y
## 1 -1.4420726
## 2  0.1441343
## 3 -1.6259807
## 4  1.9219121
## 5 -1.9324941
## 6  0.6881956

Setting up the bootstrapped model. Instructions come from the README file in the download.

indirect.mlm

Bootstrapping Accurate Indirect Effects in Multilevel Models

DESCRIPTION

This function can be used to bootstrap accurate indirect effects in multilevel mediation models where both paths that comprise the indirect effect are random. This is just a function for the boot() function from the package, “boot” (Canty & Ripley, 2015). The boot function requires the user to pass a function to it using the “statistic” argument. When you set the statistic argument to equal our function, indirect.mlm, then the boot function will conduct a multilevel mediation.

USAGE

boot( data, statistic=indirect.mlm, R, strata, y, x, mediator, group.id, covariates=NULL, random.a=T, random.b=T, random.c=T, uncentered.x=T, between.x=F, between.m=F )

indirect.mlm.summary( boot.object )

ARGUMENTS …

Applying to our data

#### Analyses ####
# Note: This will probably take a while to run! Be patient!
mediated.mlm <- boot( data=data1, statistic=indirect.mlm, R=100, strata=data1$id, 
                      y="y", x="x", mediator="m", group.id="id",
                      between.m=F, uncentered.x=F)

and then there is a summary function indirect.mlm.summary() to look at the output in the boot object.

#### Print Output ####
indirect.mlm.summary( mediated.mlm )
## #### Population Covariance ####
## Covariance of Random Slopes a and b: 0.01 [0.004, 0.026]
## 
## 
## #### Indirect Effects ####
## # Within-subject Effects
## Unbiased Estimate of Within-subjects Indirect Effect: 0.038 [0.029, 0.06]
## Biased Estimate of Within-subjects Indirect Effect: 0.029 [0.018, 0.04]
## Bias in Within-subjects Indirect Effect: 0.01 [0.004, 0.025]
## 
## 
## #### Total Effect ####
## Unbiased Estimate of Total Effect: 0.154 [0.114, 0.196]
## Biased Total Effect of X on Y (c path): 0.158 [0.115, 0.197]
## Bias in Total Effect: 0.004 [0, 0.009]
## 
## 
## #### Direct Effects ####
## Direct Effect of Predictor on Dependent Variable (c' path): 0.116 [0.067, 0.156]
## Within-subjects Effect of Predictor on Mediator (a path for group-mean centered predictor): 0.187 [0.142, 0.254]
## Within-subjects Effect of Mediator on Dependent Variable (b path for group-mean centered mediator): 0.152 [0.112, 0.193]

These are pretty close to what we got above, although the covariance between \(a_{j}\) and \(b_{j}\), \(sigma_{ajbj}\), is a bit lower than what we got above.

One of the convenient things about the function, is that it can do the person centering, and thus also provides information about the between-person associations as well.

#not run
mediated.mlm <- boot(data=data, statistic=indirect.mlm, R=100, strata=data$id, 
                      y="freldis", x="fwkstr", mediator="fwkdis", group.id="id",
                      between.m=F, uncentered.x=F )

4. Conclusion

Ok - we have attempted to push into the within-person (1-1-1) mediation models. I do not yet have much experience with these models, but see their value for examining within-person processes. There are still a variety of details that must be worked out and understood as we try to fit these models in R and publish the analyses. For example, we still have to figure out how to display the raw data and the results in ways that facilitate intuitive understanding.

For further information and discussion of moderated mediation in multilevel models check out the paper here http://quantpsy.org/pubs/bauer_preacher_gil_2006.pdf.

Good Luck! And, please be careful.