MLSEM for categorical outcomes requires numerical integration, which is quite slow using R's built-in optimizers. Thus, it will be a while before this is possible in lavaan. But there are a couple other possibilities, potentially. Can you tell us more about your model and data?

• how many subjects nested in how many clusters?
• how many indicators for (each of) how many factors, with how many thresholds?

Thanks for the comments and willingness to investigate on possibilities even if it is currently not implemented and is not planed to implement in lavaan in the nearest future. More context about my task is given below.

I have longitudinal/panel data with multiple observations over time of the same subjects. So each subject is a "cluster" in this case. There are also other clustering structure present as well e.g. some subjects belonging to the same structural group (or cluster in your terminology). The minimum I would like to do is to account somehow for the correlation within subjects (different observations of the same subject).

In my real task there are ~100k subjects with ~15-20 measurements per subject. There are ~20-30 indicators (x's) and ~5-7 factors I am considering. Most x's are already "embedded" in some way and are mapped to R (S->R), originally being ordinal or categorical. So I think x's could be treated as continuous in the worst case. Y is an endogenous (observable) variable of the primary interest and it is ordinal with 7 levels (6 underlying thresholds).
Levels of y are imbalanced meaning that there are much more observations with y = l0 compared to other levels. In the worst case binning some levels together to form a binary variable might be an option, but keeping original 7 levels would make the analysis much richer.

An illustration on the simple 2-level model is given below. I guess it could be something to start with if y would be continuous, however it is ordinal.

model <-'
level: 1
L1 =~ x1 x2 x3
L2 =~ x4 x5
y ~ L1 L2
level: 2
y ~ 0
y~~ y
'

I would be very interested in obtaining cluster-robust standard errors when using ordinal data (and WLSMV).

I am not sure if I will ever implement random effect models with ordinal data in lavaan. It would just be too slow (in R) to be practical. But Nick Rockwood is working on glavaan, which will have a backend in C . The idea of glavaan is to support many outcome types (binary, ordinal, ...) similar to generalized linear models. Eventually (but please do not ask when) he may also add support for multilevel models.

As for cluster-robust standard errors for ordinal data WLSMV: can Mplus do this? Is there any documentation about this?

Yes, I believe Mplus supports cluster-robust standard errors for ordinal data WLSMV. I'm not sure about documentation, but I created an example with manifest variables (though it appears to work with latent variables, as well), as adapted from here. See the attached files for the Mplus data file, and the input and output files for 1) a regular ordinal model and 2) an ordinal model with cluster-robust standard errors.
Mplus files.zip

Thanks. That was helpful. I have not found papers/literature about this yet, but I do have some ideas how this can be done: when we compute the 'cluster robust' standard errors, we must recompute the 'WLS.W' matrix (see lav_muthen1984) taking the clustering into account, in order to compute a cluster-robust 'Gamma' matrix, which is then used in lav_model_nvcov_robust_sem(). At least, that is my current hypothesis.

Is this of any help?
https://www.statmodel.com/download/JSM2007000746.pdf

@sda030 Thanks. This document is about twolevel SEM with WLS. Unfortunately (as always with Mplus documents), they are a bit vague on the details. It might be feasible to implement this, but I think I will focus on the 'cluster-robust' version first.