Opaque Prior Distributions in Bayesian Latent Variable Models

Illustration

An example comes from the popular “Political Democracy” model originally described by Bollen (1989), shown below in lavaan syntax (Rosseel, 2012). Intended to describe countries’ relationships between their levels of industrialization and democracy, the model contains four observed variables that are measured once during 1960 ( $y1,...,y4$) and again during 1965 ( $y5,...,y8$). The residuals of the 1960 variables are allowed to correlate with the residuals of the corresponding 1965 variables. Additionally, there exists a pair of similar variables collected during 1960 and again during 1965, leading to two more residual correlations. The full structure of the residual covariance matrix is shown in Figure 1.

As described earlier, we could elect to put a univariate prior distribution on each variance (or standard deviation or precision) in this matrix, and also on the six residual correlations that are not fixed to zero. But this is problematic because the univariate priors on correlations can yield non-positive definite correlation matrices. The specific problem that occurs depends on the software package. For example, JAGS will stop as soon as it encounters a correlation matrix that is not positive definite. This means that univariate priors on correlations cannot often be used. On the other hand, Stan will report that a non-positive definite matrix was encountered, reject it, and continue sampling. We are left with univariate priors that are collectively constrained to be positive definite. If we consider only the space of positive definite correlation matrices under these priors, then the priors are usually more informative than we originally specified. That is, the prior distributions are opaque: the analyst specifies a set of priors that are different from the implied prior distributions, which must obey the constraint of positive definiteness.

Click to enlarge

Figure 1

Political Democracy Model, Structure of Residual Covariance Matrix With Free Parameters Marked by $\times$.

model <- '
 # measurement model
  ind60 =˜ x1 + x2 + x3
  dem60 =˜ y1 + a*y2 + b*y3 + c*y4
  dem65 =˜ y5 + a*y6 + b*y7 + c*y8
 # regressions
  dem60 ˜ ind60
  dem65 ˜ ind60 + dem60
 # residual correlations
  y1 ˜˜ y5
  y2 ˜˜ y4 + y6
  y3 ˜˜ y7
  y4 ˜˜ y8
  y6 ˜˜ y8
'

How can we characterize the implied prior distributions? We could simply generate thousands of correlation matrices from the prior, then discard matrices that are not positive definite, then visualize what is left. As an example of this, we specified a Uniform( $-1,1$) prior for each free correlation in Figure 1. We then generated 100,000 correlation matrices with the desired structure, discarding the 57,818 matrices that were not positive definite. Finally, we examined the distributions of the remaining matrices, with a histogram for a single correlation parameter appearing in Figure 2. We see that the resulting distribution is no longer uniform; it is approximately symmetric around 0, with more density near 0 than near $-1$ or $1$. This is because our uniform priors did not account for the fact that correlation matrices must be positive definite.

Click to enlarge

Figure 2

Bollen Model, Implied Prior Distribution of a Single Correlation Parameter After Accounting for Positive Definite Constraints

Positive Definiteness

To technically describe the positive definite constraint here, we can simultaneously permute the rows and columns of the correlation matrix to obtain a block diagonal matrix. A block diagonal matrix is useful because the determinant of the full matrix is the product of determinants of each individual block within the matrix. This allows us to examine whether or not the full matrix is positive definite, by working with submatrices of smaller dimension.

The Cuthill-McKee algorithm (Cuthill & McKee, 1969) allows us to automatically find an appropriate, block-diagonal permutation. After carrying out this algorithm (via the netprioR package; Schmich, 2022) and permuting the rows and columns, we arrive at the matrix in Figure 3. This shows that, to keep the full matrix positive definite, we only need to worry about the $4\times 4$ block in the lower right that involves $y2$, $y4$, $y6$, and $y8$. The priors on the $(y3,y7)$ and $(y1,y5)$ correlations have no influence on the positive definiteness of the full matrix, because a $2\times 2$ matrix is positive definite for any correlation in $(-1,1)$, so that we can safely put a uniform prior on each of those correlations.

Click to enlarge

Figure 3

Political Democracy Model, Structure of Permuted Correlation Matrix With Free Parameters Marked by $\times$.

After applying traditional rules for computing matrix determinants, we can express the determinant of the $4\times 4$ block of correlations as

