I'm attempting to use latent variable modelling, as described by Hui (2016) and using the boral package in R, to explore relationships between plant species composition and environmental variables while allowing for possible interactions between species. The species data consist of presence-absence of 165 taxa over 66 survey sites.
I'd like to estimate the overall proportion of co-variation between species accounted for by the environmental variables. Hui (2016) and Warton et al. (2015) describe using the trace of the residual covariance matrix for this purpose. They compare the trace value of a model with latent variables but no covariates to that from a model with latent variables plus covariates. One expects that the trace value will be lower for the second model, with the proportional reduction being an estimate of how much of the species co-occurrence is explained by covariates. The boral package provides the function
get.residual.cor which returns, among other things, the trace value.
Attempting to apply the suggested method to my data and models gives surprising results. Alternative models with environmental covariates included are returning trace values substantially higher (e.g. 30%) than that for a base model with only latent variables. This leads me to worry that the models are somehow badly mis-posed but, if that is the case, it isn't obvious. Fitted coefficients for species against environmental variables appear to make ecological sense and posterior predictions (e.g. for species richness of field survey sites) seem reasonable.
A boral model is, behind the scenes, a Bayesian model fitted using MCMC via JAGS. I understand (in a rudimentary way) that this means the residual covariance matrix should be treated with the same caveats as any other parameter sampled by MCMC. The package documentation for function
Of course, the trace itself is random due to the MCMC sampling, and so it is not always guaranteed to produce sensible answers!
So my question is how does one know when to trust the trace value? Are there any recommended post-hoc checks? My models appear to have converged (as judged by the Geweke's diagnostic returned by the boral function) and auto-correlation of parameter samples seems acceptable. Perhaps there is a posterior sample size requirement for the trace to be reliable? I'm presently running 100,000 MCMC iterations and thinning to every 25th sample.