Efficient Bayesian Structural Equation Modeling in Stan

Structural equation models comprise a large class of popular statistical models, including factor analysis models, certain mixed models, and extensions thereof. Model estimation is complicated by the fact that we typically have multiple interdependent response variables and multiple latent variables (which may also be called random effects or hidden variables), often leading to slow and inefficient MCMC samples. In this paper, we describe and illustrate a general, efficient approach to Bayesian SEM estimation in Stan, contrasting it with previous implementations in R package blavaan (Merkle&Rosseel, 2018). After describing the approaches in detail, we conduct a practical comparison under multiple scenarios. The comparisons show that the new approach is clearly better. We also discuss ways that the approach may be extended to other models that are of interest to psychometricians.


Introduction
Structural equation models (SEMs) are commonly used in the social sciences, where it is customary to (attempt to) measure unobservable traits such as cognitive abilities, attitudes, and proficiencies. Such models provide a formal way of connecting these unobservable traits to related, observed variables (e.g., test scores, Likert responses, etc), which has led to increased popularity of SEMs. SEMs are also related to research on causality and directed acyclic graphs (e.g., Pearl 2013), to generalized linear mixed models (e.g., Bates, Mächler, Bolker, and Walker 2015;Gelman, Carlin, Stern, Dunson, Vehtari, and Rubin 2013;Stroup 2013), and to time series models (e.g., Driver, Oud, and Voelkle 2017), illustrating the models' broad applicability across disciplines.
A defining feature of the SEM framework is the ability to instantiate regressions on latent variables, as opposed to observed variables. This framework is more general than the traditional mixed modeling framework, allowing for products between latent variables and other free parameters, as might be seen in factor analysis (e.g., Bollen 1989;Merkle and Wang 2018). The generality of SEM implies that the estimation methods are relatively complex, which has historically led researchers to rely on closed-source implementations of optimization methods via software like Mplus Muthén 1998-2017), LISREL (Jöreskog and Sörbom 1997), and EQS (Bentler 2000(Bentler -2008. A small number of more recent R (R Core Team 2019) packages, including sem (Fox, Nie, and Byrnes 2017a), OpenMx (Boker, Neale, Maes, Wilde, Spiegel, Brick, Spies, Estabrook, Kenny, Bates, Mehta, and Fox 2011), and lavaan (Rosseel 2012), provide open source SEM functionality that utilize classical estimation methods including maximum likelihood or least squares.
While maximum likelihood and least squares methods are most popular, Bayesian approaches to SEM and related models have received increased recent attention (e.g., Depaoli and van de Schoot 2017;Fox 2010;Jackman 2009;Kaplan 2014;Kruschke 2011;MacCallum, Edwards, and Cai 2012;Merkle and Wang 2018;Muthén and Asparouhov 2012;van Erp, Mulder, and Oberski 2018). Researchers have specifically found the methods to be useful for estimation of complex SEMs (including, e.g., latent variable interactions; Lee, Song, and Tang 2007), for automatically handling uncertainty associated with latent variable estimation, and for scaling to high-dimensional datasets.
Despite the increased popularity of Bayesian latent variable models, coding the models via JAGS (Plummer 2003) or Stan (Stan Development Team 2017) syntax can be difficult, and the resulting sampling can be time-consuming. These issues have been partially addressed by R package blavaan (Merkle and Rosseel 2018), which uses lavaan model specification syntax and originally relied on JAGS for model estimation (via package runjags, which provides an R interface to JAGS; see Denwood 2016). Other R packages have addressed these issues for related models, including brms (Bürkner 2017) for mixed and multivariate models, ctsem (Driver et al. 2017) for time series models, edstan (Furr 2017) for item response models, and pcFactorStan (Pritikin and Furr 2019) for pairwise comparison factor models.
The original blavaan approach was similar to the brms approach for generalized linear mixed (and related) models, where JAGS code was generated at runtime from the user-specified model syntax. However, this approach became very slow for some models, forcing the user to wait hours or more for enough samples to make inferences. This reduced the viability of blavaan for applied data analysis and for simulation studies, leading us to implement Stan functionality in blavaan. The original Stan implementation was similar to the JAGS implementation, generating Stan syntax for a user-specified model and relying on package rstan (Stan Development Team 2018) for MCMC. This Stan implementation has not been formally described, which represents one contribution of the current paper. In general, though, the Stan implementation was not much faster or more efficient than the JAGS approach.
The primary contribution of this paper is to describe and illustrate a new approach to Stan SEM estimation, which greatly improves the speed and efficiency of model estimation. The approach can be flexibly applied to models in the traditional SEM framework, with general (possibly non-conjugate) prior distributions. It has been implemented in blavaan alongside the previous JAGS and Stan implementations, allowing for easy comparison across implementations.
In the sections below, we first formally define the models under consideration. We then describe the three MCMC approaches that are now implemented in blavaan: the original JAGS approach described in Merkle and Rosseel (2018), the original Stan approach that is being formally described in this paper for the first time, and the new Stan approach that is the primary focus of this paper. After describing the approaches, we explicitly discuss some problematic issues associated with estimation of SEMs via MCMC. These are issues that are often overlooked in the literature, but they are necessary to have fully functional Bayesian SEM software. Finally, we compare the three approaches via three examples, highlighting the advantages of the new Stan approach.

Model definition
The blavaan package generally relies on the lavaan representation of a structural equation model, which is based on the LISREL "all-y" representation (e.g., Jöreskog and Sörbom 1997).
Let y i be the p (continuous) observed variables associated with observation i. Then a structural equation model with m latent variables may be represented by the equations where η i is an m × 1 vector containing the latent variables; i is a p × 1 vector of residuals; and ζ i is an m × 1 vector of residuals associated with the latent variables. The vectors ν and α contain intercept parameters for the manifest and latent variables, respectively; Λ is a matrix of factor loadings; and B contains parameters that reflect directed paths between latent variables.
The residuals i and ζ i are assumed to be multivariate normal: where the associated covariance matrices are often diagonal. These assumptions imply that the marginal distribution of y (integrating out the latent variables) is multivariate normal with parameters which requires that (I − B) be invertible. The traditional LISREL framework includes additional matrices for exogenous observed variables, but these are not often utilized in lavaan. Instead, an exogenous observed variable is "upgraded" to latent variable status, where the latent variable accounts for all of the observed variable's variance (and the associated variance parameter in Θ is fixed to 0).
Many Bayesian approaches to SEM estimation rely on sampling the η i in tandem with other model parameters. This is advantageous because observed variables are often independent conditioned on the η i , so that the conditional distribution of each observed variable is a univariate normal. However, as we will see later, the sampling of the η i can have a major impact on the speed and efficiency of MCMC estimation.

MCMC approaches
In the sections below, we briefly describe the three MCMC approaches implemented in blavaan: the original JAGS approach, the original Stan approach, and the new Stan approach that is the focus of this paper.

Parameter expansion in JAGS
In previous work (Merkle and Rosseel 2018), we developed a parameter expansion approach that can be applied to SEMs for continuous data (also see Palomo, Dunson, and Bollen 2007).
The method allows researchers to place prior distributions on intuitive sets of parameters and can be generally implemented in JAGS.
The approach involves converting the model of interest (the "inferential model") to an overparameterized, equivalent model (the "working model") from which it is easier to sample. The conversion focuses on a model's covariance parameters, converting each covariance to a "phantom" latent variable. This conversion makes observed variables conditionally independent of one another (conditioned on latent variables), meaning that our likelihood involves a series of univariate distributions instead of a single multivariate distribution. Such a conversion can speed up sampling in JAGS, where computations involving the multivariate normal distribution are very slow. The full details underlying these procedures can be found in Merkle and Rosseel (2018).

Likelihood simplification in Stan
The phantom latent variable approach is not essential in Stan and, in testing, we found that the approach did not lead to gains in sampling speed or efficiency. However, we did make initial progress in Stan by capitalizing on the structure of the SEM latent variable covariance matrix Ψ. This capitalization was inspired by related work on estimating multivariate autoregressive models in Stan (Joseph 2016).
We provide an overview of this approach here. For traditional SEMs, the distribution of latent variables can typically be expressed as (see Equation (2)) Evaluation of this multivariate normal log-likelihood is time-consuming in Stan because we need to compute the inverse and determinant of the covariance matrix. However, for many models, the structure of the covariance matrix leads to simplifications. For example, the matrix B is often triangular with zeros along its diagonal (leading to a so-called "recursive" model), and the matrix Ψ is often diagonal. When both of these properties are fulfilled, we can use standard matrix properties (e.g., Petersen and Pedersen 2012) to write the determinant as a product of scalar values: Relatedly, the inverse of Σ is simplified as which completely removes the need to compute matrix inversions when Ψ is diagonal. When either B is triangular or Ψ is diagonal (but not both), we can use a subset of the above simplifications to improve sampling efficiency as much as possible. Given a specific model, package blavaan automatically determines which simplifications are available and uses them for Stan estimation. The simplifications are implemented in Stan as a custom log-probability density function. This implementation is available in blavaan via the argument target = "stanclassic".

New Stan approach
Both methods mentioned above exploit the fact that the latent variables in the model can be sampled along with other model parameters. This generally simplifies likelihood computations and allows us to immediately extend the methods to situations where observed variables have non-normal distributions. Most Bayesian approaches to SEM, and to other models with "random" parameters, sample the latent variables.
However, the sampling of latent variables greatly increases the dimension of the parameter space, which can reduce sampling speed and efficiency. Unlike JAGS, the key to fast sampling in Stan is to work with a model likelihood that is marginal over latent variables. This is somewhat unintuitive, because previous researchers have focused on the simplifications that we can gain from sampling the latent variables. One concern related to using the marginal likelihood involves our inability to make inferences about the latent variables (because they are integrated out of the likelihood). But this concern is addressed by blavaan because, conditional on the other model parameters, the latent variable posterior distribution is tractable. Thus, the latent variables can be sampled in a "generated quantities" block within the Stan syntax, even though they do not directly play a role in the MCMC sampling.
The new blavaan approach utilizes the marginal likelihood, and we have written a single Stan program that can estimate the majority of multivariate normal SEMs that a user could specify. This file is compiled once during (or before) package installation. Then, once the user specifies a model, many pieces of information about the data and about the model are passed to the compiled model, with sampling occurring immediately. It is inconvenient to enter all the required model information by hand. Thus, to complement the Stan file, we have new R code that serves as a pipeline from blavaan to the Stan model and back. The blavaan user will not notice many differences, because the commands for model specification and estimation are the same as before. However, the model is now sent to the pre-compiled Stan code by default, whereas the previous approaches wrote JAGS or Stan code at runtime.
It is worth noting that our Stan SEM file stands on its own, so that users of languages beyond R (e.g., Python) could also utilize the file if they can pass all the required data in to the Stan model. This is more challenging than it may sound due to the many pieces of data that are required, including the dimensions of all SEM matrices, the free entries of SEM matrices, equality constraints on free parameters, prior distribution parameters, and so on. Because our Stan model is precompiled, the possible models that can be estimated are restricted in two ways. First, there is some inflexibility in choice of prior distributions. For most types of model parameters, the form of each parameter's prior cannot be changed (though the specific prior hyperparameters can). For example, regression parameters (B) in blavaan currently have N(0, 10) priors by default, where the normal distribution is parameterized by standard deviation. Users can change the mean or standard deviation of this normal prior, but they cannot change the fact that the prior is normal. However, for scale parameters, users have the option to place priors on variances, standard deviations, or precisions.
The second restriction of the precompiled Stan model involves equality constraints. While our code currently allows for equality constraints within a class of parameters (e.g., loadings can be constrained equal to one another or intercepts can be constrained equal to one another), it does not allow for constraints across classes of parameters. Additionally, if users wish to set one parameter equal to a function of other parameters, that is not currently possible. However, if users desire features that are not included in our current implementation, they can take our Stan file, make the desired changes, and recompile the model. Alternatively, they could use the original MCMC methods available in blavaan, which provide more flexibility because they are not precompiled.

Challenging issues
While the above methods can be readily applied to "vanilla" models such as confirmatory factor analysis with uncorrelated factors, the SEM framework includes many model and data characteristics that require further attention for estimation. Below, we highlight three characteristics requiring special attention.

Covariance parameters
The general SEM presented earlier includes two covariance matrices with free parameters: Θ and Ψ. These matrices can include some fixed values and some free values, so prior distributions for these matrices are not straightforward. That is, there are some models for which we cannot simply place an inverse Wishart prior on the covariance matrix, nor an LKJ prior (Lewandowski, Kurowicka, and Joe 2009). Those priors were meant for unrestricted covariance/correlation matrices, not for matrices with some fixed values and some free values.
In the Stan approaches, we consequently decompose the covariance matrices into standard deviations and correlations. For example, Θ is written as where D Θ is a diagonal matrix of standard deviations and R Θ is a correlation matrix. Prior distributions are then placed individually on the free standard deviation parameters and on the free correlation parameters within the two matrices. This approach is similar to that of Barnard, McCulloch, and Meng (2000), and Liu, Zhang, and Grimm (2016) provide a comparison of this approach to the use of inverse Wisharts in the context of growth curve models.
The use of an independent prior on each free parameter can sometimes lead to a non-positive definite covariance matrix during MCMC sampling. Stan is able to reject such a covariance matrix and continue sampling, whereas JAGS will terminate. This is why we developed the parameter expansion method in JAGS: the parameter-expanded model involves diagonal covariance matrices that cannot become non-positive definite. But the non-positive definite covariance matrices still have implications for model calibration, as detailed in the simulationbased calibration study later.
We are aware of a variety of other prior distributions proposed for covariance matrices (e.g., Chung, Gelman, Rabe-Hesketh, Liu, and Dorie 2015;Consonni and Veronese 2003;Mulder and Pericchi 2018;Spezia 2018). The strategy implemented in blavaan is worthwihle because it is relatively easy to specify informative prior distributions for individual standard deviation and correlation parameters in the model. In contrast, many of the other prior distributions are proposed for convenience or due to the fact that they maintain positive definiteness, and they have less-intuitive interpretations as compared to our approach. We plan to further consider these alternative priors in the future.

Missing data
While it is often useful and desirable to directly model the missing values with the rest of the model (e.g., Merkle 2011; O'Muircheartaigh and Moustaki 1999), blavaan employs a "missing at random" approach to missing data that differs across JAGS and Stan. In JAGS, one can include NA values in the data, and JAGS will sample these missing values as if they were extra model parameters. In contrast, Stan does not allow NA values in the data, so that one must handle the missing data manually. We utilize a "full information" likelihood (e.g., Wothke 2000) in our Stan models, which is the same likelihood that is used to handle missing data in lavaan and other software that performs maximum likelihood SEM estimation. This requires some additional overhead in preparing the data to be sent to Stan, because each case's observed values must be indexed, and cases are sorted by missing data pattern to speed up computations. Missing values could also be directly sampled ("imputed") in Stan, though this functionality is not currently available.

Latent variable scaling
Structural equation models typically require some parameter constraints to achieve parameter identification, where we must "set the scale" of each latent variable. The two most popular ways to do this involve (i) fixing each latent variable's variance to 1, or (ii) fixing a single loading (parameter in λ) to 1 for each latent variable. Of these two, the latter method is most straightforward to implement in a Bayesian setting.
The former method (of fixing each latent variance to 1) is more challenging. This is because, as described by Peeters (2012), one loading per latent variable must be sign constrained to achieve global parameter identification. Otherwise, the sign of each loading may flip back and forth, with a model's regression parameters and covariance parameters potentially flipping along with the loadings. One solution to this issue involves the placement of a truncated normal prior (truncated from below at 0) on one loading per latent variable, preventing the sign changes. This solution is adopted in blavaan's JAGS approach to model estimation.
A different solution is implemented in blavaan's Stan approaches. In those approaches, the sign flipping is allowed to occur during MCMC sampling. The issue is then handled after sampling, in the "generated quantities" block. In this block, one loading per latent variable is transformed to always be positive, and the signs of associated parameters (loadings, regressions, and covariance parameters) are flipped every time the sampled value of the focal loading is negative. This approach can improve sampling efficiency because no boundary constraints are introduced in the parameter space. This approach was discussed in a thread on the Stan Discourse site (https://discourse.mc-stan.org/t/latent-factor-loadings/1483).

Applications
The estimation approaches described above are all implemented in package blavaan for general SEM estimation. These include the original JAGS approach (obtained via argument target = "jags"), the original Stan approach (target = "stanclassic"), and the new Stan approach (target = "stan").
For example, the following code specifies a model for the well-known "political democracy" data (Bollen 1989) and estimates it via each of the three approaches. This dataset includes 75 countries measured on 11 attributes, seven of which were measured in 1960 and 4 of which were measured in 1965. The intent of the model is to study relationships between countries' levels of industrialization and democracy over time.
Following model estimation, convergence diagnostics such as Rhat and effective sample size are immediately available via the summary method and blavInspect() function, and many types of plots are available via the plot method, which relies on package bayesplot (Gabry and Mahr 2019). Further examples of blavaan syntax and functionality can be found in Merkle and Rosseel (2018), noting that target="jags" was the default at the time that paper was written, while target="stan" is now the default.
In the following sections, we compare the three MCMC approaches on speed and sampling efficiency via three example models. We then study the extent to which the best method (the marginal Stan method) is calibrated.

Performance Comparison
All comparisons are carried out on a Dell desktop with a large amount of RAM, running Ubuntu Linux. We define sampling efficiency as "effective sample size per second" (ESS/s). Effective sample sizes are computed via the rstan monitor() function, and sampling time is measured after Stan model compilation. The warmup time for Stan models was fixed to 300 iterations, whereas the burn-in time for JAGS models was fixed to 1000 iterations. There is arbitrariness in the warmup and burn-in choices, so that the ESS/s metric is somewhat crude. But we think the metric is sufficient to illustrate the advantage of the new Stan approach. We examine the MCMC methods' speed and efficiency on three models, two of which are popular models often used to illustrate SEM methods. The third is a more complex model that is known to pose difficulties for the original blavaan approach. We do not conduct a full Monte Carlo study here, so our results are subject to noise. But the results are generally consistent across the models presented here as well as many others not presented, so we think they can be taken as general evidence for the approaches' relative performance. The results are also consistent with those of Yackulic, Dodrill, Dzul, Sanderlin, and Reid (in press), who study marginalization of discrete latent variables in ecological models.

Political democracy
For our first comparison, we continue with the Bollen (1989) political democracy model. The blavaan code to fit the model was shown earlier. The JAGS method was fastest here, averaging 0.55 seconds per 100 iterations. Next fastest was the new (marginal) Stan method, averaging 1.44 seconds per 100 iterations, followed by the old Stan method at 7.36 seconds per 100 iterations. But it is more important to examine the methods' sampling efficiencies (effective sample size per second), which are shown in Figure 1. A separate metric is shown for each parameter, with the parameters being numbered along the x-axis. Parameters are ordered on the x-axis by parameter type, with the details being shown in the figure caption. The figure shows that the speed of JAGS is offset by the effective sample size, so that the new Stan method is best in terms of sampling efficiency. The old Stan method exhibits efficiency similar to that of JAGS.

Holzinger and Swineford
Our second example involves a confirmatory factor analysis of the Holzinger and Swineford (1939) data. This is the version of the dataset included in package lavaan, which has 301 individuals measured on nine cognitive scales. The confirmatory factor model fit to the data includes three latent variables, each of which is associated with three observed variables. The blavaan code to specify and fit the model is R> HS.model <-visual =~x1 + x2 + x3 + textual =~x4 + x5 + x6 + speed =~x7 + x8 + x9 R> fit <-bcfa(HS.model, data = HolzingerSwineford1939) where additional arguments would typically be used to specify the number of sampling iterations, to specify priors, to specify the MCMC sampler, and so on.
In terms of speed, the JAGS method is again fastest, averaging 0.91 seconds per 100 iterations. This was followed by the marginal Stan method at 2.83 seconds per 100 iterations, then the old Stan method at 9.55 seconds per 100 iterations. The ESS/s metrics for this model are visualized in Figure 2. The graph is similar to that of the previous section, with fewer parameters in this model as compared to the last model. We see that the gold line, representing the new Stan method, is the highest for all the model parameters, being at least twice as large as the other methods' efficiencies. The JAGS and old Stan methods are again similar to one another for this example, with JAGS being better for the majority of parameters.
The sampling speed is now reversed, with the marginal Stan method at 23.4 seconds per 100 iterations, the JAGS method at 28.4 seconds per 100 iterations, and the old Stan method at 577.27 seconds per 100 iterations. If we instead compute the JAGS speed while accounting for thinning (i.e., counting only each twentieth iteration in the computations), then the JAGS speed is at 567.98 seconds per 100 iterations.
The approaches' sampling efficiencies are shown in Figure 3, where parameter ordering is again described in the figure caption. The JAGS and old Stan methods are very low for this model, with the new Stan method displaying much better efficiency and yielding useful results in a matter of minutes, as opposed to hours or days. This and related examples (not shown) convinced us to make the new Stan method the default in blavaan, replacing the original default of JAGS. The new Stan method reliably produces fast, efficient samples for a large number of models, whereas the other methods exhibit more variability in their speeds and efficiencies, and are seldom clearly better than the new Stan method.
under weak priors to the maximum likelihood estimates and standard errors from lavaan. For all three MCMC methods, the posterior means obtained under this approach typically agree with lavaan estimates to about one decimal point. The posterior standard deviations tend to be close to, but slightly larger than, the maximum likelihood standard errors. Now that there are multiple MCMC methods implemented in blavaan, we have also been able to compare MCMC methods to one another in order to verify that they were producing similar posterior distributions.
Here, we use the simulation-based calibration method proposed by Talts, Betancourt, Simpson, Vehtari, and Gelman (2018) to study the calibration of blavaan's new (marginal) Stan implementation. This involved repeatedly generating data from the model's prior distribution, fitting the model to the generated data, and examining the ranks of the posterior MCMC samples relative to samples from the prior distribution. If the MCMC algorithm is calibrated, then these ranks should be approximately uniformly distributed. Deviations from uniformity are then taken as miscalibration.

Method
Our simulation-based calibration study utilized the political democracy model presented earlier. We generated 500 datasets of size 75 from the prior predictive distribution and fit the model to each generated dataset via MCMC. The study involved two conditions that differed by the prior distributions that were used. First, we used the default prior distributions from blavaan, which are intended to be weakly informative in many situations encountered in practice. Second, we used a set of more informative prior distributions to contrast with the noninformative priors. Both sets of prior distributions are shown in Table 5.2.

Results
Rank frequencies for the blavaan default priors are shown in Figure 4. Perhaps surprisingly, the distributions are far from uniform, with peaks generally occurring near zero and one. These peaks represent posterior distributions that exhibit less variability than they should, given the non-informative priors from which we started. Clearly, the results are far from the uniformity that would be expected from a calibrated algorithm.
The non-uniformity occurs because our model has a large number of parameters with independent prior distributions, with the parameters all contributing to the model-implied covariance matrix (see Equation (6)). For the model considered here, there are many combinations of parameters that lead to a non-positive definite covariance matrix, and the MCMC sampler will avoid these combinations of parameters during sampling. In the simulation-based calibration study, this leads to posterior samples that are not calibrated with respect to the independent priors. Instead, we might say that the posteriors are calibrated with respect to regions of the prior distribution that are positive definite.
To provide evidence that the MCMC algorithm is indeed calibrated with respect to priors that maintain positive definiteness, we show the results of the informative priors in Figure 5. These frequencies are now much closer to uniform, because the information in the prior distributions now generally leads to positive definite model covariance matrices. This interplay between the informativeness of prior distributions and posterior calibration is worthy of further attention, because existing MCMC algorithms for SEM use a series of independent priors on parameters that each play a role in the model-implied covariance matrix. A researcher's priors may be more informative than expected, based solely on the fact that the model-implied covariance matrix must remain positive definite during MCMC sampling. Further, depending on the model likelihood used (marginal vs. conditional), the degree of information present in uninformative priors may vary. These results highlight the utility of blavaan for conducting detailed study of MCMC algorithms, as well as the fact that there is room to improve the default blavaan priors in the future.

Conclusion
The results in this paper show that we can improve sampling efficiency by integrating the latent variables out of the model likelihood, which is the opposite of most popular approaches to Bayesian SEM estimation (where the popular approaches are largely based on results summarized by, e.g., Lee 2007;Song and Lee 2012). We can expect the marginal sampling efficiency to be even more advantageous as sample sizes increase, because the sample size has no impact on the dimension of the parameter space here. In contrast, the dimension of the parameter space increases with sample size under conditional approaches, where latent variables count as parameters.
While the marginal approach is promising, use of the marginal likelihood leads us back to problems that frequentists often encounter in SEM. These problems include the fact that the marginal likelihood does not have a closed form when we have non-normal observed variables (e.g., ordinal variables) or when we have latent variable interactions. We think that some progress can be made here by employing other Bayesian methods, including data augmentation (e.g., Chib and Greenberg 1998) in the ordinal case. The use of data augmentation for psychometric models has been described by Fox (2010) and Fox, Mulder, and Sinharay (2017b), and such methods may be implemented in future versions of blavaan.
For situations where the marginal likelihood does not exist in closed form, it is also possible to move back to the original blavaan approaches that sample the latent variables. However, in our experience, the original approaches are even slower and less efficient in those situations (as compared to the models considered here), making them questionable for applied work. Further, even if those methods did exhibit reasonable efficiency, the marginal likelihood is generally necessary for obtaining suitable information criteria such as DIC (Spiegelhalter, Best, Carlin, and Linde 2002) or WAIC (Watanabe 2010). Merkle, Furr, and Rabe-Hesketh (2019) discuss why the marginal likelihood is preferable here, and Zhang, Tao, Wang, and Shi (2019) discuss related applications of DIC to multilevel item response models. Thus, we think that use of the new Stan approach, paired with new tricks for handling non-normal observed variables, is the most viable approach for applications to the non-normal modeling situations typically encountered in practice.    We should note that the three methods here are not the only ones that can be conceptualized in JAGS or in Stan. For example, the pre-compiled marginal Stan approach can be modified so that the latent variables are part of the model likelihood. This leads to a method that is similar to the "old Stan" method, except that it simplifies computation of the model likelihood in different ways. In limited testing, we found that this approach was somewhat closer to the JAGS approach in sampling efficiency, but still similar to the "old Stan" approach studied in this paper. Alternatively, it is possible to define a marginal approach in JAGS, but our limited testing there indicates that use of the multivariate normal distribution leads to major decreases in JAGS sampling efficiency.
Finally, our analyses here have included no use of parallelization. While between-chain parallelization is immediately available in blavaan, there have been recent advances in within-chain parallelization in Stan. These advances are promising for further improvement of Stan sampling efficiency, and we plan to include this in future versions of blavaan.