blavaan: Bayesian structural equation models via parameter expansion

This article describes blavaan, an R package for estimating Bayesian structural equation models (SEMs) via JAGS and for summarizing the results. It also describes a novel parameter expansion approach for estimating specific types of models with residual covariances, which facilitates estimation of these models in JAGS. The methodology and software are intended to provide users with a general means of estimating Bayesian SEMs, both classical and novel, in a straightforward fashion. Users can estimate Bayesian versions of classical SEMs with lavaan syntax, they can obtain state-of-the-art Bayesian fit measures associated with the models, and they can export JAGS code to modify the SEMs as desired. These features and more are illustrated by example, and the parameter expansion approach is explained in detail.


Introduction
The intent of blavaan is to implement Bayesian structural equation models (SEMs) that harness open source MCMC samplers (in JAGS; Plummer 2003) while simplifying model specification, summary, and extension. Bayesian SEM has received increasing attention in recent years, with MCMC samplers being developed for specific priors (e.g., Lee 2007;Scheines, Hoijtink, and Boomsma 1999), specific models being estimated in JAGS and BUGS (Lunn, Jackson, Best, Thomas, and Spiegelhalter 2012;Lunn, Thomas, Best, and Spiegelhalter 2000), factor analysis samplers being included in R packages bfa (Murray 2014) and MCMCpack (Martin, Quinn, and Park 2011), and multiple samplers being implemented in Mplus (Muthén and Asparouhov 2012). These methods have notable advantages over analogous frequentist methods, including the facts that estimation of complex models is typically easier (e.g., Marsh, Wen, Hau, and Nagengast 2013) and that estimates of parameter uncertainty do not rely on asymptotic arguments. Further, in addition to allowing for the inclusion of existing knowledge into model estimates, prior distributions can be used to inform underidentified models, to average across multiple models in a single run (e.g., Lu, Chow, and Loken 1999), and to avoid Heywood cases (Martin and McDonald 1975).
The traditional SEM framework (e.g., Bollen 1989) assumes an n × p data matrix Y = (y 1 y 2 . . . y n ) with n independent cases and p continuous manifest variables. Using the LISREL "all y" notation (e.g., Jöreskog and Sörbom 1997), a structural equation model with arXiv:1511.05604v2 [stat.CO] 17 Nov 2016 m latent variables may be represented by the equations y = ν + Λη + (1) where η is an m × 1 vector containing the latent variables; is a p × 1 vector of residuals; and ζ is an m × 1 vector of residuals associated with the latent variables. Each entry in these three vectors is independent of the others. Additionally, ν 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 (with the assumption that (I − B) is invertible).
In conventional SEMs, we assume multivariate normality of the and ζ vectors. In particular, with the latter equation implying multivariate normality of the latent variables. Taken together, the above assumptions imply that the marginal distribution of y (integrating out the latent variables) is multivariate normal with parameters If one is restricted to conjugate priors, then the Song and Lee (2012) or Muthén and Asparouhov (2012) procedures are often available for simple MCMC estimation of the above model. However, if one wishes to use general priors or to estimate novel models, then there is the choice of implementing a custom MCMC scheme or using a general MCMC program like BUGS, JAGS, or Stan (Stan Development Team 2014). These options are often timeconsuming and difficult to extend to further models. This is where blavaan is intended to be helpful: it allows for simple specification of Bayesian SEMs while also allowing the user to extend the original model. In addition to easing Bayesian SEM specification, the package includes a novel approach to JAGS model estimation that allows us to handle models with correlated residuals. This approach, further described below, builds on the work of many previous researchers: the approach of Lee (2007) for estimating models in WinBUGS; the approach of Palomo, Dunson, and Bollen (2007) for handling correlated residuals; and the approaches of Barnard, McCulloch, and Meng (2000) and Muthén and Asparouhov (2012) for specifying prior distributions of covariance matrices.
Package blavaan is potentially useful for the analyst that wishes to estimate Bayesian SEMs and for the methodologist that wishes to develop novel extensions of Bayesian SEMs. The package is a series of bridges between several existing R packages, with the bridges being used to make each step of the modeling easy and flexible. Package lavaan (Rosseel 2012) serves as the starting point and ending point in the sequence: the user specifies models via lavaan syntax, and objects of class blavaan make use of many lavaan functions. Consequently, writings on lavaan (e.g., Beaujean 2014;Rosseel 2012) give the user a good idea about how to use blavaan.
Following model specification, blavaan examines the details of the specified model and converts the specification to JAGS syntax. Importantly, this conversion often makes use of the parameter expansion ideas presented below, resulting in MCMC chains that quickly converge for many models. We employ package runjags (Denwood 2016) for sampling parameters and summarizing the MCMC chains. Once runjags is finished, blavaan organizes summaries and computes a variety of Bayesian fit measures.

