After some trawling through the literature, I found a solution (well, a pair of similar solutions), so I thought I would put an update here in case someone has the same problem as me somewhere in the future. I’m sure my explanation below is limited in it’s technical breadth, but hopefully it’s useful to someone.
The problem, as I’ve learned, is that the factor loadings are only identifiable down to \Lambda\Lambda^T, and so different MCMC chains can easily end up exploring different modes. Constraints on the values of \Lambda can help with identifiability - in particular, constraining the upper diagonal to zeroes forces the modes found to differ only by sign, whereas having no constraints means that the modes can also be rotations of one another. Additional constraints, such as constraining the diagonals to be positive, have been claimed to force a single signed mode, but for reasons far past my understanding can find non-trivial modes such as the ones above.
The solution(s) I have found involve post-processing the posterior draws in order to choose a single of the explored modes. The first solution (Erosheva and Curtis 2017, Dealing with reflection invariance in Bayesian factor analysis) uses the upper-diagonal=0 constraint, and returns a draws x dimension set of sign changes. The second solution (Papastamoulis and Ntzoufras 2020, [2004.05105] On the identifiability of Bayesian factor analytic models) uses no constraints on \Lambda and involves a varimax rotation step before calculating sign changes.
As I am only interested in the values of \Lambda insofar as \Lambda\Lambda^T, but do wish to interpret the values of the latent variables as an ordination, the first solution is more convenient for me as it only requires applying the sign changes to the latent variables, rather than also extracting and applying the same rotation to them. Applying this solution drops the Rhat values from somewhere in the region of 1.8 to 1 to 1.01.
I have a couple of simple follow-up questions that I haven’t quite been able to figure out that would help apply this post-processing step more easily in a workflow:
- Is it possible to modify the draws held within a CmdStanMCMC object (ideally in a copy)?
- Is it possible to extract the dimensions of \Lambda from a CmdStanMCMC object? I would rather not have to regex these out of the parameter names if possible!