# Estimating model fit in multiple mediator SEM

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)

C1 <- rbinom(10, 1, 0.5)
C2 <- rnorm(100)

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