This shows analytically how the positive definite constraint of our correlation matrix influences our univariate priors. The univariate priors are collectively constrained by the requirement that Equation (1) be greater than 0. For MCMC algorithms that sample the correlations individually, each correlation is impacted by the current values of the other correlations. For example, if the current values of $r_1$, $r_2$, and $r_3$ are .9. .1, and .1, respectively, then $r_4$ must be less than .9 in order to maintain a positive definite matrix. But if the first three correlations assume different values, then it would be possible to observe values of $r_4$ above.9.

Implications for Bayes Factors

To show how this issue becomes problematic for applied work, we consider the calculation of Bayes factors using the Bollen model. Imagine that we wish to know whether the two residual correlations involving $y2$ are necessary. To address this question, we could compute a Bayes factor comparing a model that includes those two residual correlations, to a model without those residual correlations. A computationally-cheap way to do this is the Savage-Dickey method (Dickey & Lientz, 1970; Wagenmakers et al., 2010). This involves comparing the prior distributions of the two residual covariances to their posterior distributions, focusing on the point at which those covariances equal 0.

The Savage-Dickey method is slightly more complicated than usual because the two residual correlations in question influence the possible values that other residual correlations can take (due to positive definiteness). This changes the prior distributions on other residual correlations, as we move from a model with the two focal correlations freed to a model with the two focal correlations fixed. Heck (2019) discusses that, in such a situation, we need to supply a correction term to the usual Savage-Dickey calculation that accounts for the change in prior distributions (also see Verdinelli & Wasserman, 1995).

If we ignore (or do not realize) all of the above and set Uniform priors on each individual correlation, then the prior density of the two focal correlations at 0 equals $0.25$. This corresponds to a log-density of $-1.39$. In contrast, the true joint prior density on the two correlations of the full model (with the two focal correlations free, respecting the positive definite constraint) can be approximated using the positive definite correlation matrices that we randomly generated earlier. This approximations involves a density estimation method from the ks package in R (Duong, 2022). Our approximate joint density at 0 is now 0.46, corresponding to a log-density of $-0.78$. And we have an additional correction term to account for the fact that the priors for the $y4$– $y8$ and for the $y6$– $y8$ correlations are impacted by whether or not the $y2$ residual correlations are fixed to 0.

These evaluations lead to two separate Bayes factors for the model with correlations, relative to the model without correlations. Using our incorrect priors that do not account for positive definite constraints, we obtain a log-Bayes factor of 4.66 in favor of the model with correlations. Using our priors that do account for positive definite constraints, along with the correction from Heck (2019), we obtain a log-Bayes factor of 5.54 in favor of the model with correlations. Both of these Bayes factors provide support for the model with correlations, but possibly at different levels of evidence. For example, if we subscribe to the rules of thumb offered by Kass and Raftery (1995), then these two Bayes factors lead us to conclude “strong” evidence using the incorrect prior calculation, and “very strong” evidence using the correct calculation.

We agree with you, the reader, that these cutoffs are arbitrary and that the Bayes factors do not differ by very much. But the point is that the Bayes factor systematically differ depending on how we compute prior densities. These differences will sometimes lead to different substantive conclusions in practice, with the easier computation (i.e., ignoring the positive definite constraint) being incorrect.

Solutions

While it was straightforward to visualize implied prior distributions in the Bollen model, the process becomes inefficient for correlation matrices whose dimension is larger than 3 or 4. For such correlation matrices, few of the randomly-generated matrices will be positive definite, and it will take a long time to obtain a sufficient number of positive definite matrices to describe the implied prior. Additionally, every unique structure of correlation matrix will have a unique positive definite constraint, similar to the one from Equation (1). So we desire more general solutions that do not involve random generation of correlation matrices. We describe three solutions below that differ in complexity and in the extent to which they fully solve the problem.

Illustration

To illustrate the interaction between sign indeterminacy and prior distributions, it is sufficient to consider the usual confirmatory factor model that is fit to the Holzinger and Swineford (1939) data. We suspect that most people reading this far know the dataset, which contains scores on various tests of mental ability. We are focusing on the 3-factor model that is traditionally fit to the version of the data from lavaan (Rosseel, 2012), where each factor is associated with 3 observed variables.

The Holzinger-Swineford factor model has five types of model parameters: intercepts, loadings, factor standard deviations, factor correlations, and residual standard deviations. We assign true (“population”) values to all these parameters. Intercepts receive true values of 0, factor standard deviations receive true values of 1, factor correlations receive true values of 0, and residual standard deviations receive true values of 1. Finally and importantly, loadings receive true values of $-1$.

