A Novel Bayesian Approach for Latent Variable Modeling from Mixed Data with Missing Values
Abstract
We consider the problem of learning parameters of latent variable models from mixed (continuous and ordinal) data with missing values. We propose a novel Bayesian Gaussian copula factor (BGCF) approach that is consistent under certain conditions and that is quite robust to the violations of these conditions. In simulations, BGCF substantially outperforms two stateoftheart alternative approaches. An illustration on the ‘Holzinger & Swineford 1939’ dataset indicates that BGCF is favorable over the socalled robust maximum likelihood (MLR) even if the data match the assumptions of MLR.
lightred1.0 0.7 0.7rgb
Keywords: latent variables; Gaussian copula factor model; parameter learning; mixed data; missing values
1 Introduction
In psychology, social sciences, and many other fields, researchers are usually interested in “latent” variables that cannot be measured directly, e.g., depression, anxiety, or intelligence. To get a grip on these latent concepts, one commonlyused strategy is to construct a measurement model for such a latent variable, in the sense that domain experts design multiple “items” or “questions” that are considered to be indicators of the latent variable. For exploring evidence of construct validity in theorybased instrument construction, confirmatory factor analysis (CFA) has been widely studied (Jöreskog, 1969; Castro et al., 2015; Li, 2016). In CFA, researchers start with several hypothesised latent variable models that are then fitted to the data individually, after which the one that fits the data best is picked to explain the observed phenomenon. In this process, the fundamental task is to learn the parameters of a hypothesised model from observed data, which is the focus of this paper. For convenience, we simply refer to these hypothesised latent variable models as CFA models from now on.
The most common method for parameter estimation in CFA models is maximum likelihood (ML), because of its attractive statistical properties (consistency, asymptotic normality, and efficiency). The ML method, however, relies on the assumption that observed variables follow a multivariate normal distribution (Jöreskog, 1969). When the normality assumption is not deemed empirically tenable, ML may not only reduce the accuracy of parameter estimates, but may also yield misleading conclusions drawn from empirical data (Li, 2016). To this end, a robust version of ML was introduced for CFA models when the normality assumption is slightly or moderately violated (Kaplan, 2008), but still requires the observations to be continuous. In the real world, the indicator data in questionnaires are usually measured on an ordinal scale (resulting in a bunch of ordered categorical variables, or simply ordinal variables) (Poon and Wang, 2012), in which neither normality nor continuity is plausible (Lubke and Muthén, 2004). In such cases, diagonally weighted least squares (DWLS in LISREL; WLSMV or robust WLS in Mplus) has been suggested to be superior to the ML method and is usually considered to be preferable over other methods (Barendse et al., 2015; Li, 2016).
However, there are two major issues that the existing approaches do not consider. One is the mixture of continuous and ordinal data. As we mentioned above ordinal variables are omnipresent in questionnaires, whereas sensor data are usually continuous. Therefore, a more realistic case in real applications is mixed continuous and ordinal data. A second important issue concerns missing values. In practice, all branches of experimental science are plagued by missing values (Little and Rubin, 1987), e.g., failure of sensors, or unwillingness to answer certain questions in a survey. A straightforward idea in this case is to combine missing values techniques with existing parameter estimation approaches, e.g., performing listwisedeletion or pairwisedeletion first on the original data and then applying DWLS to learn parameters of a CFA model. However, such deletion methods are only consistent when the data are missing completely at random (MCAR), which is a rather strong assumption (Rubin, 1976), and cannot transfer the sampling variability incurred by missing values to followup studies. The two modern missing data techniques, maximum likelihood and multiple imputation, are valid under a less restrictive assumption, missing at random (MAR) (Schafer and Graham, 2002), but they require the data to be multivariate normal.
Therefore, there is a strong demand for an approach that is not only valid under MAR but also works for mixed continuous and ordinal data. For this purpose, we propose a novel Bayesian Gaussian copula factor (BGCF) approach, in which a Gibbs sampler is used to draw pseudo Gaussian data in a latent space restricted by the observed data (unrestricted if that value is missing) and draw posterior samples of parameters given the pseudo data, iteratively. We prove that this approach is consistent under MCAR and empirically show that it works quite well under MAR.
The rest of this paper is organized as follows. Section 2 reviews background knowledge and related work. Section 3 gives the definition of a Gaussian copula factor model and presents our novel inference procedure for this model. Section 4 compares our BGCF approach with two alternative approaches on simulated data, and Section 5 gives an illustration on the ‘Holzinger & Swineford 1939’ dataset. Section 6 concludes this paper and provides some discussion.
2 Background
This section reviews basic missingness mechanisms and related work on parameter estimation in CFA models.
2.1 Missingness Mechanism
Following Rubin (1976), let be a data matrix with the rows representing independent samples, and be a matrix of indicators, where if was observed and otherwise. consists of two parts, and , representing observed and missing elements in respectively. When the missingness does not depend on the data, i.e., with denoting unknown parameters, the data are said to be missing completely at random (MCAR), which is a special case of a more realistic assumption called missing at random (MAR). MAR allows the dependency between missingness and observed values, i.e., . For example, all people in a group are required to take a blood pressure test at time point 1, while only those whose values at time point 1 lie in the abnormal range need to take the test at time point 2. This results in some missing values at time point 2 that are MAR.
2.2 Parameter Estimation in CFA Models
When the observations follow a multivariate normal distribution, maximum likelihood (ML) is the mostlyused method. It is equivalent to minimizing the discrepancy function (Jöreskog, 1969):
where is the vector of model parameters, is the modelimplied covariance matrix, is the sample covariance matrix, and is the number of observed variables in the model. When the normality assumption is violated either slightly or moderately, robust ML (MLR) offers an alternative. Here parameter estimates are still obtained using the asymptotically unbiased ML estimator, but standard errors are statistically corrected to enhance the robustness of ML against departures from normality (Kaplan, 2008; Muthén, 2010). Another method for continuous nonnormal data is the socalled asymptotically distribution free method, which is a weighted least squares (WLS) method using the inverse of the asymptotic covariance matrix of the sample variances and covariances as a weight matrix (Browne, 1984).
When the observed data are on ordinal scales, Muthén (1984) proposed a threestage approach. It assumes that a normal latent variable underlies an observed ordinal variable , i.e.,
where denotes the observed values of , are thresholds , and is the number of categories. The thresholds and polychoric correlations are estimated from the bivariate contingency table in the first two stages (Olsson, 1979; Jöreskog, 2005). Parameter estimates and the associated standard errors are then obtained by minimizing the weighted least squares fit function :
where is the vector of model parameters, is the modelimplied vector containing the nonredundant vectorized elements of , is the vector containing the estimated polychoric correlations, and the weight matrix is the asymptotic covariance matrix of the polychoric correlations. A mathematically simple form of the WLS estimator, the unweighted least squares (ULS), arises when the matrix is replaced with the identity matrix . Another variant of WLS is the diagonally weighted least squares (DWLS), in which only the diagonal elements of are used in the fit function (Muthén et al., 1997; Muthén, 2010), i.e.,
where is the diagonal weight matrix. Various recent simulation studies have shown that DWLS is favorable compared to WLS, ULS, as well as the MLbased methods for ordinal data (Barendse et al., 2015; Li, 2016).
3 Method
In this section, we introduce the Gaussian copula factor model and propose a Bayesian inference procedure for this model. Then, we theoretically analyze the identifiability and prove the consistency of our procedure.
3.1 Gaussian Copula Factor Model
Definition 1 (Gaussian Copula Factor Model).
Consider a latent random (factor) vector , a response random vector and an observed random vector , satisfying
(1)  
(2)  
(3) 
with a correlation matrix over factors, a matrix of factor loadings (), residuals with , the standard deviation of , the cumulative distribution function (CDF) of the standard Gaussian, and the pseudoinverse of a CDF . Then this model is called a Gaussian copula factor model.

