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.
set.seed(13012) library(lavaan) # 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) fitmeasures(fit) OUTPUT: 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