Bayesian SEM
As noted above, the approach of Lee (2007; see also Song and Lee 2012) involves the use of conjugate priors on the parameters from Equations (1) to (4). Specific priors include inverse gamma distributions on variance parameters, inverse Wishart distributions on unrestricted covariance matrices (typically covariances of exogenous latent variables), and normal distributions on other parameters.
Importantly, Lee assumes that the manifest variable covariance matrix Θ and the endogenous latent variable covariance matrix Ψ are diagonal. This assumption of diagonal covariance matrices is restrictive in some applications. For example, researchers often wish to fit models with correlated residuals, and it has been argued that such models are both necessary and under-utilized (Cole, Ciesla, and Steiger 2007). Correlated residuals pose a difficult problem for MCMC methods because they often result in covariance matrices with some (but not all) off-diagonal entries equal to zero. In this situation, we cannot assign an inverse Wishart prior to the full covariance matrix because this does not allow us to fix some off-diagonal entries to zero. Further, if we assign a prior to each free entry of the covariance matrix and sample the parameters separately, we can often obtain covariance matrices that are not positive definite. This can lead to numerical problems, resulting in an inability to carry out model estimation. Finally, it is possible to specify an equivalent model that possesses the required, diagonal matrices. However, the setting of prior distributions can be unclear here: if the analyst specifies prior distributions for her model of interest, then it may be cumbersome to translate these into prior distributions for the equivalent model.
To address the issue of non-diagonal covariance matrices, Muthén and Asparouhov (2012) implemented a random walk method that is based on work by Chib and Greenberg (1998). This method samples free parameters of the covariance matrix via Metropolis-Hastings steps. While the implementation is fast and efficient, it does not allow for some types of equality constraints because parameters are updated in blocks (either all parameters in a block must be constrained to equality, or no parameters in a block can be constrained). Further, the method is unreliable for models involving many latent variables. Consequently, Muthén and Asparouhov (2012) implemented three other MCMC methods that are suited to different types of models.
In our initial work on package blavaan, we developed methodology for fitting models with free residual covariance parameters. We sought methodology that (i) would work in JAGS, (ii) was reasonably fast and efficient (by JAGS standards), and (iii) allowed for satisfactory specification of prior distributions. In the next section, we describe the resulting methodology. The methodology can often be used to reliably estimate posterior means to the first decimal place on the order of minutes (often, less than two minutes on desktop computers, but it also depends on model complexity). This will never beat compiled code, but it is typically better than alternative JAGS parameterizations that rely on multivariate distributions.