The model is also defined in Murray et al. (2013), but the authors restrict the factors to be independent of each other while we allow for their interactions. Our model is a combination of a Gaussian factor model (from to ) and a Gaussian copula model (from to ). The first part allows us to model the latent concepts that are measured by multiple indicators, and the second part provides a good way to model diverse types of variables (depending on in Equation 3, can be either continuous or ordinal). Figure 1 shows an example of the model. Note that we allow the special case of a factor having a single indicator, e.g., , because this allows us to incorporate other (explicit) variables (such as age and income) into our model. In this special case, we set and , thus .
In the typical design for questionnaires, one tries to get a grip on a latent concept through a particular set of welldesigned questions (MartínezTorres, 2006; Byrne, 2013), which implies that a factor (latent concept) in our model is connected to multiple indicators (questions) while an indicator is only used to measure a single factor, as shown in Figure 1. This kind of measurement model is called a pure measurement model (Definition 8 in Silva et al. (2006)). Throughout this paper, we assume that all measurement models are pure, which indicates that there is only a single nonzero entry in each row of the factor loadings matrix . This inductive bias about the sparsity pattern of is fully motivated by the typical design of a measurement model.
In what follows, we transform the Gaussian copula factor model into an equivalent model that is used for inference in the next subsection. We consider an integrated dimensional random vector , which is still multivariate Gaussian, and obtain its covariance matrix
(4) 
and precision matrix
(5) 
Since is diagonal and only has one nonzero entry per row, contains many intrinsic zeros. The sparsity pattern of such can be represented by an undirected graph , where whenever by construction. Then, a Gaussian copula factor model can be transformed into an equivalent model controlled by a single precision matrix , which in turn is constrained by , i.e., .
Definition 2 (Wishart Distribution).
Given an undirected graph , a zeroconstrained random matrix has a Wishart distribution, if its density function is
with the space of symmetric positive definite matrices with offdiagonal elements whenever , the number of degrees of freedom, a scale matrix, the normalizing constant, and the indicator function (Roverato, 2002).
The Wishart distribution is the conjugate prior of precision matrices that are constrained by a graph (Roverato, 2002). That is, given the Wishart prior, i.e., and data drawn from , the posterior for is another Wishart distribution:
When the graph is fully connected, the Wishart distribution reduces to a Wishart distribution (Murphy, 2007). Placing a Wishart prior on is equivalent to placing an inverseWishart on , a product of multivariate normals on , and an inversegamma on the diagonal elements of . With a diagonal scale matrix and the number of degrees of freedom equal to the number of factors plus one, the implied marginal densities between any pair of factors are uniformly distributed between (Barnard et al., 2000).
3.2 Inference for Gaussian Copula Factor Model
We first introduce the inference procedure for complete mixed data and incomplete Gaussian data respectively, based on which the procedure for mixed data with missing values is then derived. From this point on, we use to denote the correlation matrix over the response vector .
Mixed Data without Missing Values
For a Gaussian copula model, Hoff (2007) proposed a likelihood that only concerns the ranks among observations, which is derived as follows. Since the transformation is nondecreasing, observing implies a partial ordering on , i.e., lies in the space restricted by :
Therefore, observing suggests that must be in
Taking the occurrence of this event as the data, one can compute the following likelihood Hoff (2007)
Following the same argumentation, the likelihood in our Gaussian copula factor model reads
which is independent of the margins .
For the Gaussian copula factor model, inference for the precision matrix of the vector can now proceed via construction of a Markov chain having its stationary distribution equal to , where we ignore the values for and in our samples. The prior graph is uniquely determined by the sparsity pattern of the loading matrix and the residual matrix (see Equation 5), which in turn is uniquely decided by the pure measurement models. The Markov chain can be constructed by iterating the following three steps:

