Mediation — practical

Path mediation then latent SEM mediation with bootstrap BCAs

Author

Affiliation

George Michaelides

Norwich Business School, University of East Anglia

Published

July 7, 2026

Learning objectives

By the end of this tutorial you will be able to:

label regression paths in lavaan with * and define the indirect effect with :=
fit a path mediation model with bootstrap standard errors and extract BCa confidence intervals
extend the mediation model to a full latent-variable SEM with three factors
report mediation results in APA format, including the indirect effect and its confidence interval

Setup

Setup code (you have seen this before)

library(haven)
library(lavaan)
library(tidySEM)

The mediation model

We test whether positive affect (PA) mediates the relationship between job control and creativity:

\(X\) = job control (con)
\(M\) = positive affect (pa)
\(Y\) = creativity (creat)

We fit this model in two forms:

Path mediation — manifest variables (mean scores), as a direct extension of Tutorial 09
Latent SEM mediation — latent variables with item-level indicators, as an extension of Tutorial 10

Both use the same := syntax to define the indirect effect.

Worked example 1: Path mediation with manifest variables

Computing the manifest variables

ps$con   <- (ps$con1 + ps$con3 + ps$con6 + ps$con8 + ps$con9 + ps$con10) / 6
ps$pa    <- (ps$pa1 + ps$pa2 + ps$pa3) / 3
ps$creat <- (ps$creat1 + ps$creat2 + ps$creat3) / 3

Specifying the model with labelled paths

The key innovation is labelling the regression paths with * so we can refer to them in := definitions:

model_med_path <- '
  creat ~ b * pa + c * con
  pa    ~ a * con

  indirect := a * b
  direct   := c
  total    := c + (a * b)
'

Reading the syntax:

creat ~ b * pa + c * con — the path from pa to creat is labelled b; the path from con to creat is labelled c
pa ~ a * con — the path from con to pa is labelled a
indirect := a * b — defines the indirect effect as the product \(ab\)
direct := c — defines the direct effect as \(c'\)
total := c + (a * b) — defines the total effect as \(c' + ab\)

Fitting with bootstrap standard errors

To test the indirect effect, we use bootstrap confidence intervals. This requires two changes:

Switch the estimator from MLR to ML (MLR is incompatible with se = "boot")
Request bootstrap standard errors with se = "boot" and bootstrap = 1000

set.seed(1234)
fit_med_path <- sem(model_med_path, data = ps,
                    estimator = "ML", se = "boot", bootstrap = 1000)

Estimator constraint with bootstrapping

When you request se = "boot", lavaan requires estimator = "ML" — the robust MLR estimator cannot be combined with bootstrapping. This is because the bootstrap builds its own empirical sampling distribution, making the Satorra–Bentler scaling correction unnecessary. The point estimates are the same under ML and MLR; only the standard errors differ, and for the indirect effect we use the bootstrap standard errors.

Reading the output

summary(fit_med_path, standardized = TRUE, rsquare = TRUE)

## lavaan 0.6-21 ended normally after 1 iteration
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                         5
## 
##                                                   Used       Total
##   Number of observations                           332         343
## 
## Model Test User Model:
##                                                       
##   Test statistic                                 0.000
##   Degrees of freedom                                 0
## 
## Parameter Estimates:
## 
##   Standard errors                            Bootstrap
##   Number of requested bootstrap draws             1000
##   Number of successful bootstrap draws            1000
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   creat ~                                                               
##     pa         (b)    0.185    0.042    4.399    0.000    0.185    0.235
##     con        (c)   -0.091    0.069   -1.324    0.186   -0.091   -0.073
##   pa ~                                                                  
##     con        (a)    0.229    0.081    2.816    0.005    0.229    0.144
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .creat             0.690    0.046   15.162    0.000    0.690    0.945
##    .pa                1.150    0.072   15.970    0.000    1.150    0.979
## 
## R-Square:
##                    Estimate
##     creat             0.055
##     pa                0.021
## 
## Defined Parameters:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##     indirect          0.042    0.018    2.294    0.022    0.042    0.034
##     direct           -0.091    0.069   -1.324    0.186   -0.091   -0.073
##     total            -0.049    0.073   -0.668    0.504   -0.049   -0.039

The output now includes a Defined Parameters section at the bottom, showing the estimates for indirect, direct, and total. These are the derived parameters we defined with :=.

Extracting bootstrap confidence intervals

For clean BCa confidence intervals, use parameterEstimates() with boot.ci.type = "bca.simple":

pe_med <- parameterEstimates(fit_med_path, standardized = TRUE,
                             boot.ci.type = "bca.simple")
pe_med[pe_med$lhs %in% c("indirect", "direct", "total"),
       c("lhs", "est", "se", "z", "pvalue", "ci.lower", "ci.upper")]

##        lhs    est    se      z pvalue ci.lower ci.upper
## 7 indirect  0.042 0.018  2.294  0.022    0.012    0.085
## 8   direct -0.091 0.069 -1.324  0.186   -0.226    0.043
## 9    total -0.049 0.073 -0.668  0.504   -0.195    0.096

The columns:

est — the estimate of the defined parameter
se — the bootstrap standard error
z, pvalue — the Wald statistic and p-value (based on the bootstrap SE)
ci.lower, ci.upper — the BCa 95% confidence interval

The key row is indirect. If the BCa 95% CI for the indirect effect does not include zero, the indirect effect is significant at the .05 level — we have evidence of mediation.

Interpreting the results

Read the three defined parameters:

indirect (\(ab\)): the mediated effect of con on creat through pa. If the CI excludes zero, mediation is present.
direct (\(c'\)): the direct effect of con on creat, controlling for pa. If this is not significant, mediation is full; if it is significant, mediation is partial.
total (\(c' + ab\)): the total effect of con on creat. This may or may not be significant — a non-significant total effect does not rule out mediation (see Tutorial 11).

The total effect can be non-significant even when mediation is present

If the indirect effect is positive and the direct effect is negative (or vice versa), they can cancel each other out, producing a non-significant total effect. This is called inconsistent mediation or a suppressor effect. Always examine the indirect effect directly, not just the total effect.

Worked example 2: Latent SEM mediation

Now we extend the mediation model to use latent variables for all three constructs. This is the full SEM version, combining the measurement model from Tutorial 10 with the mediation structure from Tutorial 11.

The model

# Build parcels for job control (6 items -> 3 parcels)
ps$p1 <- (ps$con9 + ps$con10) / 2
ps$p2 <- (ps$con1 + ps$con8)  / 2
ps$p3 <- (ps$con3 + ps$con6)  / 2

model_med_latent <- '
  # Measurement model
  con   =~ p1 + p2 + p3
  pa    =~ pa1 + pa2 + pa3
  creat =~ creat1 + creat2 + creat3

  # Structural model with labelled paths
  creat ~ b * pa + c * con
  pa    ~ a * con

  # Defined parameters
  indirect := a * b
  direct   := c
  total    := c + (a * b)
'

Fitting with bootstrap

set.seed(1234)
fit_med_latent <- sem(model_med_latent, data = ps,
                      estimator = "ML", se = "boot", bootstrap = 1000)

Reading the output

summary(fit_med_latent, standardized = TRUE, fit.measures = TRUE, rsquare = TRUE)

## lavaan 0.6-21 ended normally after 28 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of model parameters                        21
## 
##                                                   Used       Total
##   Number of observations                           332         343
## 
## Model Test User Model:
##                                                       
##   Test statistic                                25.528
##   Degrees of freedom                                24
##   P-value (Chi-square)                           0.378
## 
## Model Test Baseline Model:
## 
##   Test statistic                              1289.618
##   Degrees of freedom                                36
##   P-value                                        0.000
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    0.999
##   Tucker-Lewis Index (TLI)                       0.998
## 
## Loglikelihood and Information Criteria:
## 
##   Loglikelihood user model (H0)              -3551.266
##   Loglikelihood unrestricted model (H1)      -3538.503
##                                                       
##   Akaike (AIC)                                7144.533
##   Bayesian (BIC)                              7224.441
##   Sample-size adjusted Bayesian (SABIC)       7157.828
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.014
##   90 Percent confidence interval - lower         0.000
##   90 Percent confidence interval - upper         0.047
##   P-value H_0: RMSEA <= 0.050                    0.965
##   P-value H_0: RMSEA >= 0.080                    0.000
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.028
## 
## Parameter Estimates:
## 
##   Standard errors                            Bootstrap
##   Number of requested bootstrap draws             1000
##   Number of successful bootstrap draws            1000
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   con =~                                                                
##     p1                1.000                               0.687    0.772
##     p2                0.810    0.070   11.616    0.000    0.557    0.797
##     p3                0.898    0.073   12.330    0.000    0.617    0.794
##   pa =~                                                                 
##     pa1               1.000                               0.956    0.806
##     pa2               1.046    0.079   13.300    0.000    1.000    0.780
##     pa3               1.090    0.078   14.017    0.000    1.042    0.843
##   creat =~                                                              
##     creat1            1.000                               0.745    0.753
##     creat2            1.047    0.079   13.193    0.000    0.780    0.818
##     creat3            1.127    0.083   13.513    0.000    0.840    0.855
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   creat ~                                                               
##     pa         (b)    0.205    0.053    3.896    0.000    0.264    0.264
##     con        (c)   -0.102    0.070   -1.442    0.149   -0.094   -0.094
##   pa ~                                                                  
##     con        (a)    0.227    0.090    2.530    0.011    0.163    0.163
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .p1                0.320    0.039    8.293    0.000    0.320    0.404
##    .p2                0.178    0.021    8.358    0.000    0.178    0.364
##    .p3                0.224    0.030    7.431    0.000    0.224    0.370
##    .pa1               0.492    0.075    6.565    0.000    0.492    0.350
##    .pa2               0.645    0.073    8.892    0.000    0.645    0.392
##    .pa3               0.444    0.067    6.653    0.000    0.444    0.290
##    .creat1            0.423    0.042   10.059    0.000    0.423    0.432
##    .creat2            0.300    0.038    7.825    0.000    0.300    0.330
##    .creat3            0.259    0.040    6.488    0.000    0.259    0.268
##     con               0.472    0.062    7.662    0.000    1.000    1.000
##    .pa                0.890    0.103    8.607    0.000    0.973    0.973
##    .creat             0.517    0.063    8.148    0.000    0.930    0.930
## 
## R-Square:
##                    Estimate
##     p1                0.596
##     p2                0.636
##     p3                0.630
##     pa1               0.650
##     pa2               0.608
##     pa3               0.710
##     creat1            0.568
##     creat2            0.670
##     creat3            0.732
##     pa                0.027
##     creat             0.070
## 
## Defined Parameters:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##     indirect          0.047    0.022    2.082    0.037    0.043    0.043
##     direct           -0.102    0.070   -1.442    0.149   -0.094   -0.094
##     total            -0.055    0.075   -0.734    0.463   -0.051   -0.051

The output now has all the sections you met in Tutorial 10 (Latent variables, Regressions, Variances, R-Square) plus the Defined Parameters section at the bottom.

Bootstrap confidence intervals

pe_med_latent <- parameterEstimates(fit_med_latent, standardized = TRUE,
                                    boot.ci.type = "bca.simple")
pe_med_latent[pe_med_latent$lhs %in% c("indirect", "direct", "total"),
              c("lhs", "est", "se", "z", "pvalue", "ci.lower", "ci.upper",
                "std.all")]

##         lhs    est    se      z pvalue ci.lower ci.upper std.all
## 25 indirect  0.047 0.022  2.082  0.037    0.011    0.098   0.043
## 26   direct -0.102 0.070 -1.442  0.149   -0.241    0.041  -0.094
## 27    total -0.055 0.075 -0.734  0.463   -0.208    0.091  -0.051

Compare the indirect effect from the latent SEM to the indirect effect from the path mediation. The latent SEM estimate should be larger in absolute value — disattenuation removes the measurement error from all three constructs, giving an unbiased estimate of the mediated effect.

Fit indices

The latent SEM is overidentified (unlike the path mediation, which was just-identified), so we can examine fit:

fitMeasures(fit_med_latent, c("chisq", "df", "pvalue",
                              "cfi", "tli",
                              "rmsea", "rmsea.ci.lower", "rmsea.ci.upper",
                              "srmr"))

##          chisq             df         pvalue            cfi            tli 
##         25.528         24.000          0.378          0.999          0.998 
##          rmsea rmsea.ci.lower rmsea.ci.upper           srmr 
##          0.014          0.000          0.047          0.028

If the fit indices meet the Hu and Bentler cutoffs, we have evidence that the model — both its measurement and structural parts — is consistent with the data.

Drawing the mediation model

tidySEM::graph_sem(fit_med_latent, std = TRUE)

Path diagram for the latent-variable mediation model. Job control (con) is a latent factor measured by three parcels; positive affect (pa) and creativity (creat) are latent factors measured by three items each. The indirect effect of con on creat through pa is a * b.

Writing up the results

An APA-style write-up of the latent SEM mediation:

Structural equation modelling was used to test whether positive affect (PA) mediated the relationship between job control and creativity. Job control was modelled as a latent factor measured by three parcels (p1, p2, p3); PA and creativity were modelled as latent factors measured by their three respective items. The structural paths were estimated with maximum likelihood and 1,000 bootstrap replications; bias-corrected and accelerated (BCa) confidence intervals are reported for the indirect effect.

The model fit the data well, \(\chi^2(24) = ...\), CFI = …, TLI = …, RMSEA = …, 90% CI […, …], SRMR = …. The indirect effect of job control on creativity through PA was \(ab = ...\), 95% BCa CI […, …], indicating significant mediation. The direct effect of job control on creativity was \(c' = ...\), \(p = ...\), suggesting [full/partial] mediation. The total effect was \(c = ...\), \(p = ...\). Positive affect thus [fully/partially] mediated the relationship between job control and creativity.

(Replace the ... with the actual values from the output above.)

Reporting bootstrap CIs

When you report a bootstrap CI, state:

The number of bootstrap replications (e.g. 1,000)
The type of CI (BCa is preferred)
The estimate and its CI
Whether the CI includes zero

Do not report a p-value for the indirect effect as your primary test — the bootstrap CI is the test. The p-value in the parameterEstimates() output is based on the bootstrap SE and is less accurate than the CI itself.

Common errors

lavaan WARNING: bootstrap only supported with estimator = ML

You requested se = "boot" but left estimator = "MLR" (or did not specify). Change to estimator = "ML". The point estimates are identical; only the standard errors change.

`set.seed()` not called

Without set.seed(), the bootstrap results will differ every time you run the model. Always set a seed before fitting a bootstrap model so your results are reproducible. We use set.seed(1234) throughout this series.

Very wide bootstrap CIs

If the bootstrap CI for the indirect effect is very wide, it may indicate a small sample, a weak indirect effect, or high variability in the bootstrap samples. Increasing the number of bootstrap replications (e.g. to 5,000) will narrow the CI slightly but will not fix a fundamentally weak effect.

Negative variance (Heywood case) in the latent SEM

If the latent SEM produces a Heywood case — a negative variance for one of the latent factors or indicators — check the measurement model’s fit on its own (fit a CFA first). A poorly-fitting measurement model can distort the structural estimates.

Exercises

Fit the path mediation model from the tutorial (manifest variables, ML with se = "boot", bootstrap = 1000). Set set.seed(1234). Extract the BCa confidence interval for the indirect effect. Is the indirect effect significant? Is the mediation full or partial?
Fit the latent SEM mediation model from the tutorial. Compare the indirect effect estimate to the path mediation estimate from Exercise 1. Is the latent SEM estimate larger? Why?
Add a second mediator to the path mediation model: negative affect (neg). The model becomes: creat ~ pa + neg + con, pa ~ con, neg ~ con. Define two specific indirect effects (indirect_pa := a1 * b1 and indirect_neg := a2 * b2) and the total indirect effect. Fit with bootstrap. Which mediator shows a significant indirect effect?
Using the latent SEM mediation model from Exercise 2, write a complete APA-style results paragraph. Report the model fit, the indirect effect with its BCa CI, the direct effect, and the total effect. State whether the mediation is full or partial.

Solutions

Exercise 1

set.seed(1234)
fit_ex1 <- sem(model_med_path, data = ps,
               estimator = "ML", se = "boot", bootstrap = 1000)
pe1 <- parameterEstimates(fit_ex1, boot.ci.type = "bca.simple")
pe1[pe1$lhs %in% c("indirect", "direct", "total"),
    c("lhs", "est", "ci.lower", "ci.upper")]

##        lhs    est ci.lower ci.upper
## 7 indirect  0.042    0.012    0.085
## 8   direct -0.091   -0.226    0.043
## 9    total -0.049   -0.195    0.096

Check whether the ci.lower and ci.upper for the indirect row exclude zero. If they do, the indirect effect is significant. Then check whether the direct row’s CI includes zero — if it does, the mediation is full; if not, it is partial.

Exercise 2

set.seed(1234)
fit_ex2 <- sem(model_med_latent, data = ps,
               estimator = "ML", se = "boot", bootstrap = 1000)
pe2 <- parameterEstimates(fit_ex2, boot.ci.type = "bca.simple")
pe2[pe2$lhs %in% c("indirect", "direct", "total"),
    c("lhs", "est", "ci.lower", "ci.upper")]

##         lhs    est ci.lower ci.upper
## 25 indirect  0.047    0.011    0.098
## 26   direct -0.102   -0.241    0.041
## 27    total -0.055   -0.208    0.091

The latent SEM indirect effect should be larger in absolute value than the path mediation indirect effect. This is disattenuation — the latent variables remove measurement error, producing an unbiased estimate of the mediated effect.

Exercise 3

ps$neg <- (ps$na1 + ps$na2 + ps$na3) / 3

model_2med <- '
  creat ~ b1 * pa + b2 * neg + c * con
  pa    ~ a1 * con
  neg   ~ a2 * con

  indirect_pa  := a1 * b1
  indirect_neg := a2 * b2
  total_indirect := (a1 * b1) + (a2 * b2)
  direct        := c
  total         := c + (a1 * b1) + (a2 * b2)
'

set.seed(1234)
fit_2med <- sem(model_2med, data = ps,
                estimator = "ML", se = "boot", bootstrap = 1000)
pe_2med <- parameterEstimates(fit_2med, boot.ci.type = "bca.simple")
pe_2med[pe_2med$lhs %in% c("indirect_pa", "indirect_neg",
                           "total_indirect", "direct", "total"),
        c("lhs", "est", "ci.lower", "ci.upper")]

##               lhs    est ci.lower ci.upper
## 10    indirect_pa  0.060    0.022    0.123
## 11   indirect_neg -0.064   -0.137   -0.018
## 12 total_indirect -0.004   -0.071    0.056
## 13         direct -0.055   -0.186    0.077
## 14          total -0.060   -0.204    0.079

Check the CIs for indirect_pa and indirect_neg. The mediator whose CI excludes zero shows a significant specific indirect effect.

Exercise 4

Using the output from Exercise 2, write the paragraph following the template in the “Writing up the results” section above. Replace each ... with the actual value from the parameterEstimates() and fitMeasures() output. State whether the indirect effect’s BCa CI excludes zero (significant mediation) and whether the direct effect is significant (full vs. partial mediation).

Summary & next tutorial

You can now fit mediation models in lavaan — both path mediation with manifest variables and full latent-variable SEM mediation — using labelled paths and the := operator to define the indirect effect. You can extract BCa bootstrap confidence intervals and report the results in APA format. This completes the core lavaan toolkit: you can specify, estimate, and report CFA, path analysis, SEM, and mediation models. In the next tutorial we turn to the question of reporting — how to export lavaan results to publication-quality tables and figures, with a focus on the texreg package.

References

Hayes, A. F. (2018). Introduction to mediation, moderation, and conditional process analysis: A regression-based approach.

Kelloway, E. K. (2015). Using mplus for structural equation modeling (2nd ed.). SAGE Publications.

Rosseel, Y. (2012). Lavaan: An r package for structural equation modeling. Journal of Statistical Software, 48(2), 1–36. https://doi.org/10.18637/jss.v048.i02

Sobel, M. E. (1990). Effect analysis and causation in linear structural equation models. Sociological Methods & Research, 19(1), 45–71.

Learning objectives

Setup

The mediation model

Worked example 1: Path mediation with manifest variables

Computing the manifest variables

Specifying the model with labelled paths

Fitting with bootstrap standard errors

Reading the output

Extracting bootstrap confidence intervals

Interpreting the results

Worked example 2: Latent SEM mediation

The model

Fitting with bootstrap

Reading the output

Bootstrap confidence intervals

Fit indices

Drawing the mediation model

Writing up the results

Common errors

set.seed() not called

Very wide bootstrap CIs

Negative variance (Heywood case) in the latent SEM

Exercises

Summary & next tutorial

Further reading

References

`set.seed()` not called