Using these true values, we generated a dataset of 1,000 observations and re-fit the 3-factor model back to the data (where true values were treated as unknown). We used common, non-informative priors for the parameters of the estimated model, which are currently the defaults in blavaan (Merkle et al., 2021; Merkle & Rosseel, 2018).

The resulting posterior distributions of the loadings appear in Figure 4. These distributions are centered near $+1$, with the distributions being fully on the positive side of the space. So we have a situation where the true loading values were $-1$, we used noninformative priors that exhibit little influence on the posterior, and the posterior distributions of loading values are nowhere near the true values.

Click to enlarge

Figure 4

Estimated Posterior Distributions of Loadings, Holzinger-Swineford 3-Factor Model

At this point, readers might object that this just illustrates sign indeterminacy. And posterior inferences about factor loadings are not generally impacted here. But the example highlights that, for loadings, a noninformative prior centered at 0 usually ignores the identification constraint that was chosen. That is, when employing noninformative priors, we are usually attempting to say that we have no idea about the loadings’ values, or to avoid influencing the results of the model estimation. But the prior ignores the fact that (i) we typically fix a loading to be positive for identification, and (ii) observed variables are usually positively correlated, so that we can expect other loadings to be positive. So our original priors, which were intended to be noninformative, may actually state that an implausible part of the parameter space is plausible (also see Gelman et al., 2017; Seaman et al., 2012). Further, as mentioned earlier, the prior distribution of relabelled factor loadings generally differs from the prior distribution of un-relabelled factor loadings, which can impact post-estimation computations that rely on evaluation of the prior.

Implications for MCMC Validation

Sign indeterminacy also complicates MCMC algorithm validation, which is used to ensure that MCMC samplers are working correctly. The MCMC validation process is difficult even without sign indeterminacy, because randomness is inherent in MCMC sampling. This means that we cannot simply examine whether the posterior means and standard deviations match other samplers to many decimal places. While it is possible to obtain analytic posterior distributions for certain models, analytic results are the exception instead of the rule for models estimated via MCMC.

Simulation-Based Calibration (SBC; Modrák et al., 2022; Talts et al., 2018) is an MCMC validation method that has received recent attention and implementation. Given a model of interest (including likelihood and priors), simulation-based calibration can be described in four steps:

If the MCMC sampler is working correctly, then the posteriors from Step 3 should look like the priors from which we started. We can graphically examine this idea by comparing the posterior means from Step 3 to the parameter values from Step 1; we should see an identity line when plotting the parameter values against the posterior means.

Using the same model from the previous section, we used the SBC package (Kim et al., 2022) to conduct simulation-based calibration under two sets of priors, with 1000 simulated data sets each. Set 1 was exactly the same as the noninformative priors from the previous section. Set 2 was also similar to the previous section, differing only in the priors for the loadings: instead of Normal(0, 100), the priors for loadings were Normal(1, 1/16) (where the second number is a variance). This is an informative prior reflecting the belief that all loadings should be similar to one another (recall that a single loading is being fixed to 1 for identification).

Results for the Set 1 and Set 2 priors are shown in Figure 5 and Figure 6, respectively. In both figures, the parameter values simulated from the prior are on the x-axis, and the posterior means estimated from the artificial data are on the y-axis. Each panel represents a factor loading (there are 6 free loadings in the model). Each point represents a replication, and “correct” MCMC algorithms should lead to points that are crowded around the blue diagonal line.

Click to enlarge

Figure 5

Prior Parameter Values Versus Posterior Means for N(0,100) Priors on Loadings

Click to enlarge

Figure 6

Prior Parameter Values Versus Posterior Means for N(1, $\frac{1}{16}$) Priors on Loadings

Figure 5 shows that, for noninformative priors on factor loadings, the points follow a V or an X pattern. The V pattern occurs for loadings that are constrained to be positive during estimation, while the X pattern occurs for the remaining loadings. This means that the estimation algorithm often recovers the true parameter value, but it also often flips the sign of the parameter value. While this is not a problem if we realize the sign indeterminacy issue, it is easy to overlook in the context of simulation-based calibration. For example, typical SBC summaries involve cumulative distributions, which would lead us to conclude that there are problems with the MCMC algorithm, and which would not provide clear clues about sign indeterminacy. See Merkle et al. (2021) for some discussion of similar analyses.