Sample : ;
Since each coordinate directly depends on only one factor, i.e., such that , we can sample each of them independently through . 
Sample : ;

Sample : .
Gaussian Data with Missing Values
Suppose that we have Gaussian data consisting of two parts, and , denoting observed and missing values in respectively. The inference for the correlation matrix of in this case can be done via the socalled data augmentation technique that is also a Markov chain Monte Carlo procedure and has been proven to be consistent under MAR (Schafer, 1997). This approach iterates the following two steps to impute missing values (Step 1) and draw correlation matrix samples from the posterior (Step 2):

;

.
Mixed Data with Missing Values
For the most general case of mixed data with missing values, we combine the procedures of Sections 3.2.1 and 3.2.2 into the following fourstep inference procedure:

;

;

;

.
A Gibbs sampler that achieves this Markov chain is summarized in Algorithm 1 and implemented in R.
By iterating the steps in Algorithm 1, we can draw correlation matrix samples over the integrated random vector , denoted by . The mean over all the samples is a natural estimate of the true , i.e.,
(6) 
Based on Equations (4) and (6), we obtain estimates of the parameters of interests:
(7)  
We refer to this procedure as a Bayesian Gaussian copula factor approach (BGCF).
3.3 Theoretical Analysis
Identifiability of Without additional constraints, is nonidentifiable (Anderson and Rubin, 1956). More precisely, given a decomposable matrix , we can always replace with and with to obtain an equivalent decomposition , where is a invertible matrix. Since only has one nonzero entry per row in our model, can only be diagonal to ensure that has the same sparsity pattern as (see Lemma 1 in Appendix). Thus, from the same , we get a class of solutions for , i.e., , where can be any invertible diagonal matrix. In order to get a unique solution for , we impose two sufficient identifying conditions: 1) restrict to be a correlation matrix; 2) force the first nonzero entry in each column of to be positive. See Lemma 2 in Appendix for proof. Condition 1 is implemented via line 31 in Algorithm 1. As for the second condition, we force the covariance between a factor and its first indicator to be positive (line 27), which is equivalent to Condition 2. Note that these conditions are not unique; one could choose one’s favorite conditions to identify , e.g., setting the first loading to 1 for each factor. The reason for our choice of conditions is to keep it consistent with our model definition where is a correlation matrix.
Identifiability of and Under the two conditions for identifying , factor loadings and residual variances are also identified except for the case in which there exists one factor that is independent of all the others and this factor only has two indicators. For such a factor, we have 4 free parameters (2 loadings, 2 residuals) while we only have 3 available equations (2 variances, 1 covariance), which yields an underdetermined system. See Lemmas 3 and 4 in Appendix for detailed analysis. Once this happens, one could put additional constraints to guarantee a unique solution, e.g., by setting the variance of the first residual to zero. However, we would recommend to leave such an independent factor out (especially in association analysis) or study it separately from the other factors.
Under sufficient conditions for identifying , , and , our BGCF approach is consistent even with MCAR missing values. This is shown in Theorem 1, whose proof is provided in Appendix.
Theorem 1 (Consistency of the BGCF Approach).
Let be independent observations drawn from a Gaussian copula factor model. If is complete (no missing data) or contains missing values that are missing completely at random, then
where , , and are parameters learned by BGCF, while , , and are the true ones.
4 Simulation Study
In this section, we compare our BGCF approach with alternative approaches via simulations.
4.1 Setup
Model specification Following typical simulation studies on CFA models in the literature (YangWallentin et al., 2010; Li, 2016), we consider a correlated 4factor model in our study. Each factor is measured by 4 indicators, since Marsh et al. (1998) concluded that the accuracy of parameter estimates appeared to be optimal when the number of indicators per factor was four and marginally improved as the number increased. The interfactor correlations (offdiagonal elements of the correlation matrix over factors) are randomly drawn from [0.2, 0.4], which is considered a reasonable and empirical range in the applied literature (Li, 2016). For the ease of reproducibility, we construct our as follows. set.seed(12345) C < matrix(runif(4^2, 0.2, 0.4), ncol=4) C < (C*lower.tri(C)) + t(C*lower.tri(C)) diag(C) < 1 In the majority of empirical research and simulation studies (DiStefano, 2002), reported standardized factor loadings range from 0.4 to 0.9. For facilitating interpretability and again reproducibility, each factor loading is set to 0.7. Each corresponding residual variance is then automatically set to 0.51 under a standardized solution in the population model, as done in (Li, 2016).
Data generation Given the specified model, one can generate data in the response space (the in Definition 1) via Equations (1) and (2). When the observed data (the in Definition 1) are ordinal, we discretize the corresponding margins into the desired number of categories. When the observed data are nonparanormal, we set the in Equation (3) to the CDF of a distribution with degrees of freedom df. The reason for choosing a distribution is that we can easily use df to control the extent of nonnormality: a higher df implies a distribution closer to a Gaussian. To fill in a certain percentage of missing values (we only consider MAR), we follow the procedure in Kolar and Xing (2012), i.e., for , : is missing if .
Evaluation metrics We use average relative bias (ARB) and root mean squared error (RMSE) to examine the parameter estimates, which are defined as
where and represent the estimated and true values respectively. An ARB value less than 5% is interpreted as a trivial bias, between 5% and 10% as a moderate bias, and greater than 10% as a substantial bias (Curran et al., 1996). Note that ARB describes an overall picture of average bias, that is, summing up bias in a positive and a negative direction together. A smaller absolute value of ARB indicates better performance on average.
4.2 Ordinal Data without Missing Values
In this subsection, we consider ordinal complete data since this matches the assumptions of the diagonally weighted least squares (DWLS) method, in which we set the number of ordinal categories to be 4. We also incorporate the robust maximum likelihood (MLR) as an alternative approach, which was shown to be empirically tenable when the number of categories is more than 5 (Rhemtulla et al., 2012; Li, 2016). See Section 2 for details of the two approaches.
Before conducting comparisons, we first check the convergence property of the Gibbs sampler used in our BGCF approach. Figure 2 shows the RMSE of estimated interfactor correlations (left panel) and factor loadings (right panel) over 100 iterations for a randomlydrawn sample with sample size . We see quite a good convergence of the Gibbs sampler, in which the burnin period is only around 10. More experiments done for different numbers of categories and different random samples show that the burnin is less than 20 on the whole across various conditions.
Now we evaluate the three involved approaches. Figure 3 shows the performance of BGCF, DWLS, and MLR over different sample sizes , providing the mean of ARB (left panel) and the mean of RMSE with 95% confidence interval (right panel) over 100 experiments. From Figure (a)a, interfactor correlations are, on average, trivially biased (within two dashed lines) for all the three methods that in turn give indistinguishable RMSE regardless of sample sizes. From Figure (b)b, MLR moderately underestimates the factor loadings, and performs worse than DWLS w.r.t. RMSE especially for a larger sample size, which confirms the conclusion in previous studies (Barendse et al., 2015; Li, 2016). Most importantly, our BGCF approach outperforms DWLS in learning factor loadings especially for small sample sizes, even if the experimental conditions entirely match the assumptions of DWLS.
4.3 Mixed Data with Missing Values
In this subsection, we consider mixed nonparanormal and ordinal data with missing values, since some latent variables in realworld applications are measured by sensors that usually produce continuous but not necessarily Gaussian data. The 8 indicators of the first 2 factors (4 per factor) are transformed into a distribution with , which yields a slightlynonnormal distribution (skewness is 1, excess kurtosis is 1.5) (Li, 2016). The 8 indicators of the last 2 factors are discretized into ordinal with 4 categories.
One alternative approach in such cases is DWLS with pairwisedeletion (PD), in which heterogeneous correlations (Pearson correlations between numeric variables, polyserial correlations between numeric and ordinal variables, and polychoric correlations between ordinal variables) are first computed based on pairwise complete observations, and then DWLS is used to estimate model parameters. A second alternative concerns the full information maximum likelihood (FIML) (Arbuckle, 1996; Rosseel, 2012), which first applies an EM algorithm to impute missing values and then uses MLR to learn model parameters.
Figure 4 shows the performance of BGCF, DWLS with PD, and FIML for over different percentages of missing values . First, despite a good performance with complete data () DWLS (with PD) deteriorates significantly with an increasing percent of missing values especially for factor loadings, while BGCF and FIML show quite good scalability. Second, our BGCF approach overall outperforms FIML: indistinguishable for interfactor correlations but better for factor loadings.
Two more experiments are provided in Appendix. One concerns incomplete ordinal data with different numbers of categories, showing that BGCF is substantially favorable over DWLS (with PD) and FIML for learning factor loadings, which becomes more prominent with a smaller number of categories. Another one considers incomplete nonparanormal data with different extents of deviation from a Gaussian, which indicates that FIML is rather sensitive to the deviation and only performs well for a slightlynonnormal distribution while the deviation has no influence on BGCF at all. See Appendix for more details.
5 Application to Realworld Data
In this section, we illustrate our approach on the ‘Holzinger & Swineford 1939’ dataset (Holzinger and Swineford, 1939), a classic dataset widely used in the literature and publicly available in the R package lavaan (Rosseel, 2012). The data consists of mental ability test scores of 301 students, in which we focus on 9 out of the original 26 tests as done in Rosseel (2012). A latent variable model that is often proposed to explore these 9 variables is a correlated 3factor model shown in Figure 5, where we rename the observed variables to “Y1, Y2, …, Y9” for simplicity in visualization and to keep it identical to our definition of observed variables (Definition 1). The interpretation of these variables is given in the following list.

