I have just created my first mediation model using sem() with the lavaan package in R. I am using a bootstrapping with 5000 resamples and BCA to calculate the confidence intervals at a 0.9 level.

Now, while researching the web for the right procedure, some practioners added covariances between the mediators into the SEM model, others do not (e.g. Example Covariances included Example Covariances not included). Literature on this topic is (at least for me) very confusing.

When exactly is it needed to include these covariances of the mediators and does it harm the model to not include them?

I tried including them and my CIs got more robust, however my model fit is not given anymore as I end up with negative degrees of freedom when adding every inbetween mediator interaction.

My model includes 4 mediators: v_13, v_36, v_3, and v_5

```
mediation_model <- '
#Regressions Meditation BC -> Aff -> Exp
v_13 ~ n1 * BC
v_36 ~ l2 * v_13 + l1 * BC + ExpMan
#Regressions Meditation Exp -> Disco -> Rating
v_3 ~ d3 * v_36 + d1 * BC + ExpMan + v_30 + PerfMan + d2 * v_13
v_6 ~ b3 * v_36 + b4 * v_3 + b1 * BC + ExpMan + v_30 + PerfMan + b2 * v_13 + b5 * v_5
#Regressions Mediation Exp -> Disco -> Sat
v_5 ~ g3 * v_36 + g4 * v_3 + g1 * BC + ExpMan + v_30 + PerfMan + g2 * v_13
#indirect effects
ACME1 := n1 * l2
ACME2 := d3 * b4
ACME3 := d3 * g4
indESP := g3 * b5
indEDSP := d3 * g4 * b5
indBSP := g1*b5
indBDP := d1*b4
indBDSP := d1*g4*b5
indBEP := l1*b3
indBEDP := l1*d3*b4
indBEDSP := l1*d3*g4*b5
indBESP := l1*g3*b5
indBAP := n1*b2
indBASP := n1*g2*b5
indBADP := n1*d2*b4
indBADSP := n1*d2*g4*b5
indBAEP := n1*l2*b3
indBAESP := n1*l2*g3*b5
indBAEDP := n1*l2*d3*b4
indBAEDSP := n1*l2*d3*g4*b5
indTotal := indBSP + indBDP + indBDSP + indBEP + indBEDP + indBEDSP + indBESP + indBAP + indBASP + indBADP + indBADSP + indBAEP + indBAESP + indBAEDP + indBAEDSP
TotalEffectAFE := ACME1 + l1
'
fit <- sem(model = mediation_model, data = data_all, se = "boot", bootstrap = 5000)
summary(fit, fit.measures = T, modindices = F)
parameterEstimates(fit, boot.ci.type = "bca.simple", level = 0.9)
```

## Best Answer

The modeled implication of fixing them to 0 is that 100% of the reason(s) that they covary are sharing common causes, all of which are among their predictors in your model. If that does not sound like a reasonable assumption, then you should estimate the residual covariances (FYI, the standardized counterpart in

`std.all`

is their partial correlation). If you want to treat is as an empirical question, then comparing the fit of models with(out) fixing them to zero provides a test of the $H_0$ that they are zero.