Figure 6, on the other hand, is closer to what we would hope to see from a correct MCMC algorithm. In this figure, the points generally form a cloud around the blue diagonal line, implying that the posterior means match the parameter values that generated the data (note that the axis limits have changed, as compared to the previous figure). The informative priors on factor loadings work similarly to use of truncated normal priors, where the signs of the loadings were generally restricted to be positive. But informative priors differ from the truncated normals in that they allow for the possibility of negative loadings.

Solutions

If researchers are aware of how their software handles sign indeterminacy, then they can potentially avoid the issues described here and successfully employ noninformative prior distributions. Barring that, we advise researchers to explicitly consider the loadings’ expected signs when setting prior distributions. In many models, we expect that all the observed variables corresponding to a factor will have the same direction of relationship with that factor. Additionally, a single loading is often set to 1 for identification. In a situation like this, it is often reasonable to place priors on the free loadings that have a mean of 1 and a standard deviation of, say, .5. These priors look very informative at first glance, as compared to, say, a Normal prior with a mean of 0 and a variance of 10,000. But the suggested priors better represent what the researcher knows about the signs of the loadings, combined with the fact that some loadings are being fixed to 1.

Instead of fixing a single loading to 1 for identification, researchers may fix the latent variance to 1. In this case, as described by Peeters (2012), one loading per latent variable must be constrained to be positive (or negative) in order to achieve identification. If that loading has a Normal prior whose support includes the negative (positive) reals, then we again face a conflict where our stated prior is not the implied prior. The implied prior for the sign-constrained loading may become truncated normal or truncated $t$, depending on other details of the MCMC implementation (see Ghosh & Dunson, 2009), while the priors on remaining loadings are as stated. While a truncated distribution arises from the identification constraint here, we think that truncated priors on all loadings should be avoided, because they prevent the loadings’ posterior distributions from overlapping with zero. The relabeling strategies described by Erosheva and Curtis (2017) should be preferred.

There are alternative prior distributions that avoid the issue and/or that make it easier to specify informative prior distributions, though they are not readily available in popular software. Graves and Merkle (2022) studied prior distributions on ratios of factor loadings, as well as prior specification under effect coding (Little et al., 2006). Ratios of loadings avoid the sign-switching issue (e.g., the ratio of two negative loadings remains positive), and effect coding can make it easier to specify informative prior distributions. Additionally, in the case of exploratory factor analysis, Heaps (2022) recently proposed a prior distribution for the common factor covariance matrix. That is, factor analysis models typically imply a covariance matrix of the form $\mathbf{\Lambda }{\mathbf{\Lambda }}^{\mathrm{\prime }}+\mathbf{\Psi }$, where $\mathbf{\Lambda }$ is the factor loading matrix with many fewer columns than rows. Heaps (2022) proposes to place a matrix normal prior on $\mathbf{\Lambda }{\mathbf{\Lambda }}^{\mathrm{\prime }}$, which encodes researcher knowledge about shared variation in the observed variables. The term $\mathbf{\Lambda }{\mathbf{\Lambda }}^{\mathrm{\prime }}$ remains invariant across factor loadings’ signs and rotations, so it avoids the issues described here. The proposed priors can also shrink sets of loadings towards zero, which is similar to confirmatory factor analysis. It remains to be seen how this prior could be translated to more general SEMs, where the model-implied covariance matrix becomes more complex.

Example

Say that a researcher fits a factor analysis model to a set of 4-category ordinal variables, and that she specifies a Normal(0,5) prior on all threshold parameters in the model. Because there are four categories per variable, we require three order-constrained thresholds per variable. We wish to know what these priors look like, after accounting for the order constraints.

We consider two ways that we could translate Normal(0,5) priors to three ordered parameters (also see Padgett et al., 2023). First, we could imagine drawing three separate variates from the normal distribution, then ordering the variates to obtain ordered thresholds. Second, we could imagine drawing the first (lowest) threshold from a Normal(0,5), then adding a Lognormal(0,5) variate to that threshold in order to obtain the second threshold. Once we have obtained the second threshold, we could add another Lognormal variate to obtain the third threshold. Lognormal variables can only take positive values, so we are guaranteed to have an ordered set of thresholds under this approach.

For both of these translations, the threshold parameters’ prior distributions differ from the Normal(0,5) distribution that the researcher originally declared. We expand on this point below, separately for the two methods.