Y1: Visual perception;

Y2: Cubes;

Y3: Lozenges;

Y4: Paragraph comprehension;

Y5: Sentence completion;

Y6: Word meaning;

Y7: Speeded addition;

Y8: Speeded counting of dots;

Y9: Speeded discrimination straight and curved capitals.

The summary of the 9 variables in this dataset is provided in Table 1, showing the number of unique values, skewness, and (excess) kurtosis for each variable. From the column of uniques values, we notice that the data are approximately continuous. The average of ‘absolute skewness’ and ‘absolute excess kurtosis’ over the 9 variables are around 0.40 and 0.54 respectively, which is considered to be slightly nonnormal (Li, 2016). Therefore, we choose MLR as the alternative to be compared with our BGCF approach, since these conditions match the assumptions of MLR.
Variables  Unique Values  Skewness  Kurtosis 

Y1  35  0.26  0.33 
Y2  25  0.47  0.35 
Y3  35  0.39  0.89 
Y4  20  0.27  0.10 
Y5  25  0.35  0.54 
Y6  40  0.86  0.84 
Y7  97  0.25  0.29 
Y8  84  0.53  1.20 
Y9  129  0.20  0.31 
We run our Bayesian Gaussian copula factor approach on this dataset. The learned parameter estimates are shown in Figure 5, in which interfactor correlations are on the bidirected edges, factor loadings are in the directed edges, and unique variance for each variable is around the selfreferring arrows. The parameters learned by the MLR approach are not shown here, since we do not know the ground truth so that it is hard to conduct a comparison between the two approaches.
In order to compare the BGCF approach with MLR quantitatively, we consider answering the question: “What is the value of when we observe the values of the other variables, denoted by , given the population model structure in Figure 5?”.
This is a regression problem but with additional constraints to obey the population model structure. The difference from a traditional regression problem is that we should learn the regression coefficients from the modelimplied covariance matrix rather than the sample covariance matrix over observed variables.

