How to plot individual case residuals (ICRs) from lavaan model?

Question

I was reading this article on saturated models in structural equation modeling (SEM) and it appears they describe a tool for checking how well the model fits data even if the model fit indices indicate a "perfect fit" under raw data misspecification. They simply take the individual case residuals (ICRs) and plot them against a variable of interest or parameter not included in the original SEM. They can be seen below, which are just the scatterplots of these checks:

However, I was curious how one can obtain this information in packages like lavaan or otherwise, as the examples they give are in Stata. I am aware one can get residuals by running resid(fit) or something similar with lavResiduals, but these typically only give you residual matrices. Is there a way of obtaining ICRs in R?

Answer 1

Is there a way of obtaining ICRs in R?

Not directly, but there is a relatively new method for obtaining casewise predicted values in SEM:

de Rooij, M., Karch, J. D., Fokkema, M., Bakk, Z., Pratiwi, B. C., & Kelderman, H. (2023). SEM-based out-of-sample predictions. Structural Equation Modeling, 30(1), 132-148. https://doi.org/10.1080/10705511.2022.2061494

This is implemented in the new lavPredictY() function. Once you save predicted values, you can easily subtract those from observed values to get casewise residuals. Jarrett Byrnes recently shared some code for his own lavResidualsY() function that could help:

https://github.com/yrosseel/lavaan/issues/269