Parameter expansion
General parameter expansion approaches to Bayesian inference are described by Gelman (2004Gelman ( , 2006, with applications to factor analysis models detailed by Ghosh and Dunson (2009). Our approach here is related to that of Palomo et al. (2007), who employ phantom latent variables (they use the term pseudo-latent variable) to simplify the estimation of models with non-diagonal Θ matrices. This results in a working model that is estimated, with the working model parameters being transformed to the inferential model parameters from Equations (1) and (2). The use of phantom latent variables in SEM has long been known (e.g., Rindskopf 1984), though the Bayesian approach involves new issues associated with prior distribution specification. This is further described below.

Overview
Assuming v nonzero entries in the lower triangle of Θ (i.e., v residual covariances), Palomo et al. (2007) take the inferential model residual vector and reparameterize it as where Λ D is a p × v matrix with many zero entries, D is a v × 1 vector of phantom latent variables, and * is a p × 1 residual vector. This approach is useful because, by carefully choosing the nonzero entries of Λ D , both Ψ D and Θ * are diagonal. This allows us to employ an approach related to Lee (2007), avoiding high-dimensional normal distributions in favor of univariate normals.
Under this working model parameterization, the inferential model covariance matrix Θ can be re-obtained via We can also handle covariances between latent variables in the same way, as necessary. Assuming m latent variables with w covariances, we have where the original covariance matrix Ψ is re-obtained via: This approach to handling latent variable covariances can lead to slow convergence in models with many (say, > 5) correlated latent variables. In these cases, we have found it better (in JAGS) to sample the latent variables from multivariate normal distributions. Package blavaan attempts to use multivariate normal distributions in situations where it can, reverting to the above reparameterization only when necessary.
To choose nonzero entries of Λ D (with the same method applying to nonzero entries of B E ), we define two v × 1 vectors r and c. These vectors contain the respective row numbers and column numbers of the nonzero, lower-triangular entries of Θ. For j = 1, . . . , v, the nonzero entries of Λ D then occur in column j, rows r j and c j . Palomo et al. (2007) set these nonzero entries equal to 1, which can be problematic if the covariance parameters' posteriors are negative or overlap with zero. This is because the only remaining free parameters are variances, which can only be positive. Instead of 1s, we free the parameters in Λ D so that they are functions of inferential model parameters. This is further described in the next section. First, however, we give an example of the approach.

Example
Consider the political democracy example of Bollen (1989), which includes eleven observed variables that are hypothesized to arise from three latent variables. The inferential structural equation model affiliated with these data appears as a path diagram (with variance parameters omitted) at the top of Figure 1. This model can be written in matrix form as y 1 y 2 y 3 y 4 y 5 y 6 y 7 y 8 y 9 y 10 y 11 e 1 e 2 e 3 e 4 e 5 e 6 e 7 e 8 e 9 e 10 e 11 Additionally, the latent residual covariance matrix Ψ is diagonal and the observed residual covariance matrix is θ 7,11 θ 9,11 θ 11,11 where the blank entries all equal zero. The six unique, off-diagonal entries of Θ make the model difficult to estimate via Bayesian methods.
To ease model estimation, we follow Palomo et al. (2007) in reparameterizing the model to the working model displayed at the bottom of Figure 1. The working model is now y 1 y 2 y 3 y 4 y 5 y 6 y 7 y 8 y 9 y 10 y 11 The covariance matrix associated with * , Θ * , is now diagonal, which makes the model easier to estimate via MCMC. The latent variable residual, ζ, is maintained as before because its covariance matrix was already diagonal. In general, we only reparameterize residual covariance matrices that are neither diagonal nor unconstrained.
As mentioned earlier, the difference between our approach and that of Palomo et al. (2007) is that we estimate the loadings with D subscripts, whereas Palomo et al. (2007) fix these loadings to 1. This allows us to obtain posteriors on residual covariances that overlap with zero or that become negative, as necessary. Estimation of these loadings comes with a cost, however, in that the prior distributions of the inferential model do not immediately translate into prior distributions of the working model. This problem and a solution are further described in the next section.

Priors on covariances
Due to the reparameterization of covariances, prior distributions on the working model parameters require care in their specification. Every inferential model covariance parameter is potentially the product of three parameters: two representing paths from a phantom latent variable to the variable of interest, and one precision (variance) of the phantom latent variable. These three parameters also impact the inferential model variance parameters. We carefully restrict these parameters to arrive at prior distributions that are both meaningful and flexible.
We place univariate prior distributions on model variance and correlation parameters, with the distributions being based on the inverse Wishart or on other prior distributions described in the literature. This is highly related to the approach of Barnard et al. (2000), who separately specify prior distributions on the variance and correlation parameters comprising a covariance matrix. Our approach is also related to Muthén and Asparouhov (2012), who focus on marginal distributions associated with the inverse Wishart. A notable difference is that we always use proper prior distributions, whereas Muthén and Asparouhov (2012) often use improper prior distributions.
To illustrate the approach, we refer to Figure 2. This figure uses path diagrams to display two observed variables with correlated residuals (left panel) along with their parameter-expanded representation in the working model (right panel). Our goal is to specify a sensible prior distribution for the inferential model's variance/covariance parameters in the left panel, using the parameters in the right panel. The covariance of interest, θ 12 , has been converted to a normal phantom latent variable with variance ψ D and two directed paths, λ 1 and λ 2 . To implement an approach related to that of Barnard et al. (2000) in JAGS, we set where ρ 12 is the correlation associated with the covariance θ 12 , and sign(a) equals 1 if a > 0 and −1 otherwise.
Using this combination of parameters, we need only set prior distributions on the inferential model parameters θ 11 , θ 22 , and ρ 12 , with θ 12 being obtained from these parameters in the usual way. If we want our priors to be analogous to an inverse Wishart prior with d degrees of freedom and diagonal scale matrix S, we can set univariate priors of where Beta (−1,1) is a beta distribution with support on (−1, 1) instead of the usual (0, 1) and p is the dimension of S (typically, the number of observed variables). These priors are related to those used by Muthén and Asparouhov (2012), based on results from Barnard et al. (2000). They are also the default priors for variance/correlation parameters in blavaan, with d = (p + 1) and S = I. We refer to this parameterization option as "srs", reflecting the fact that we are dissecting the covariance parameters into standard deviation and correlation parameters. Barnard et al. (2000) avoid the inverse gamma priors on the variance parameters, instead opting for log-normal distributions on the standard deviations that they judged to be more numerically stable. Package blavaan allows for custom priors, so that the user can employ the log-normal or others on precisions, variances, or standard deviations (this is illustrated later). This approach is generally advantageous because it allows for flexible specification of prior distributions on covariance matrices, including those with fixed zeros or those where we have different amounts of prior information on certain variance/correlation parameters within the matrix.
While we view the "srs" approach as optimal for dealing with covariance parameters here, there exist similar alternative approaches in JAGS. For example, we can treat the phantom latent variables similarly to the other latent variables in the model, assigning prior distributions in the same manner. This alternative approach, which we call "fa" (because the priors resemble traditional factor analysis priors), involves the following default prior distributions Working model Figure 2: Example of phantom latent variable approach to covariances. on the working model from Figure 2: with the inferential model parameters being obtained in the usual way: The main disadvantage of the "fa" approach is that the implied prior distributions on the inferential model parameters are not of a common form. Thus, it is difficult to introduce informative prior distributions for the inferential model variance/covariance parameters. For example, the prior on the inferential covariance (θ 12 ) is the product of two normal prior distributions and an inverse gamma, which can be surprisingly informative. To avoid confusion here, we do not allow the user to directly modify the priors on the working parameters λ 1 , λ 2 , and ψ D under the "fa" approach. The priors chosen above are approximately noninformative for most applications, and the user can further modify the exported JAGS code if he/she desires. Along with the prior distribution issue, the working model parameters are not identified by the likelihood because each inferential covariance parameter is the product of the three working parameters. This can complicate Laplace approximation of the marginal likelihood (which, as described below, is related to the Bayes factor), because the approximation works best when the posterior distributions are unimodal.
In blavaan, the user is free to specify "srs" or "fa" priors for all covariance parameters in the model. In the example section, we present a simple comparison of these two options' relative speeds and efficiencies. Beyond these options, the package attempts to identify when it can sample latent variables directly from a multivariate normal distribution, which can improve sampling efficiency. This most commonly happens when we have many exogenous latent variables that all covary with one another. In this case, we place an inverse Wishart prior distribution on the latent variable covariance matrix.
While we prefer the "srs" approach to covariances between observed variables, the "fa" approach may be useful for complex models that run slowly. In these cases, it could be useful to sacrifice precise prior specification in favor of speed.

Model fit and comparison
Package blavaan supplies a variety of statistics for model evaluation and comparison, available via the fitMeasures() function. This function is illustrated later in the examples, and specific statistics are described in Appendix A. For model evaluation, blavaan supplies posterior predictive checks of the model's log-likelihood (e.g., Gelman, Carlin, Stern, and Rubin 2004). For model comparision, it supplies a variety of statistics including the Deviance Information Criterion (DIC; Spiegelhalter, Best, Carlin, and Linde 2002), the (Laplace-approximated) log-Bayes factor, the Widely Applicable Information Criterion (WAIC; Watanabe 2010), and the leave-one-out cross-validation statistic (LOO; e.g., Gelfand 1996). The latter two statistics are computed by R package loo (Vehtari, Gelman, and Gabry 2015b), using output from blavaan.
Calculation of the information criteria is somewhat complicated by the fact that, during JAGS model estimation, we condition on the latent variables η i , i = 1, . . . , n. That is, our model utilizes likelihoods of the form L(ϑ, η i |y i ), where ϑ is the vector of model parameters and η i is a vector of latent variables associated with individual i. Because the latent variables are random effects, we must integrate them out to obtain L(ϑ|Y ), the likelihood by which model assessment statistics are typically computed. This is easy to do for the models considered here, because the integrated likelihood continues to be multivariate normal. However, we cannot rely on JAGS to automatically calculate the correct likelihood, so that blavaan calculates the likelihood separately after parameters have been sampled in JAGS. Now that we have provided background information on the models implemented in blavaan, the next section further describes how the package works. We then provide a series of examples.

Overview of blavaan
Readers familiar with lavaan will also be familiar with blavaan. The main functions are the same as the main lavaan functions, except that they start with the letter 'b'. For example, confirmatory factor analysis models are estimated via bcfa() and structural equation models are estimated via bsem(). These two functions call the more general blavaan() function with a specific arrangement of arguments. The blavaan model specification syntax is nearly the same as the lavaan model specification syntax, and many functions are used by both packages.
As compared to lavaan, there are a small number of new features in blavaan that are specific to Bayesian modeling: prior distribution specification, export/use of JAGS syntax, convergence diagnostics, and specification of initial values. We discuss these topics in the context of a simple, one-factor confirmatory model. The model uses the Holzinger and Swineford (1939) data included with lavaan, with code for model specification and estimation appearing below.

Prior distribution specification
Package blavaan defines a default prior distribution for each type of model parameter (see Equations (1) and (2)). These default priors can be viewed via which includes a separate default prior for eight types of parameters. Default prior distributions are placed on precisions instead of variances, so that the letter i in itheta and ipsi stands for "inverse." Most of these parameters correspond to notation from Equations (1) to (4), with the exception of rho and ibpsi. The rho prior is used only for the "srs" option, as it is for correlation parameters associated with covariances in Θ or Ψ. The ibpsi prior is used for the covariance matrix of latent variables that are all allowed to covary with one another (commonly, blocks of exogenous latent variables). This block prior can improve sampling efficiency and reduce autocorrelation.
For most parameters, the user is free to declare prior distributions using any prior distribution that is defined in JAGS. Changes to default priors can be supplied with the dp argument. Further, the modifiers [sd] and [var] can be used to put a prior on a standard deviation or variance parameter, instead of the corresponding precision parameter. For example, > fit <-bcfa(model, data = HolzingerSwineford1939, + dp = dpriors(nu = "dnorm(4,1)", itheta = "dunif(0,20)[sd]")) sets the default priors for (i) the manifest variables' intercepts to normal with mean 4 and precision 1, and (ii) the manifest variables' standard deviations to uniform with bounds of 0 and 20.
Priors associated with the rho and ibpsi parameters are less amenable to change. The default rho prior is a beta(1,1) distribution with support on (−1, 1); beta distributions with other parameter values could be used, but this distribution must be beta for now. The default prior on blocks of precision parameters is the Wishart with identity scale matrix and degrees of freedom equal to the dimension of the scale matrix plus one. The user cannot easily make changes to this prior distribution via the dp argument. Changes can be made, however, by exporting the JAGS syntax and manually editing it. This is further described in the next section.
In addition to priors on classes of model parameters, the user may wish to set the prior of a specific model parameter. This is accomplished by using the prior() modifier within the model specification. For example, > model <-' visual =~x1 + prior("dnorm(1,1)")*x2 + x3 ' sets a specific prior distribution for the loading parameter from the visual factor to x2. Priors set in this manner (via the model syntax) take precedent over the default priors set via the dp argument. Additionally, if one specifies a covariance parameter in the model syntax, the prior() modifier is for the correlation associated with the covariance, as opposed to the covariance itself. The prior should be a distribution with support on (0, 1) (typically the beta distribution), which is automatically converted to an analogous distribution that has support on (−1, 1).

JAGS syntax
Users may find that some desired features are not currently implemented in blavaan; such features may be related to, e.g., the use of specific types of priors or the handling of discrete variables. In these situations, we allow the user to export the JAGS syntax so that they can implement new features themselves. We believe this will be especially useful for researchers who wish to develop new types of models: these researchers can specify a traditional model in blavaan that is related to the one that they want to implement, and they can then obtain the JAGS syntax and data for the traditional model. This JAGS syntax should ease implementation of the novel model.
To export the JAGS syntax, users can employ the jagfile argument that was used in the bcfa command above. When jagfile is set to TRUE, the syntax will be written to the lavExport folder; if a character string is supplied, blavaan will use the character string as the folder name.
When the syntax is exported, two files are written within the folder. The first, sem.jag, contains the model specification in JAGS syntax. The second, semjags.rda, is a list named jagtrans that contains the data and initial values in JAGS format, as well as labels associated with model parameters. These pieces can be used to run the model "manually" via, e.g., package rjags (Plummer 2015a) or runjags (Denwood 2016). The example below illustrates how to load the JAGS data (along with initial values and parameter labels) and estimate the exported model via runjags.
> load("lavExport/semjags.rda") > fit <-run.jags("lavExport/sem.jag", monitor = jagtrans$coefvec$jlabel, + data = jagtrans$data, inits = jagtrans$inits) If the user modifies the JAGS file, then the monitor and inits arguments might require modification (or, for automatic initial values from JAGS, the user could omit the inits argument). The data should not require modification unless the user adds or removes observed variables from the model.
We also note that the dimensions of the parameter vectors/matrices in the JAGS syntax generally match the dimensions that we would expect from Equations (1) and (2), with a third dimension being added for multiple group models. These matrix dimensions (and the specific nonzero entries within each matrix) are obtained directly from lavaan.

Convergence diagnostics
Package blavaan offers two methods for monitoring chain convergence, specified as the "auto" or "manual" options to the convergence argument. Regardless of which option is used, the default number of chains is three and can be modified via the n.chains argument.
Under "auto" convergence, chains are sampled for an initial period in order to (i) achieve convergence as determined by the Potential Scale Reduction Factor (PSRF; Gelman and Rubin 1992) and (ii) determine the number of samples necessary to obtain precise posterior estimates. Following this initial period, the chains are further sampled for the determined number of iterations. This automatic assessment comes directly from the autorun.jags() function of package runjags (Denwood 2016); see there for further detail.
Under "manual" convergence, the user can specify the desired number of adaptation, burnin, and sample iterations via arguments of the same name (with defaults of 1,000 adaptation iterations, 4,000 burnin iterations, and 10,000 sample iterations). Adaptation iterations are those where samplers can modify themselves to increase sampling efficiency, whereas burnin iterations are discarded samples obtained from "fixed" samplers (e.g., Plummer 2015b). A warning is issued if any PSRF is greater than 1.2, though some users may desire a more stringent criterion that is closer to 1.0. Access to the parameters' PSRF values is obtained via > blavInspect(fit, "psrf") Additionally, there is a plot method that includes the familiar time series plots, autocorrelation plots, and more. For example, autocorrelation plots for the first four free parameters are obtained via > plot(fit, 1:4, "autocorr") where parameter numbers correspond to the ordering of coef(fit). This plot functionality comes from package runjags, with plotting options being found in the help file for plot.runjags().

Initial values
In many Bayesian SEMs, poorly-chosen initial values can lead to extremely slow convergence. Package blavaan provides multiple options for supplying initial values that vary from chain to chain and/or that improve chain convergence. These are obtained via the inits argument. Under option "prior" (default), starting values are random draws from each parameter's prior distribution. However, to aid in convergence, loading and regression parameters are all required to be positive and close to 1, and correlation parameters are required to be close to zero. Under option "jags", starting values are automatically set by JAGS.
Additional options for starting values are obtained as byproducts of lavaan. For example, the options "Mplus" and "simple" are analogous to the options for the lavaan start argument, where the same initial values are used for all chains. Individual initial values can also be set via the start() modifier to the lavaan syntax, and these initial values will be used for all chains.
Finally, it is also possible to supply a full set of user-defined starting values for each chain. This is most easily accomplished by estimating the model once in order to obtain the list of initial values. The initial values of a fitted model can be obtained via blavInspect(); i.e.,

> myinits <-blavInspect(fit, "inits")
The object myinits is a list whose length equals n.chains, with each entry being another list that contains initial values for the model's parameter vector. It may be helpful to call blavInspect() with the "jagnames" argument, in order to examine correspondence between blavaan parameter names and JAGS parameter names. Finally, after editing, a call to blavaan with the argument inits = myinits will estimate a model using the custom starting values. Now that we have seen some features that are novel to blavaan, we will further illustrate the package by example.

Applications
In the applications below, we first illustrate some general features involving prior distribution specification and model fit measures. We then illustrate the study of measurement invariance via Bayesian models, which involves multiple across-group parameter constraints. Finally, we provide some advanced examples involving the direct modification of exported JAGS code and the use of informative prior distributions.
> plot(fit, 1:4, "trace") As described earlier in the paper, users can also control the handling of residual covariances via the cp argument. For the particular model considered here, we compared the sampling speed and efficiency (excluding posterior predictive checks) of the "srs" and "fa" options using a four-core Dell Precision laptop running Ubuntu linux. We observed that the "srs" option took 76 seconds and the "fa" option took 65 seconds, illustrating the speed advantage of the "fa" approach. However, the ratio of effective sample sizes for the two approaches averaged 1.44 in favor of the "srs" approach (where the average is taken across parameters), implying that "srs" leads to more efficient sampling. In our experience, these results can vary depending on the specific model of interest, but it generally appears that "fa" is somewhat faster and "srs" somewhat more efficient.
Once we have sampled for the desired number of iterations, we can make use of lavaan functions in order to summarize and study the fitted model. Users may primarily be interested in summary(), which organizes results in a manner similar to the classical lavaan models. While the results look similar to those that one would see in lavaan, there are a few differences that require explanation. First, the top of the output includes two model evaluation measures: a Laplace approximation of the marginal log-likelihood and the posterior predictive p-value (see Appendix A). Second, the "Parameter Estimates" section contains many new columns. These are the posterior mean (Post.Mean), the posterior standard deviation (Post.SD), a 95% highest posterior density interval (HPD.025 and HPD.975), the potential scale reduction factor for assessing chain convergence (PSRF; Gelman and Rubin 1992), and the prior distribution associated with each model parameter (Prior). Users can optionally obtain posterior medians, posterior modes, and effective sample sizes (number of posterior samples drawn, adjusted for autocorrelation). These can be obtained by supplying the logical arguments postmedian, postmode, and neff to summary(). where, as previously mentioned, the WAIC and LOO statistics are computed with the help of package loo. Other lavaan functions, including parameterEstimates() and parTable(), similarly work as expected.

Measurement invariance
Verhagen and Fox (2013; see also Verhagen, Levy, Millsap, and Fox 2016) describe the use of Bayesian methods for studying the measurement invariance of sets of items or scales. Briefly, measurement invariance means that the items or scales measure diverse groups of individuals in a fair manner. For example, two individuals with the same level of mathematics ability should receive the same score on a mathematics test (within random noise), regardless of the individuals' backgrounds and demographic variables. In this section, we illustrate a Bayesian measurement invariance approach using the Holzinger and Swineford (1939) data. This section also illustrates how we can estimate a model in blavaan using pieces of a fitted lavaan object. This may be of interest to advanced users who would prefer to write/edit a lavaan "parameter table" instead of using model specification syntax. Additionally, the use of a lavaan object provides a path from Mplus to blavaan, via the mplus2lavaan() function from lavaan.
Following model estimation, we can immediately compare the models via fitMeasures().
For example, focusing on the first two models, we obtain As judged by the posterior predictive p-value (ppp), neither of the two models provide a good fit to the data. In practice, we might stop our measurement invariance study there. However, for the purpose of demonstration, we can also examine the information criteria. Based on these, we would prefer the second model (due to the smaller DIC, WAIC, and LOOIC values).

Informative prior distributions
We have seen that blavaan allows users to specify prior distributions for individual model parameters and for classes of model parameters. In this section, we illustrate how informative prior distributions can be set for factor analysis parameters in the context of a specific dataset, where priors are set based on the parameters' interpretations.
We consider a one-factor, three-indicator model fit to a subset of data from a stereotype threat study (Wicherts, Dolan, and Hessen 2005). The data are available in R package psychotools (Zeileis, Strobl, Wickelmaier, Abou El-Komboz, and Kopf 2015), and the model specifies a general "ability" factor underlying three academic tests (of verbal ability, numerical ability, and abstract reasoning). The test scores are sums of the number of items answered correctly, so they are bounded from below at 0 and from above at 16 (verbal), 14 (numerical), and 18 (abstract reasoning). While these bounds technically violate the normality assumption of the factor analysis model, it is customary to apply this model to test scores (especially in the current situation, where the test scores are unimodal with modes in the middle of the score ranges). A factor analysis model with default (non-informative) prior distributions and automatic convergence assessment can be fit to the data via > library("psychotools") > data("StereotypeThreat") > st <-subset(StereotypeThreat, ethnicity == "majority") > model <-'ability =~abstract + verbal + numerical' > bfit <-bcfa(model, data = st, convergence = "auto") Instead of using the default prior distributions, we can tailor prior distributions to this particular dataset by simultaneously considering the data and each model parameter's interpretation. First, because the test scores are all bounded, we know the maximum possible standard deviation on each test (obtained if half the participants scored a 0 and half score the maximum). These are 8 for the verbal test, 7 for the numerical test, and 9 for the abstract reasoning test. Thus, we place uniform priors on the manifest variables' standard deviations with lower bounds at 0 and upper bounds at these maxima: > model <-c(model, + 'abstract~~prior ("dunif(0,9)[sd]") * abstract + verbal~~prior ("dunif(0,8)[sd]") * verbal + numerical~~prior("dunif(0,7)[sd]") * numerical') Next, we consider the factor variance. Because the loading associated with the abstract reasoning test will be fixed at 1, the factor variance can be interpreted as the variability in the abstract reasoning test that can be attributed to the ability factor. Because the three tests were chosen to be highly correlated (and because the total standard deviation of abstract reasoning is 9 or less), we use a uniform prior with an upper bound of 6 here. Thus, even in the case where we observe the maximum standard deviation of 9, the factor could still account for 2/3 of the variability in the ability test.
> model <-c(model, + 'ability~~prior("dunif(0,4.5)[sd]") * ability') Next, we consider the manifest variables' intercept parameters, which reflect average scores on each of the three tests. These tests are known historically to result in average scores near the middle of the range (Wicherts et al. 2005), so our prior distributions are truncated normals centered at the middle of each test's range (with the truncation points at the minimum and maximum possible scores on each test): > model <-c(model, + 'abstract~prior("dnorm(9,.25)T(0,18)") * 1 + verbal~prior ("dnorm(8,.25)T(0,16)") * 1 + numerical~prior("dnorm(7,.25)T(0,14)") * 1') Finally, we consider the factor loadings. We have no knowledge that the ability factor will influence one test more than others, so we might expect the two free loadings to be around 1 (because the fixed loading equals 1). We also expect positive loadings, because ability should positively influence all three tests. Thus, we supply uniform(0,3) priors to the dp argument of bcfa(): > bfitInform <-bcfa(model, data = st, + dp = dpriors(lambda = "dunif(0,3)"), + convergence = "auto") The two models' fit measures are generally similar and not shown. We compare the two models' posterior means and standard deviations in Table 1. The differences are unlikely to impact substantive conclusions, but two of them are noteworthy. First, the factor variance ψ 1 is larger under the model with informative priors, likely because the informative prior (uniform(0,4.5) on the standard deviation) placed more density on larger values of the standard deviation. We observe a similar phenomenon with the residual variance of the numerical test (θ 3 ). Second, the posterior means and standard deviations of the loadings (λ 2 and λ 3 ) are somewhat smaller under the informative priors. This is likely related to the larger variance estimates.
This example could be used as the start of a larger analysis of posterior sensitivity to prior distributions. For the purpose of automation, prior distributions could be entered directly into the model's parameter table (obtained via, e.g., parTable(bfit)) and the model subsequently re-estimated, similar to what was done in the measurement invariance example.

Extensions of JAGS syntax
Finally, we provide an example involving use of the JAGS syntax to estimate novel models. Consider the robust factor analysis model of Zhang, Li, and Liu (2014), which essentially replaces the factor analysis model's normal distributions with t distributions. In particular, each factor η arises from a t df (0, 1) distribution, and each residual j arises from a t df (0, ψ j ) distribution. The degrees of freedom, df, is a free parameter and is shared by all the t distributions. To our knowledge, there exists no software to readily estimate this model, and Zhang et al. (2014) implemented their own Gibbs sampler (making use of analytic posterior distributions under conjugate priors). Here, we show how the JAGS syntax from blavaan can be modified to easily estimate this model. For illustration, we again apply a one-factor model to the Wicherts et al. (2005) data from the previous section.
We begin by using blavaan to export JAGS syntax for a simple, one-factor model.  By default, the exported JAGS code is written to the file lavExport/sem.jag. Upon opening that file, we must make two edits (see Figure 4). First, the normal distributions (dnorm()) are replaced with t distributions (dt()). Second, we need to specify a prior distribution for the degrees of freedom parameter. Similar to Zhang et al. (2014), we place a flat prior on the inverse degrees of freedom. The lower bound is nonzero to aid in convergence; this is justified by realizing that, as the degrees of freedom increase, the t distribution becomes the normal distribution. Once we surpass, say, 200 degrees of freedom, we can accept that larger values are practically equivalent to 200.
The modified syntax from Figure 4 could then be estimated manually, making use of the exported data and adding a monitor for the df parameter. For example, if the exported files are saved as lavExport/robustfa.jag and lavExport/robustfa.rda, the model can be estimated via > load("lavExport/robustfa.rda") > fit <-run.jags("lavExport/robustfa.jag", + monitor = c(jagtrans$monitors, "df"), + data = jagtrans$data, inits = jagtrans$inits) While many model extensions will be more complicated than the one considered here, this example illustrates blavaan's potential use for statistical researchers who are developing new models: instead of writing JAGS syntax entirely from scratch, these researchers may use blavaan to obtain syntax for a basic model that is similar to the desired model. In many cases, this will simplify implementation of the desired model.

Conclusion
As described throughout the paper, package blavaan combines many existing tools so that it is easy to estimate multivariate normal SEMs via open-source software. The package can be useful to applied researchers who need an expanded set of prior distributions at their disposal, freeing them from the need to learn the intricacies of an MCMC package. It can also be useful to methodological researchers, freeing them from the need to code their own JAGS models from scratch. The idea of separating model specification from MCMC coding seems generally useful beyond SEM, and others are indeed making progress in this direction. For example, package rstanarm (Gabry and Goodrich 2016) allows users to estimate generalized linear mixed models via Stan, using lme4 (Bates, Mächler, Bolker, and Walker 2015) syntax. The coding of complex models in BUGS, JAGS, or Stan is often tedious even for experienced users, so that progress here could improve most users' workflows.
There are a variety of additional models that we plan to support in future versions of blavaan, including models with latent interactions, with ordinal variables, and with mixture and multilevel components. These models have been addressed in the literature (e.g., Asparouhov and Muthén 2010;Lee 2001, 2012), and JAGS examples for estimating specific instances of these models is available (e.g., Cho, Preacher, and Bottge 2015; Merkle and Wang in press). Work remains, however, to sample from these models efficiently and to specify the models via lavaan's syntax. The sampling efficiency of Stan may be useful for at least some of these models, and we plan to explore this in the future.
In summary, blavaan is currently useful for estimating many types of multivariate normal SEMs, and the JAGS export feature allows researchers to extend the models in any fashion desired. As additional features are added, we hope that the package keeps pace with lavaan as an open set of tools for SEM estimation and study.

Acknowledgments
A. Details on model assessment statistics In the subsections below, we review the model evaluation and comparison metrics that blavaan supplies.

Posterior predictive checks
As a measure of the model's absolute fit, blavaan computes a posterior predictive p-value that compares observed likelihood ratio test (χ 2 ) statistics to likelihood ratio test statistics generated from the model's posterior predictive distribution (also see Muthén and Asparouhov 2012;Scheines et al. 1999). The likelihood ratio test statistic is generally computed via where ϑ sat is a "saturated" parameter vector that perfectly matches the observed mean and covariance matrix.
2. Generate artificial data Y rep from the model, assuming parameter values equal to ϑ s .

Record whether or not LRT
The posterior predictive p-value is then the proportion of the time (out of S draws) that the posterior predictive LRT statistic is larger than the observed LRT statistic. Values closer to .5 indicate a model that fits the observed data, while values closer to 0 indicate the opposite. For use in practice, Muthén and Asparouhov (2012) use a threshold of .05 (analogous to a frequentist α level) and present some evidence that, as compared to the frequentist likelihood ratio test, the posterior predictive p is less sensitive to minor model misspecifications and exhibits better performance at small sample sizes. On the other hand, Hoijtink and van de Schoot (2016) show that the posterior predictive p-value misbehaves in situations involving highly informative prior distributions (as are sometimes used for shrinkage/regularization). They recommend a prior-posterior predictive p-value that may be implemented in future versions of blavaan.
In blavaan, the posterior predictive p-value is currently impractical when there are missing data. This is because it is impractical to compute the saturated log-likelihood in the presence of missing data. We need to use, e.g., the EM algorithm to do this, and this needs to be done (S + 1) times (once for the saturated likelihood of the observed data and once for the saturated likelihood of each of the S artificial datasets). With complete data, we can more simply calculate the saturated log-likelihood using the sample mean vector and covariance matrix. The argument test="none" can be supplied to blavaan in order to bypass these slow calculations when there are missing data, and future versions may sample the missing observations in order to avoid the issue.

DIC
The DIC is given as DIC = −2 log L(θ|Y ) + 2 efp DIC , whereθ is a vector of posterior parameter means (or other measure of central tendency), L(θ|Y ) is the model's likelihood at the posterior means, and efp DIC is the model's effective number of parameters, calculated as efp DIC = 2 log L(θ|Y ) − log L(ϑ|Y ) .
The latter term, log L(ϑ|Y ), is obtained by calculating the log-likelihood for each posterior sample and then averaging.
Many readers will be familiar with the automatic calculation of DIC within programs like BUGS and JAGS. As we described above, the automatic DIC is not what users typically desire because the likelihood conditions on the latent variables η i . This greatly increases the effective number of parameters and may result in poor inferences (for further discussion of this issue, see Millar 2009). As previously mentioned, blavaan avoids the automatic DIC computation in JAGS, calculating its own likelihoods after model parameters have been sampled.

WAIC and LOO
The WAIC and LOO are asymptotically equivalent measures that are advantageous over DIC while also being more difficult to compute than DIC. Many computational difficulties are overcome in developments by Vehtari, Gelman, and Gabry (2015a), with their methodology being implemented in package loo (Vehtari et al. 2015b). As input, loo requires casewise loglikelihoods associated with a set of posterior samples: log L(ϑ s |y i ), s = 1, . . . , S; i = 1, . . . , n.
Like DIC, both of these measures seek to characterize a model's predictive accuracy (generalizability). The definition of WAIC looks similar to that of DIC: where the first term is related to log-likelihoods of observed data and the second term involves an effective number of parameters. Both terms now involve log-likelihoods associated with individual observations, however, which help us estimate the model's expected log pointwise predictive density (a measure of predictive accuracy).
The first term, the log pointwise predictive density of the observed data (lppd), is estimated via where S is the number of posterior draws and f (y i |ϑ s ) is the density of observation i with respect to the parameters sampled at iteration s. The second term, the effective number of parameters, is estimated via var s (log f (y i |ϑ)), where we compute a separate variance for each observation i across the S posterior draws.
The LOO measure estimates the predictive density of each individual observation from a cross-validation standpoint; that is, the predictive density when we hold out one observation at a time and use the remaining observations to update the prior. We can use analytical results to estimate this quantity based on the full-data posterior, removing the need to reestimate the posterior while sequentially holding out each observation. These results lead to the estimate where the w s i are importance sampling weights based on the relative magnitude of individual i's density function across the S posterior samples. Package loo smooths these weights via a generalized Pareto distribution, improving the estimate.
We can also estimate the effective number of parameters under the LOO measure by comparing LOO to the lppd that was used in the WAIC calculation. This gives us efp LOO = lppd + LOO/2, where these terms come from Equations (9) and (10), respectively. Division of the latter term by two (and addition, vs. subtraction) offsets the multiplication by −2 that occurs in (10). For further detail on all these measures, see Vehtari et al. (2015a).

Bayes factor
The Bayes factor between two models is defined as the ratio of the models' marginal likelihoods (e.g., Kass and Raftery 1995). Candidate model 1's marginal likelihood can be written as where ϑ 1 is a vector containing candidate model 1's free parameters (excluding the latent variables η i ) and f () is a probability density function. The marginal likelihood of candidate model 2 may be written in the same manner. The Bayes factor is then with values greater than 1 favoring Model 1.
Because the integral from (11) is generally difficult to calculate, Lewis and Raftery (1997) describe a Laplace-Metropolis estimator of the Bayes factor (also see Raftery 1993). This estimator relies on the Laplace approximation to integrals that involve a natural exponent. Such integrals can be written as exp(h(u))∂u ≈ (2π) Q/2 |H * | 1/2 exp(h(u * )), where h() is a function of the vector u, Q is the length of u, u * is the value of u that maximizes h, and H * is the inverse of the information matrix evaluated at u * . As applied to marginal likelihoods, we take h(ϑ) = log(f (ϑ)f (Y |ϑ)) so that our solution to Equation (11) is where ϑ * is a vector of posterior means (or other posterior estimate of central tendency), obtained from the MCMC output. In package blavaan, we follow Lewis and Raftery (1997) and compute the logarithm of this approximation for numerical stability.
The Bayes factor is sensitive to choice of prior distribution (e.g., Liu and Aitkin 2008;Vanpaemel 2010), so researchers using the Bayes factor are advised to carefully consider their priors. Kass and Raftery (1995) provide popular rules of thumb for interpreting the log-Bayes factor, though the extent to which these rules are meaningful in any specific application is unclear. In general, as the log-Bayes factor increases from 0, we gain increasing support for