For MLR, we first learn the model parameters on the training set, from which we extract the linear regression intercept and coefficients of on . Then we predict the value of based on the values of . See Algorithm 2 for pseudo code of this procedure.

For BGCF, we first estimate the correlation matrix over response variables (the in Definition 1) and the empirical CDF of on the training set. Then we draw latent Gaussian data given and , i.e., . Lastly, we obtain the value of from via , i.e., . See Algorithm 3 for pseudo code of this procedure. Note that we iterate the prediction stage (lines 78) for multiple times in the actual implementation to get multiple solutions to , then the average over these solutions is taken as the final predicted value of . This idea is quite similar to multiple imputation.
The mean squared error (MSE) is used to evaluate the prediction accuracy, where we repeat a 10fold cross validation for 10 times (thus 100 MSE estimates totally). Also, we take as the outcome variable alternately while treating the others as predictors (thus 9 tasks totally). Figure 6 provides the results of BGCF and MLR for all the 9 tasks, showing the mean of MSE with a standard error represented by error bars over the 100 estimates. We see that BGCF outperforms MLR for Tasks 5 and 6 although they perform indistinguishably for the other tasks. The advantage of BGCF over MLR is encouraging, considering that the experimental conditions match the assumptions of MLR. More experiments are done (not shown) after we make the data moderately or substantially nonnormal, suggesting that BGCF is significantly favorable to MLR, as expected.
6 Summary and Discussion
In this paper, we proposed a novel Bayesian Gaussian copula factor (BGCF) approach for learning parameters of CFA models that can handle mixed continuous and ordinal data with missing values. We analyzed the separate identifiability of interfactor correlations , factor loadings , and residual variances , since different researchers may care about different parameters. For instance, it is sufficient to identify for researchers interested in learning causal relations among latent variables (Silva and Scheines, 2006; Silva et al., 2006; Cui et al., 2016), with no need to worry about additional conditions to identify and . Under sufficient identification conditions, we proved that our approach is consistent for MCAR data and empirically showed that it works quite well for MAR data.
In the experiments, our approach outperforms DWLS even under the assumptions of DWLS. Apparently, the approximations inherent in DWLS, such as the use of the polychoric correlation and its asymptotic covariance, incur a small loss in accuracy compared to an integral approach like the BGCF. When the data follow from a more complicated distribution and contain missing values, the advantage of BGCF over its competitors becomes more prominent. Another highlight of our approach is that the Gibbs sampler converges quite fast, where the burnin period is rather short. To further reduce the time complexity, a potential optimization of the sampling process is available (Kalaitzis and Silva, 2013).
There are various generalizations to our inference approach. While our focus in this paper is on the correlated factor models, it is straightforward to extent the current procedure to other class of latent models that are often considered in CFA, such as bifactor models and secondorder models, by simply adjusting the sparsity structure of the prior graph . Also, one may concern models with impure measurement indicators, e.g., a model with an indicator measuring multiple factors or a model with residual covariances (Bollen, 1989), which can be easily solved with BGCF by changing the sparsity pattern of and . Another line of future work is to analyze standard errors and confidence intervals while this paper concentrates on the accuracy of parameter estimates. Our conjecture is that BGCF is still favorable because it naturally transfers the extra variability incurred by missing values to the posterior Gibbs samples: we indeed observed a growing variance of the posterior distribution with the increase of missing values in our simulations. On top of the posterior distribution, one could conduct further studies, e.g., causal discovery over latent factors (Silva et al., 2006; Cui et al., 2018), regression analysis (as we did in Section 5), or other machine learning tasks.
Appendix A: Proof of Theorem 1
Theorem 1 (Consistency of the BGCF Approach).
Let be independent observations drawn from a Gaussian copula factor model. If is complete (no missing data) or contains missing values that are missing completely at random, then
where , , and are parameters learned by BGCF, while , , and are the true ones.
Proof.
If is the response vector’s covariance matrix, then its correlation matrix is , where is a diagonal matrix containing the diagonal entries of . We make use of Theorem 1 from Murray et al. (2013) to show the consistency of . Our factoranalytic prior puts positive probability density almost everywhere on the set of correlation matrices that have a factor decomposition. Then, by applying Theorem 1 in Murray et al. (2013), we obtain the consistency of the posterior distribution on the response vector’s correlation matrix for complete data, i.e.,
(8) 
where is the space restricted by observed data, and is a neighborhood of the true parameter . When the data contain missing values that are completely at random (MCAR), we can also directly obtain the consistency of by again using Theorem 1 in Murray et al. (2013), with an additional observation that the estimation of ordinary and polychoric/polyserial correlations from pairwise complete data is still consistent under MCAR. That is to say, the consistency shown in Equation (8) also holds for data with MCAR missing values.
From this point on, to simplify notation, we will omit adding the tilde to refer to the rescaled matrices , , and . Thus, from now on refers to the correlation matrix of the response vector. and refer to the scaled factor loadings and noise variance respectively.
The Gibbs sampler underlying the BGCF approach has the posterior of (the correlation matrix of the integrated vector ) as its stationary distribution. contains , the correlation matrix of the response random vector, in the upper left block and in the lower right block. Here is the correlation matrix of factors, which implicitly depends on the Gaussian copula factor model from Definition 1 of the main paper via the formula . In order to render this decomposition identifiable, we need to put constraints on , , . Otherwise, we can always replace with and with , where is any invertible matrix, to obtain the equivalent decomposition . However, we have assumed that follows a particular sparsity structure in which there is only a single nonzero entry for each row. This assumption restricts the space of equivalent solutions, since any has to follow the same sparsity structure as . More explicitly, maintains the same sparsity pattern if and only if is a diagonal matrix (Lemma 1).
By decomposing , we get a class of solutions for and , i.e., and , where can be any invertible diagonal matrix. In order to get a unique solution for , we impose two identifying conditions: 1) we restrict to be a correlation matrix; 2) we force the first nonzero entry in each column of to be positive. These conditions are sufficient for identifying uniquely (Lemma 2). We point out that these sufficient conditions are not unique. For example, one could replace the two conditions with restricting the first nonzero entry in each column of to be one. The reason for our choice of conditions is to keep it consistent with our model definition where is a correlation matrix. Under the two conditions for identifying , factor loadings and residual variances are also identified except for the case in which there exists one factor that is independent of all the others and this factor only has two indicators. For such a factor, we have 4 free parameters (2 loadings, 2 residuals) while we only have 3 available equations (2 variances, 1 covariance), which yields an underdetermined system. Therefore, the identifiability of and relies on the observation that a factor has a single or at least three indicators if it is independent of all the others. See Lemmas 3 and 4 for detailed analysis.
Now, given the consistency of and the unique smooth map from to , , and , we obtain the consistency of the posterior mean of the parameter , , and , which concludes our proof.
∎
Lemma 1.
If is a factor loading matrix with only a single nonzero entry for each row, then will have the same sparsity pattern if and only if is diagonal.
Proof.
() We prove the direct statement by contradiction. We assume that has an offdiagonal entry that is not equal to zero. We arbitrarily choose that entry to be . Due to the particular sparsity pattern we have chosen for , there exists such that and , i.e., the unique factor corresponding to the response is . However, we have , which means has a different sparsity pattern from . We have reached a contradiction, therefore is diagonal.
() If is diagonal, i.e., , then . This means that , so the sparsity pattern is preserved. ∎
Lemma 2 (Identifiability of ).
Given the factor structure defined in Section 3 of the main paper, we can uniquely recover from if 1) we constrain to be a correlation matrix; 2) we force the first element in each column of to be positive.
Proof.
Here we assume that the model has the stated factor structure, i.e., that there is some , , and such that . We then show that our chosen restrictions are sufficient for identification using an argument similar to that in Anderson and Rubin (1956).
The decomposition constitutes a system of equations:
(9)  
where , and is the map from a response variable to its corresponding factor. Looking at the equation system in (9), we notice that each factor correlation term , appears only in the equations corresponding to response variables indexed by and such that and or vice versa. This suggests that we can restrict our analysis to submodels that include only two factors by considering the submatrices of that only involve those two factors. To be more precise, the idea is to look only at the equations corresponding to the submatrix , where is the preimage of under . Indeed, we will show that we can identify each individual correlation term corresponding to pairs of factors only by looking at these submatrices. Any information concerning the correlation term provided by the other equations is then redundant.
Let us then consider an arbitrary pair of factors in our model and the corresponding submatrices of , , , and . (The case of a single factor is trivial.) In order to simplify notation, we will also use , , , and to refer to these submatrices. We also reindex the two factors involved to and for simplicity. In order to recover the correlation between a pair of factors from , we have to analyze three separate cases to cover all the bases (see Figure 7 for examples concerning each case):