Reordering

When we draw three values from the Normal distribution and then order them, the act of ordering influences the resulting prior distributions. The specific distributions can be described via statistical theory on order statistics. For our example, the Normal(0,5) priors translate into the following probability density functions (pdfs) for individual thresholds:

$$\begin{array}{rl}p(g_1)& =3\times \varphi (g_1/5)\times [1-\mathrm{\Phi }(g_1/5){]}^2\\ p(g_2)& =6\times \varphi (g_2/5)\times \mathrm{\Phi }(g_2/5)\times [1-\mathrm{\Phi }(g_1/5)]\\ p(g_3)& =3\times \varphi (g_2/5)\times \mathrm{\Phi }(g_2/5{)}^2,\end{array}$$

where $\varphi ()$ is the standard normal pdf and $\mathrm{\Phi }()$ is the standard normal cdf. These distributions are visualized in Figure 7, with the Normal(0,5) distribution overlayed for comparison. The figure shows that the prior for the lower threshold ( $g_1$) is centered below 0, while the prior for the upper threshold ( $g_3$) is centered above 0. None of the three distributions matches the Normal(0,5), despite the fact that the researcher declared a Normal(0,5) prior for all three parameters.

Click to enlarge

Figure 7

Priors of Threshold Parameters Whose Stated Priors are Normal(0,5), as Implied by Order Constraints

Lognormals

The blavaan package uses the lognormal approach mentioned previously. The Normal(0,5) prior goes on the lowest threshold parameter per item. Then, Normal(0,5) priors are placed on log-differences between subsequent parameters. For example, considering the three thresholds from the example, we would have

$$\begin{array}{rl}{g}_1& \sim \text{Normal}(0,5)\\ \log (g_2-g_1)& \sim \text{Normal}(0,5)\\ \log (g_3-g_2)& \sim \text{Normal}(0,5),\end{array}$$

where we could alternatively say that the differences between thresholds follow Lognormal priors.

Just like the previous section, the above priors can be translated to priors on individual thresholds. It is obvious that the prior for $g_1$ continues to be the stated prior, which is Normal(0,5). But the priors for $g_2$ and $g_3$ involve the sum of a normal distribution and Lognormal distribution(s). The resulting distributions can be written as

$$\begin{array}{rl}p(g_2)& ={\displaystyle {\int }_{-\mathrm{\infty }}^{g_2}\varphi ((g_1-\mu )/\sigma )\times \text{logN}(g_2-g_1,\mu ,\sigma )\mathrm{\partial }{g}_1}\\ p(g_3)& ={\displaystyle {\int }_{-\mathrm{\infty }}^{g_3}{\int }_{-\mathrm{\infty }}^{g_2}\varphi ((g_1-\mu )/\sigma )\times \text{logN}(g_2-g_1,\mu ,\sigma )\times \text{logN}(g_3-g_2,\mu ,\sigma )\mathrm{\partial }{g}_1\mathrm{\partial }{g}_2,}\end{array}$$

where $\text{logN}(x,\mu ,\sigma )$ is the density function of the Lognormal distribution with mean $\mu$ and standard deviation $\sigma$ (both on the log scale), evaluated at $x$.

While the priors for $g_2$ and $g_3$ do not have nice forms, we can numerically approximate the integrals to visualize them. These distributions are shown in Figure 8, which is arranged similarly to Figure 7. As stated before, the lowest threshold, $g_1$, has the stated Normal(0,5) distribution. The peaks of the distributions of the remaining thresholds are closer to that of $g_1$, as compared to the distributions in Figure 7. The distributions of $g_2$ and $g_3$ are also skewed to the right, reflecting the fact that we are adding a positive variate to the distribution of $g_1$.

Click to enlarge

Figure 8

Priors of Threshold Parameters Whose Stated Priors are Normal(0,5), as Implied by Placing Priors on Log-Differences

Solution

Unlike the previous issues with positive definite constraints and sign constraints, researchers do not have to consider changing their priors in order to address order constraints. The main solution is to be aware of the fact that, if one places univariate prior distributions on a set of order-constrained parameters, then some translation will take place to ensure that the parameters are ordered correctly. And this translation will influence the implied prior distribution of each parameter. It is worthwhile to understand how this is handled by one’s software, especially for prior predictive assessments, Bayes factor calculation, and simulation-based calibration methods.