Structural equation models without circular definition of latent variables

Question

I'm trying to perform Bayesian structural equation modeling in Python and PyMC3, but I think the problem is similar for most probabilistic progamming languages, include JAGS, Stan, etc.

SEMs are often defined using equations like these:

$\begin{equation} \boldsymbol{y}=\boldsymbol{\nu}+\boldsymbol{\Lambda}\boldsymbol{\eta}+\boldsymbol{\epsilon}\\ \boldsymbol{\eta}=\boldsymbol{\alpha}+\boldsymbol{B}\boldsymbol{\eta}+\boldsymbol{\zeta}\\ \boldsymbol{\epsilon}\sim\mathrm{MvNormal}\left(0, \boldsymbol{\Theta}\right)\\ \boldsymbol{\zeta}\sim\mathrm{MvNormal}\left(0, \boldsymbol{\Psi}\right) \end{equation}$

As you can see, the latent variables in $\boldsymbol{\eta}$ are defined in terms of themselves. The model makes sense because $\boldsymbol{B}$ is a sparse, invertible, non-symmetric matrix that defines paths between latent variables, but I'm not sure how to code this in practice. Can the second equation be re-written/expanded to so that the definition is no longer "circular" and therefore more amenable to coding?

I realize that under SEM one usually marginalizes out $\boldsymbol{\eta}$ and so the problem is moot. However, I think it should be possible to code it up with $\boldsymbol{\eta}$ explicitly sampled.

Answer 1

Assuming $(\mathbf{I} - \mathbf{B})$ is invertible, we can rewrite the $\boldsymbol{\eta}$ equation as follows:

$$ \begin{align} \boldsymbol{\eta} &= \boldsymbol{\alpha}+\mathbf{B}\boldsymbol{\eta}+\boldsymbol{\zeta} \\ &\Updownarrow \\ \boldsymbol{\eta} - \mathbf{B}\boldsymbol{\eta} &= \boldsymbol{\alpha}+\boldsymbol{\zeta} \\ &\Updownarrow \\ (\mathbf{I} - \mathbf{B})\boldsymbol{\eta} &= \boldsymbol{\alpha}+\boldsymbol{\zeta} \\ &\Updownarrow \\ \boldsymbol{\eta} &= (\mathbf{I} - \mathbf{B})^{-1} (\boldsymbol{\alpha}+\boldsymbol{\zeta}) \\ \end{align} $$

Where $\mathbf{I}$ is an identity matrix with the appropriate shape.