The two factors are not correlated, i.e., . (There are no restrictions on the number of response variables that the factors can have.)

The two factors are correlated, i.e., , and each has a single response, which implies that and .

The two factors are correlated, i.e., , but at least one of them has at least two responses.

Case 1: If the two factors are not correlated (see the example in the left panel of Figure 7), this fact will be reflected in the matrix . More specifically, the offdiagonal blocks in , which correspond to the covariance between the responses of one factor and the responses of the other factor, will be set to zero. If we notice this zero pattern in , we can immediately determine that .
Case 2: If the two factors are correlated and each factor has a single associated response (see the middle panel of Figure 7), the model reduces to a Gaussian Copula model. Then, we directly get since we have put the constraints if has a single indicator .
Case 3: If at least one of the factors (w.l.o.g., ) is allowed to have more than one response (see the example in the right panel of Figure 7), we arbitrarily choose two of these responses. We also require one response variable corresponding to the other factor (). We use , and to denote the loadings of these response variables, where . From Equation (9) we have:
Since we are in the case in which , which automatically implies that , we can divide the last two equations to obtain . We then multiply the result with the first equation to get . Without loss of generality, we can say that is the first entry in the first column of , which means that . This means that we have uniquely recovered and .
We can also assume without loss of generality that is the first entry in the second column of , so . If has at least two responses, we use a similar argument to the one before to uniquely recover . We can then use the above equations to get . If has only one response, then , which means that , so again is uniquely recoverable and we can obtain from the equations above.
Thus, we have shown that we can correctly determine only from in all three cases. By applying this approach to all pairs of factors, we can uniquely recover all pairwise correlations. This means that, given our constraints, we can uniquely identify from the decomposition of .
∎
Lemma 3 (Identifiability of ).
Given the factor structure defined in Section 3 of the main paper, we can uniquely recover from if 1) we constrain to be a correlation matrix; 2) we force the first element in each column of to be positive; 3) when a factor is independent of all the others, it has either a single or at least three indicators.
Proof.
Compared to identifying , we need to consider another case in which there is only one factor or there exists one factor that is independent of all the others (the former can be treated as a special case of the latter). When such a factor only has a single indicator, e.g., in the left panel of Figure 7, we directly identify because of the constraint . When the factor has two indicators, e.g., in the left panel of Figure 7, we have four free parameters (, , , and ) while we can only construct three equations from (, , and ), which cannot give us a unique solution. Now we turn to the threeindicator case, as shown in Figure 8. From Equation (9) we have:
We then have , which has a unique solution for together with the second constraint , after which we can naturally get the solutions to and . For the other cases, the proof follows the same line of reasoning as Lemma 2.

