I'm building a multiple mediation model in lavaan. My research question concerns to what extent a relationship between a continuous X and a binary Y is mediated by two continuous mediators M1 and M2. I have large observational datasets (N > 5000) and I include a bunch of controls (5-6 continuous or categorical variables, the latter recoded as dummies), which may be reasonably assumed to affect both the mediators and the outcome and thus could bias the estimates. Now, following Preacher and Hayes (2008) I allow the mediator residuals' covariance to be freely estimated. This also makes theoretical sense, because there could be things which influence both, that I cannot control for.

My issue is that my model appears to be just identified i.e. saturated and thus with no degrees of freedom, tests of model fit are not estimated. Neither can I rely on comparative fit indicies (AIC, BIC), because of my categorical DV.

The question of course is, what to do.
- I've thought about trying to justify a saturated model on theoretical grounds, but I don't think I could get away with this.
- I've thought about fixing the mediator residual covariance to 0 just for the sake of estimating model fit, but it affects my estimates substantially and is an unreasonable restriction.
- I've thought about fixing constrain two 'unimportant' coefficients to be equal. I've got a few categorical variables rescaled as a series of dummies. I could claim something like let's assume people living in small vs mid sized towns are the same. I don't think this should affect my models, but this seems awfully arbitrary.

I'm writing my analysis for a social science audience. I've tried to look at published multiple mediator SEM models, but they do not really discuss overidentifying restrictions (?), which is strange because no matter how I look at this, saturation seems to be inherent to these kidns of models. Thanks for the help.


# independent variable of interest
X <- rnorm(100)

# add two controls 
C1 <- rbinom(10, 1, 0.5) 
C2 <- rnorm(100)

# add mediators
M1 <- 0.5*X + 0.1*C1 + rnorm(100)
M2 <- 0.35*X + 0.1*C1 + 0.3*C2 + rnorm(100)

# linear combinations of variables
z <-  0.5* M1 + 1*M2 + 0.5 * X + 2*C1 + 0.25 *C2 + rnorm(100) 

pr <- 1/(1+exp(-z)) # inverse link function 

# generate outcome of interest
Y <- rbinom(100,1,pr)

# combine to df
Data <- data.frame(X = X, Y = Y, M1 = M1, M2 = M2, C1 = C1, C2 = C2)

# specify SEM model
multipleMediation <- '
        # outcome model 
        Y ~ b1 * M1 + b2 * M2 + c * X + z1 * C1 + z2 * C2

        # mediator models 
        M1 ~ a1 * X + z3 * C1 + z4 * C2
        M2 ~ a2 * X + z5 * C1 + z6 * C2

        # allow mediator residuals to covary
        M1 ~~ M2

        # estimate indirect & total effects, contrast & proportion mediated 
        indirect1 := a1 * b1
        indirect2 := a2 * b2
        total    := c + (a1 * b1) + (a2 * b2)
        contrast := indirect2-indirect1
        prop1 := indirect1/total
        prop2 := indirect2/total
fit <- sem(model = multipleMediation, 
           ordered = "Y",
           estimator = "dwls", 
           se = "bootstrap", 
           bootstrap = 10,
           data = Data)

          npar                fmin               chisq                  df              pvalue 
         17.000               0.000               0.000               0.000                  NA 
 baseline.chisq         baseline.df     baseline.pvalue                 cfi                 tli 
         68.855              12.000               0.000               1.000               1.000 
           nnfi                 rfi                 nfi                pnfi                 ifi 
          1.000               1.000               1.000               0.000               1.000 
            rni               rmsea      rmsea.ci.lower      rmsea.ci.upper        rmsea.pvalue 
          1.000               0.000               0.000               0.000                  NA 
            rmr          rmr_nomean                srmr        srmr_bentler srmr_bentler_nomean 
          0.000               0.000               0.000               0.000               0.000 
    srmr_bollen  srmr_bollen_nomean          srmr_mplus   srmr_mplus_nomean               cn_05 
          0.000               0.000               0.000               0.000               1.000 
          cn_01                 gfi                agfi                pgfi                 mfi 
          1.000               1.000               1.000               0.000               1.000  

1 Answer 1


This is an old question but I think it could be useful for some users. Every time I talked to SEM experts, the answer is that there is no a "one-way" solution. As long as you a transparent on each step explaining pros and cons of each model tested you can present a saturated model, the point is that I would not only present this model. You could test and discuss different model structures. For example, I would save a df by omitting the direct relationship between X and the outcome variabile, but I would also present the model that does not include the residual correlation between the mediators.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.