This post is from Bob.
Mitzi and I were swotting up on structural equation models (SEM) for our class this past Monday at the Modern Modeling and Methods (M3) conference at Fordham University. It was a lot of fun and now I think I understand SEM notation. I really like these applied conferences and this was a group of psychometrician, econometricians, and sociometricians. Many if not most of them thought about models in terms of SEM, so we thought we should figure it out. But I was left with a concern you may be able to help me sort out.
The example
The first worked example in Ken Bollen’s seminal 1979 textbook on SEM is a study of how industrialization relates to democracy. It comes from his paper,
- Bollen, Kenneth A. (1979). “Political Democracy and the Timing of Development.” American Sociological Review, 44(4).
and was reprised in his book
- Bollen, Kenneth A. (1989). Structural Equations with Latent Variables. Wiley.
I had the pleasure of sitting across from Ken at the invited speakers dinner at the conference, so I’m glad I looked into SEM before that. Good news for the SEM devotees—he released a completely revised guide to SEM a few months ago.
- Bollen, Kenneth A. 2026. Elements of Structural Equation Models. Cambridge University Press.
The data and parameters
The data consists of eleven covariates (called “indicators” in SEM) for each of 75 countries. Four of the covariates are related to democracy in 1960 (y1, y2, y3, y4), the same four measurements were taken again again in 1965 (y5, y6, y7, y8) , and there were three measurements of industrialization in 1960 (x1, x2, x3).
The SEM model the original researcher came up with here assumes three latent scalars per country, industrialization in 1960 (IND60), level of democracy in 1960 (DEM60), and level of democracy in 1965 (DEM65). These latent parameters are related in the following way: democracy in 1960 is a regression on industrialization in 1960, and democracy in 1965 is a regression on both democracy in 1960 and industrialization in 1960.
The covariates are then modeled like a seemingly unrelated regression in econometrics. The four democracy 1965 parameters are treated as regressions on the latent level of democracy in 1965, and similarly for the democracy in 1960, and industrialization in 1960.
Rather than independent errors, a SEM model explicitly indicates with arrows which pairs of observations are allowed to have non-zero correlation in the covariance matrix for the observations. The three industrialization observations are assumed to have zero correlation—there are no arrows between any of the three measurements in the SEM diagram. Each of the four measurements in 1960 is assumed to covary with the same measurement taken in 1965. In addition, the second and fourth measurement in each year are assumed to be correlated with each other, which leads to a box-like structure.
The SEM diagram
Here are the arrows in the diagram, where I’m not using their standard LISREL notation, but writing them in R expression syntax to indicate what is regressed on what. In their graphical notation, just replace ~ with <-. All three latent variables and all eleven measurements are indexed by country.
IND60 DEM60 ~ IND60 DEM65 ~ DEM60, IND60 x1, x2, x3 ~ IND60 y1, y2, y3, y4 ~ DEM60 y5, y6, y7, y8 ~ DEM65
The covariance structure is indicated by stating which pairs of measurements are modeled with non-zero correlation. The first four just pair the measurements of the same thing across 1960 and 1965.
y1 <-> y5 y2 <-> y6 y3 <-> y7 y4 <-> y8
The last pair of correlations are within 1960 and within 1965.
y2 <-> y4 y6 <-> y8
Together, these induce an odd box structure, where y2 is correlated with y6 and y4, both of which are correlated with y8, but y2 and y8 are assumed to have zero correlation.
y2 <-> y6 ^ ^ | | v v y4 <-> y8
Stan implementation
We didn’t get this far in my half of the class, so I will share here the Stan Playground example where I fit Bollen’s example (you can get the data and the Stan model through the Playground link:
It gets the right answer compared to lavaan/blavaan, which is nice. In the Stan code, xi is IND60 and eta1, eta2 are DEM60, DEM65. The relation among the latent parameters are modeled directly as regressions. The correlations among the observations are modeled using soft zeroing, where I just put a tight prior around zero on the structural zero elements, because Stan doesn’t give you a good way of setting up structural zeroes in a covariance matrix (Sean Pinkney or Ben Goodrich might know how to do this?).
This makes me curious how the lavaan package in R manages this. There’s a Bayesian version of lavaan built on top of Stan, blavaan. The first example right at the top of the home pages for both the lavaan and blavaan is Bollen’s democracy model. I guess it’s like the Scottish lip cancer data set for spatial modeling or Fisher’s iris data for regressions.
My questions
Consider a simple diagram among measurements like the following.
x <-> y y <-> z
This says there can be non-zero correlation between A/B and also between B/C, but the correlation between A/C is zero. It’s a simplified case of the box we saw in the actual example. These arrows implies the correlation matrix looks as follows.
| 1 rho[x,y] 0 | | rho[x,y] 1 rho[y,z] | = Omega | 0 rho[y,z] 1 |
Given that the correlation matrix Omega must be positive definite, this limits the range of rho[x,y] and rho[y,z]. For example, we can’t have rho[x,y] = rho[y,z] = 0.9, or rho[x,z] would have to be greater than zero to maintain positive definiteness.
Q1: Why doesn’t SEM instead say that the correlation rho[x,z] is just the minimum value it can be given rho[x,y] and rho[y,z]? I’m suggesting that we instead treat the above diagram as implying no additional correlation between x and z other than that implied by the correlation between x and y and the correlation between y and z? That is, why try to shrink rho[x,z] all the way to zero? From the text, it feels like the motivation is to enforce zero correlation in the model. But all this is doing is simplifying regressions—it won’t actually enforce zero correlation among the measurements that are modeled with zero correlation. I wished I’d asked Ken this question at dinner, but I’ll ping him about this blog post and hopefully get a response.
Of course, in the pragmatic Bayesian workflow, we’d use posterior predictive checks to evaluate whether there’s unmodeled correlation between x and z.
Q2: I’m also curious what Andrew and others think about enforcing structural zeroes in correlation between measurements as opposed to just estimating a dense covariance matrix and inspecting where the correlations fall.
One term to look up Bob is “Heywood case”. I have only seen examples for negative variance estimates, see https://stats.stackexchange.com/a/656803/1036, but it would not surprise me if the estimated correlation matrix is not feasible as well in some circumstances.
I am not up to date on the auto-identify what the path model could be. But Bollen’s work on “vanishing tetrads” kind of fits where you are going with Q number 2 (observed correlations are very inconsistent with the specified SEM).
Can you take a full covariance matrix and multiply it by a mask, pytorch style to block the gradients? Of course, then you have these free parameters independent of the data dragging on the sampler, but an appropriate prior at least stops them from becoming space filling. Or else, can you just take a vector of parameters and manually fill them into an otherwise zeroed out lower triangular matrix which is then taken to be the cholesky factorization?
> can you just take a vector of parameters and manually fill them into an otherwise zeroed out lower triangular matrix which is then taken to be the cholesky factorization?
I’ve tried this approach myself with mixed success. I say “mixed” because (a) I’m still not sure I implemented it well and (b) I still ran into major issues with the sampler getting stuck in bad parts of the parameter space. Compounding (b), I found it hard to specify a prior on the entries of the Cholesky factor that wouldn’t lead to asymmetric constraints on their values (ie, the ordering of variables in the matrix changed the posterior modes).
I’m not saying it can’t be done, of course, just that I had a hard time myself!
The problem is how to keep the resulting matrix positive definite. Masking a positive definite matrix with zeroes does not guarantee positive definiteness.
Bob:
I haven’t been following the details, but for marginal posterior mode estimation, we use this multivariate zero-avoiding prior so as to guarantee a positive definite estimate of the covariance matrix: https://sites.stat.columbia.edu/gelman/research/published/chung_cov_matrices.pdf
> But all this is doing is simplifying regressions—it won’t actually enforce zero correlation among the measurements that are modeled with zero correlation.
Doesn’t that sound like saying that setting sigma=1 for a normally distributed variable won’t actually enforce unit variance in measurements that are modeled as normal(mu,1)? Maybe I’m misunderstanding everything but I don’t think that a model will ever enforce anything on observations.
Yes, this is all I’m saying—nothing deep here.
I probably should’ve focused on what we can do in practice. I’d suggest posterior predictive checks on correlations of the data. If the data’s correlated, but the model assumes zero correlation, then posterior predictive inference will generate data without correlations that won’t match the actual data.
Those correlation parameters in question would be viewed as correlated residuals. The latent variables are also contributing to the correlations between observed variables. That is, the observed correlations between y1-y8 would primarily be explained by the latent variables DEM60 and DEM65. Then the y1 y5 correlation (and others) specifies that the residuals of those two observed variables are correlated after we account for the latent variables. And a residual correlation being fixed to 0 is saying that two observed variables are uncorrelated conditional on the latent variables.
This strikes me as similar to a simple regression model where we specify that the residuals are uncorrelated. And nonzero residual correlations remind me of the structured covariance matrices sometimes used in GLMM, as summarized here:
https://documentation.sas.com/doc/en/statug/15.2/statug_glimmix_syntax24.htm#statug.glimmix.gmxtabcovstructexamples
Thanks, Ed. This is getting at what I’m asking. I was asking people at the conference, and from what I could tell from the books and documentation and tutorials, the intention really was to set covariance to zero when there’s no double-headed arrow linking variables. I just wanted to make sure that was the intention versus just saying that there’s no additional correlation beyond what’s implied from other variables.
In the a <-> b and b <-> c case, insisting on zero correlation between an and c limits what the correlations between a/b and b/c can be. I was really expecting the minimum correlation between a and c given the other correlations to be the target.
Bob:
Regarding Q2 at the end of your post (“about enforcing structural zeroes in correlation between measurements as opposed to just estimating a dense covariance matrix and inspecting where the correlations fall”), my answer is that enforcing structural zeroes can be thought of as a quick-and-dirty form of prior distribution. You’re fitting a restricted model in settings where an unregularized estimate of the full covariance matrix would be unacceptably noisy. Or you’re fitting a restricted model because it is easier to interpret.
Right. This is one of those cases where the boundary between prior and likelihood feels blurry.
In a simpler example, imagine I’m fitting a multivariate normal to data. That’s a likelihood and I can take flat priors on mean and covariance. Or I can write the likelihood as a product of univariate normals and that’s a different likelihood I can also fit with flat priors on mean and covariance. Alternatively, I can stick to the multivariate normal and make an assumption that the covariances are all zero, which feels like a very strong prior on covariance.
I was expecting you to say that we should fit the full model and see if the correlation really is zero in the data. That’s what I recall you saying when we discussed priors that would restrict autoregressive priors to stationary parameterizations.
There is always a more general model (until there is no model). Instead of fitting a multivariate normal to data you could fit a copula of generalized gaussians and see if it’s really a multivariate gaussian.
SEM’s power comes from how it translates a causal path diagram to a covariance matrix for all variables, meaning you can model each observation (vector) as multivariate normal. This is why people have been estimating huge, complicated models for decades. Once you translate a path model to a covariance matrix, you just apply the same optimization and inference algorithms. This is also why anything beyond simple paths between continuous variables (e.g. ordered data; interactions between latent variables) ends up being a much more complicated problem. The former is like Bayesian inference with conjugate priors (easy to estimate); the latter with non-conjugate priors.
You could update your Stan model in a couple ways. You could construct the full covariance matrix and model X and Y simultaneously. (As Ed notes, the model doesn’t include a residual covariance between y2 and y8 but still implies a non-zero correlation through their shared variance from Ind60 and Dem60). If you had any missing data, you’d just need to select the observed rows/columns (FIML in SEM lingo; Marginal Maximum Likelihood in IRT; factorized likelihood elsewhere). You could also marginalize out the latent variables complete rather than estimate them in the model.
You can enforce a zero correlation in this simple example using the Cholesky approach. The basic idea is that you work backwards to the independent sources of variation, from which you can create a sparse or constrained covariance matrix. I believe this is what blavaan does under the hood in some situations. Here’s an example.
Stan Playground
Thanks, Simon!
Where I’m getting confused with the intent of the notation (not the statistics) is when there’s a double-headed arrow between A,B and another double-headed arrow between B,C, but no such arrow between A,C. That is consistent with zero correlation between A,C, but if that’s zero, then it limits how much correlation there can be in the A,B and B,C cases.
Indeed missing data’s very easy for multivariate normal models—just select the relevant sub-covariance matrix and you’re done. I’ll have to ask Ben Goodrich how blavaan works—he’s a genius with Cholesky factors (he’s the one who designed the LKJ prior on Cholesky factors for Stan and also our base encodings). I’m much more used to seeing sparse precision matrices and sparse Cholesky factors for those, which arise with spatio-temporal models.
To hint at an answer to Q1, there is a considerable literature on identification of SEMs that essentially hinges on whether elements of a covariance matrix are exactly zero, in the sense that if you were given the true covariance matrix with this or that element being exactly zero, could you solve uniquely for the structural parameters of interest using algebra?
No true Bayesian would think like that and you can get a unique posterior distribution in other ways, which Stan can often draw from efficiently.
Also, the old literature basically uses the observed covariance matrix like a sufficient statistic in the normal case, rather than passing the observed data. I did a Stan presentation in 2018 that corresponds to how Ken’s model would look like in other software that SEM people use that later became part of blavaan:
https://github.com/bgoodri/LERSIL
Ben says, “Also, the old literature basically uses the observed covariance matrix like a sufficient statistic in the normal case, rather than passing the observed data.”
Very clever. I’m always surprised at the normal distribution tricks that people get up to! I tried something like this with Andrew when I was trying to show Stan’s optimization could fit large-scale linear models faster than lm() built into R, but he didn’t like that I set the error scale to 1 for the fit, then took the sample residual variance for sigma in generated quantities. I think he thought it felt like cheating!
Posts like this keep me up at night. As a management scholar, many papers use SEM. However, it’s painfully clear that the discussion would baffle them. “SEM gave us the answer. They must be right.” Sometimes the diagnostics are clear that results are questionable if not invalid. Has anyone created a how-to guide to debunking (or verifying) SEM results. How many potentially high impact studies are actually wrong?
p.s. having started with LISREL (yes, I;m old) I *am* happy that SEM can be powerful.
TIA!
The second half of Bollen’s 1989 book is all about verifying SEM results. Arguably, the whole book is just workflow for SEM. As such, it has all sorts of great insights about working with data that don’t depend on SEM per se.
I think the biggest questions arise from the Bayesian approach to SEM (especially when moving a bit more into IRT-land with hierarchical priors), not SEM in general. Since the literature tends to focus on models identified by fixing some parameters to exact values (usually 0 or 1), there’s a disconnect between theoretical and statistical models. Furthermore, the applied researcher may have some domain expertise where they have a strong prior about some subset of elements of the covariance matrix but not others and it can be difficult to formalize this in a multivariate prior.
> Rather than independent errors, a SEM model explicitly indicates with arrows which pairs of observations are allowed to have non-zero correlation in the covariance matrix for the observations.
I found you description difficult to understand as you seem to mix “errors” with “observations” in your text. I think you meant “pairs of observations are allowed to have non-zero correlation in the covariance matrix for measurement errors (the difference between the value observed and the value predicted by the model from the value estimated for the latent variable).”
> Q1: Why doesn’t SEM instead say that the correlation rho[x,z] is just the minimum value it can be given rho[x,y] and rho[y,z]?
Because “It is likely that the error in measuring an indicator in 1960 is correlated with the error in measuring the same indicator in 1965. In addition the 1960 measures of y_2 and y_4 and the 1965 ones of y_6 and y_8 are from the same data source. The corresponding errors of measurement might be positively correlated due to systematic biases present in the source. Therefore the (4,2), (5,1), (6,2), (7,3), (8,4), and (8,6) off-diagonal elements of SIGMA_epsilon may be nonzero.” But for other pairs of indicators there is no reason to assume that there is any correlation in the errors. The thing that SEM does is not imposing that some errors are uncorrelated but exactly the opposite: it’s relaxing the usual independence assumption and allowing for correlations in measurement errors where appropriate introducing structure in the covariance matrix.
> I’m suggesting that we instead treat the above diagram as implying no additional correlation between x and z other than that implied by the correlation between x and y and the correlation between y and z?
I thing that you may be looking at this from the wrong angle. That should not happen in this model and if that puts an upper limit on those (joint) correlations that’s a feature and not a bug. Say that x = measurement error of A in 1960, y = measurement error of A in 1965 and z = measurement error of B in 1965. The measurement errors are disturbances that disrupt the relation between the latent and observed variables. There may be a systematic error in “A” measurements that introduces a correlation between x and y. There may be a systematic error in “1965” measurements that introduces a correlation between y and z. The model allows for that, but also considers that there is nothing in common between measuring A in 1960 and measuring B in 1965. Intuitively it seems that if the “A or B” effect on measurement errors is much larger than the “1960 or 1965” effect on measurement errors then the correlation between x and y will be very high but the correlation between y and z won’t.
I make a lot of mistakes like this, but this time I really did mean correlations on errors. It’s confusing because you can bounce back and forth, so I may not have been consistent.
We have something like x = some-regression + epsilon, where the error term epsilon ~ normal(0, Sigma) gets a multivariate normal distribution. Alternatively, you can get rid of the epsilon and write that as the values of your being correlated, x ~ normal(some-regression, Sigma).
I understand what SEM does and how it expands on independent errors. What I was trying to ask with Q1 is why they take the correlation to be zero, rather than the minimum correlation implied by other correlations taken to be non-zero. That is, why not the minimum correlation implied by positive-definiteness rather than zero? I’m still not sure they don’t because I find the language very hard to read.
> Rather than independent errors, a SEM model explicitly indicates with arrows which pairs of observations are allowed to have non-zero correlation in the covariance matrix for the observations.
> I make a lot of mistakes like this, but this time I really did mean correlations on errors.
Not sure if you are agreeing with me, but I do agree with you. I think that you really meant correlations on errors when you wrote “observations” at the end of that sentence.
> We have something like x = some-regression + epsilon, where the error term epsilon ~ normal(0, Sigma) gets a multivariate normal distribution. Alternatively, you can get rid of the epsilon and write that as the values of your being correlated, x ~ normal(some-regression, Sigma).
I don’t understand what you’re writing. Let’s say that there are five countries and the true latent values are DEM60 = [1 2 3 4 5] and DEM65 = [2 3 3 3 4]. Then we make three very simple measurements x = DEM60 + error_x, y = DEM60 + error_y and z = DEM65 + error_z. Let’s say that the three errors are uncorrelated and have standard deviation 0.01 and the observations look like x = [1.00 2.00 3.01 3.99 5.00], y = [1.02 2.00 3.01 4.00 5.00], z= [2.01 1.99 3.00 4.00 4.03]
The correlation between x and y is 1.00 and the correlation between x or y and z is 0.95. The correlation between the observations is not zero!
What would be zero is the correlation between the errors because they were assumed independent. (The realised correlation here is not zero but it would be close to zero if we had many countries and it will be zero on expectation.)
> What I was trying to ask with Q1 is why they take the correlation to be zero, rather than the minimum correlation implied by other correlations taken to be non-zero.
Because the minimum correlation implied by the other correlations has to be zero. You cannot set whatever value you want for the correlation between x and y and the correlation between y and z. They cannot be both equal to one. The correlation between errors will stay zero unless there is a relationship between them.
In the example above the model allows for some correlation between x and y, maybe there was a big systematic error in the measurements for the first country so the first component of the error_x and error_y is 1 instead of normal(0,0.01). Then they will have correlation close to 1. Let’s introduce now some correlation between y and z, maybe those measurements are made by the same person who really loves the fifth country and gives it one extra point so the last component of error_y and error_z is 1 instead of normal(0,0.01). Now the correlation between error_x and error_y is down to 0.6, the correlation between error_y and error_z is also 0.6 and the correlation between error_x and error_z is still 0. It cannot be affected indirectly by the systematic errors we introduced.
To be fair in they example above I could have said that the person gives an extra point to the first country rather than the last and then the correlation would have been 1 between the three error_x, error_y and error_z. But that would be more “luck” (or rather a careful choice) than a systematic effect. If there was a random effect affecting the measurement of x and y and a random effect affecting the measurement of y and z they would be in principle independent. Imagine a crossed hierarchical model where for each observation there may be multiple additive sources of variance but if two observations do not share any of those random effects they will be uncorrelated.
> The three industrialization observations are assumed to have zero correlation—there are no arrows between any of the three measurements in the SEM diagram.
I was curious to see what the correlation between the measurements actually looked like and it seems that the three industrialization observations are the most correlated of them all.
corrplot::corrplot(cor(cSEM::PoliticalDemocracy))
https://imgur.com/a/KrNwLFa
1. You asked what (b)lavaan do to handle structural zeroes in the covariance matrices. As of a few years ago, in practical terms, you must manually set the values for the elements of the matrix. In other words, you can build the entire matrix item-wise and place spike priors on individual elements or try to treat them as fixed (you can omit them!). This doesn’t work very well in my experience. Priors on the entire matrix tend to perform much better, but then some models will simply not be “identified” in the typical SEM way of defining the term or it may be nearly impossible to incorporate domain expertise into the development of priors (since the distribution is not flexible enough).
2. The biggest issue with blavaan relates to this exact issue: it is possible to specify priors on individual elements of the covariance matrix that jointly make a positive definite matrix (almost) impossible. This can mean the priors you specify may have unexpected, strange consequences: a prior that at first may seem “weak” may end up being incredibly informative because of these other constraints, since seemingly independent priors are actually related by virtue of the positive definite constraint. If you’re “testing a model”, then this is quite problematic as your priors along may not imply a positive definite matrix, so it’s not clear what model you are “testing.” The software will reject non-PD “solutions” and keep only PD solutions that cannot readily be anticipated. See here: https://d-nb.info/1304520641/34 (p.231+ gets it across).
I have essentially given up on classic SEM in Bayesian modeling for this reason. Other software (MPlus has some sort of Bayesian SEM I believe) has a similar issue. Some IRT stuff has similar issues as well.