∎
Lemma 4 (Identifiability of ).
Given the factor structure defined in Section 3 of the main paper, we can uniquely recover from if 1) we constrain to be a correlation matrix; 2) when a factor is independent of all the others, it has either a single or at least three indicators.
Proof.
We conduct our analysis case by case. For the case where a factor has a single indicator, we trivially set . For the case in Figure 8, it is straightforward to get from (the same for and ). Another case we need to consider is Case 3 in Figure 7, where we have (see analysis in Lemma 2), based on which we obtain . By applying this approach to all single factors or pairs of factors, we can uniquely recover all elements of .∎
Appendix B: Extended Simulations
This section continues the experiments in Section 4 of the main paper, in order to check the influence of the number of categories for ordinal data and the extent of nonnormality for nonparanormal data.
B1: Ordinal Data with Different Numbers of Categories
In this subsection, we consider ordinal data with various numbers of categories , in which the sample size and missing values percentage are set to and respectively. Figure 9 shows the results obtained by BGCF (Bayesian Gaussian copula factor), DWLS (diagonally weighted least squares) with PD (pairwise deletion), and FIML (full information maximum likelihood), providing the mean of ARB (average relative bias) and the mean of RMSE (root mean squared error) with 95% confidence interval over 100 experiments for (a) interfactor correlations and (b) factor loadings. In the case of two categories, FIML underestimates factor loadings dramatically, DWLS obtains a moderate bias, while BGCF just gives trivial bias. With an increasing number of categories, FIML gets closer and closer to BGCF, but still BGCF is favorable.
B2: Nonparanormal Data with Different Extents of Nonnormality
In this subsection, we consider nonparanormal data, in which we use the degrees of freedom of a distribution to control the extent of nonnormality (see Section 5.1 of the main paper for details). The sample size and missing values percentage are set to and respectively, while the degrees of freedom varies .
Figure 10 shows the results obtained by BGCF, DWLS with PD, and FIML, providing the mean of ARB (left panel) and the mean of RMSE with 95% confidence interval (right panel) over 100 experiments for (a) interfactor correlations and (b) factor loadings. The major conclusion drawn here is that, while a nonparanormal transformation has no effect on our BGCF approach, FIML is quite sensitive to the extent of nonnormality, especially for factor loadings.
Footnotes
 The code including those used in simulations and realworld applications is provided in https://github.com/cuiruifei/CopulaFactorModel.
