Residuals and Homoscedascity in SEM with lavaan I intend on investigating longitudinal effects through a cross-lagged panel design, facilitating the lavaan package in R. Currently, I have set up a simple Model given the scheme and code below, in accordance to the lavaan tutorial, as well as various lecture notes by the package's author. I am using a robust method as my input data is not on all variables normally distributed - lavaan is suited to handle this, there should be no problem. There are no latent variables, as there is only one observed variable each, so it is an autoregression, really.

library(lavaan)

model <- '
  # linear part
    A2 ~ A1 # a12
    A3 ~ A2 # a23
    B2 ~ B1 # b12
    B3 ~ B2 # b23
  # cross-lagged part
    A2 ~ B1 # d12
    A3 ~ B2 # d23
    B2 ~ A1 # c12
    B3 ~ A2 # c23
  # variance
    A1 ~~ A1
    A2 ~~ A2
    A3 ~~ A3
    B1 ~~ B1
    B2 ~~ B2
    B3 ~~ B3
  # covariance
    A1 ~~ B1
    A2 ~~ B2
    A3 ~~ B3
'

fit <- sem(model, data = Data, estimator = "MLR")
summary(fit, fit.measures = T)

Now even without constraints (such as a = a12 = a23 and such) the CFI is above 0.95 and P-value > 0.01 (not great, but it's a start). I suppose this means the autoregression did work.
So for checking for homoscedascity, I would need the residual plot (residuals vs. fitted), given a predicted dependent variable from my model distribution in compairison to residuals of my acutal data. I don't really know how I would go about generating a model distribution like that given in my CLP-model, let alone the fact that it is bivariate. Or is homoscedasticity no longer a necessary assumption for this kind of model?
Could anyone please explain this to me and maybe even help me solve this?
 A: Robust methods in SEM are sandwich estimators, and the standard errors are robust to violation of the homoscedasticity assumption. (You can try this by generating simple data and replicating a t-test in regression 
Here's some code:
library(lavaan)
set.seed(2020)
d <- data.frame(x = c(rep(0, 100), rep(1, 900)), 
                y1 = rnorm(1000))
d$y2 <- ifelse(d$x == 0, d$y1 / 10, d$y1)

with(d, t.test(y1 ~ x))
with(d, t.test(y1 ~ x, var.equal = TRUE))

with(d, t.test(y2 ~ x))
with(d, t.test(y2 ~ x, var.equal = TRUE))


fit1 <- "y1 ~ x"
summary(lavaan::sem(fit1, data = d, estimator = "ML"))
summary(lavaan::sem(fit1, data = d, , estimator = "MLR"))

fit2 <- "y2 ~ x"
summary(lavaan::sem(fit2, data = d, estimator = "ML"))
summary(lavaan::sem(fit2, data = d, , estimator = "MLR"))

Just looking at the p-values (I won't paste all the output).
First two t-tests, (where homogeneity of variance is not violated), and I do unequal (Welch's) t-test, which does not assume HoV, and equal variance, which does: 0.20 and 0.17 - pretty similar. Second two, where HoV is violated: 0.14, 0.60. Big difference.
First two p-values from lavaan, with ML and MLR: 0.14, 0.20 - pretty much the same. Second two, with HoV violated: 0.61, 0.14. Same as the t-test. 