References
 Anderson, T.W., Rubin, H.: Statistical inference in factor analysis. In: Proceedings of the Third Berkeley Symposium on Mathematical Statistics and Probability, Volume 5: Contributions to Econometrics, Industrial Research, and Psychometry, pp. 111–150. University of California Press, Berkeley, Calif. (1956)
 Arbuckle, J.L.: Full information estimation in the presence of incomplete data. Advanced structural equation modeling: Issues and techniques 243, 277 (1996)
 Barendse, M., Oort, F., Timmerman, M.: Using exploratory factor analysis to determine the dimensionality of discrete responses. Struct. Equ. Modeling. 22(1), 87–101 (2015)
 Barnard, J., McCulloch, R., Meng, X.L.: Modeling covariance matrices in terms of standard deviations and correlations, with application to shrinkage. Stat. Sinica. pp. 1281–1311 (2000)
 Bollen, K.: Structural equations with latent variables. NY Wiley (1989)
 Browne, M.W.: Asymptotically distributionfree methods for the analysis of covariance structures. Brit. J. Math. Stat. Psy. 37(1), 62–83 (1984)
 Byrne, B.M.: Structural equation modeling with EQS: Basic concepts, applications, and programming. Routledge (2013)
 Castro, L.M., Costa, D.R., Prates, M.O., Lachos, V.H.: Likelihoodbased inference for Tobit confirmatory factor analysis using the multivariate Studentt distribution. Stat. Comput. 25(6), 1163–1183 (2015)
 Cui, R., Groot, P., Heskes, T.: Copula PC algorithm for causal discovery from mixed data. In: Joint European Conference on Machine Learning and Knowledge Discovery in Databases, pp. 377–392. Springer (2016)
 Cui, R., Groot, P., Heskes, T.: Learning causal structure from mixed data with missing values using Gaussian copula models. Stat. Comput. (2018)
 Curran, P.J., West, S.G., Finch, J.F.: The robustness of test statistics to nonnormality and specification error in confirmatory factor analysis. Psychol. Methods. 1(1), 16 (1996)
 DiStefano, C.: The impact of categorization with confirmatory factor analysis. Struct. Equ. Modeling. 9(3), 327–346 (2002)
 Hoff, P.D.: Extending the rank likelihood for semiparametric copula estimation. Ann. Stat. pp. 265–283 (2007)
 Holzinger, K.J., Swineford, F.: A study in factor analysis: the stability of a bifactor solution. Suppl. Educ. Monogr. 48 (1939)
 Jones, B., Carvalho, C., Dobra, A., Hans, C., Carter, C., West, M.: Experiments in stochastic computation for highdimensional graphical models. Stat. Sci. pp. 388–400 (2005)
 Jöreskog, K.G.: A general approach to confirmatory maximum likelihood factor analysis. Psychometrika 34(2), 183–202 (1969)
 Jöreskog, K.G.: Structural equation modeling with ordinal variables using LISREL. Tech. rep., Technical report, Scientific Software International, Inc., Lincolnwood, IL (2005)
 Kalaitzis, A., Silva, R.: Flexible sampling of discrete data correlations without the marginal distributions. In: Advances in Neural Information Processing Systems, pp. 2517–2525 (2013)
 Kaplan, D.: Structural equation modeling: Foundations and extensions, vol. 10. Sage Publications (2008)
 Kolar, M., Xing, E.P.: Estimating sparse precision matrices from data with missing values. In: International Conference on Machine Learning (2012)
 Lancaster, G., Green, M.: Latent variable techniques for categorical data. Stat. Comput. 12(2), 153–161 (2002)
 Li, C.H.: Confirmatory factor analysis with ordinal data: Comparing robust maximum likelihood and diagonally weighted least squares. Behav. Res. Methods. 48(3), 936–949 (2016)
 Little, R.J., Rubin, D.B.: Statistical analysis with missing data (1987)
 Lubke, G.H., Muthén, B.O.: Applying multigroup confirmatory factor models for continuous outcomes to likert scale data complicates meaningful group comparisons. Struct. Equ. Modeling. 11(4), 514–534 (2004)
 Marsh, H.W., Hau, K.T., Balla, J.R., Grayson, D.: Is more ever too much? the number of indicators per factor in confirmatory factor analysis. Multivar. Behav. Res. 33(2), 181–220 (1998)
 MartínezTorres, M.R.: A procedure to design a structural and measurement model of intellectual capital: an exploratory study. Informa. Manage. 43(5), 617–626 (2006)
 Murphy, K.P.: Conjugate Bayesian analysis of the Gaussian distribution. def 1(22), 16 (2007)
 Murray, J.S., Dunson, D.B., Carin, L., Lucas, J.E.: Bayesian Gaussian copula factor models for mixed data. J. Am. Stat. Assoc. 108(502), 656–665 (2013)
 Muthén, B.: A general structural equation model with dichotomous, ordered categorical, and continuous latent variable indicators. Psychometrika 49(1), 115–132 (1984)
 Muthén, B., du Toit, S., Spisic, D.: Robust inference using weighted least squares and quadratic estimating equations in latent variable modeling with categorical and continuous outcomes. Psychometrika (1997)
 Muthén, L.: Mplus user’s guide, (Muthén & Muthén, los angeles). Mplus User’s Guide,(Muthén & Muthén, Los Angeles) (2010)
 Olsson, U.: Maximum likelihood estimation of the polychoric correlation coefficient. Psychometrika 44(4), 443–460 (1979)
 Poon, W.Y., Wang, H.B.: Latent variable models with ordinal categorical covariates. Stat. Comput. 22(5), 1135–1154 (2012)
 Rhemtulla, M., BrosseauLiard, P.É., Savalei, V.: When can categorical variables be treated as continuous? A comparison of robust continuous and categorical SEM estimation methods under suboptimal conditions. Psychol. Methods. 17(3), 354 (2012)
 Rosseel, Y.: lavaan: An R package for structural equation modeling. J. Stat. Softw. 48(2), 1–36 (2012)
 Roverato, A.: Hyper inverse Wishart distribution for nondecomposable graphs and its application to Bayesian inference for Gaussian graphical models. Scan. J. Stat 29(3), 391–411 (2002)
 Rubin, D.B.: Inference and missing data. Biometrika pp. 581–592 (1976)
 Schafer, J.L.: Analysis of incomplete multivariate data. CRC press (1997)
 Schafer, J.L., Graham, J.W.: Missing data: our view of the state of the art. Psychol. Methods. 7(2), 147 (2002)
 Silva, R., Scheines, R.: Bayesian learning of measurement and structural models. In: International Conference on Machine Learning, pp. 825–832 (2006)
 Silva, R., Scheines, R., Glymour, C., Spirtes, P.: Learning the structure of linear latent variable models. J. Mach. Learn. Res. 7(Feb), 191–246 (2006)
 YangWallentin, F., Jöreskog, K.G., Luo, H.: Confirmatory factor analysis of ordinal variables with misspecified models. Struct. Equ. Modeling. 17(3), 392–423 (2010)