Unexpected result with simulated data?

[miguel-porto]

Hello, I feel like this is an unexpected result given the simulated data I'm trying. Here's a totally reproducible example. Please see the comments within the script.

# Generate a species matrix where sp1 and sp2 are highly correlated, whereas sp3 is
# independent from the others
nsamples = 1000
sp1 = runif(nsamples, 0, 1)
sp = data.frame(cbind(
    "sp1" = ifelse(sp1>0.7, 1, 0)
    , "sp2" = ifelse(sp1 + rnorm(nsamples, 0, 0.05) > 0.7, 1, 0))
    , "sp3" = ifelse(runif(nsamples, 0, 1)>0.7, 1, 0))

cor(sp)
# sp1 and sp2 are highly correlated (~0.90)
# sp3 is independent from the others

table(sp[,1], sp[,2])
# Again confirm association between sp1 and sp2

# Make a unit-level random effect
random = data.frame(uni = factor(paste0("P",1:nsamples)))
rownames(random) = random[,1]
rownames(sp) = random[,1]

# Model data
md = as.HMSCdata(Y = as.matrix(sp), Random = random)

# Fit model
m = hmsc(md, family="probit", niter=5000, nburn=1000, thin=10)

# Plot mixing object
mix = as.mcmc(m, parameters = "paramLatent")
plot(mix)
# Huh? All the latent parameters are centered on zero...
# I would expect that sp1 and sp2, because of being highly correlated, would show similar
# latent parameters, different from zero.

apply(corRandomEff(m, cor = TRUE)[, , , 1], 1:2, mean)
# This confirms the unexpected result: the high correlation of sp1 and sp2 is not captured
# in the random effect! (or very slightly captured)
# It seems the model wasn't able to extract the latent factor to which sp1 and sp2
# respond in consonance...
# Moreover, these mean correlations considerably fluctuate between runs with the same
# data; e.g. sp1 x sp2 varied between 0.24 and 0.35 correlation

Am I doing something wrong in model fitting so that it's actually not converging? Pardon if this is not the right place to post this question... Thanks, all the best, Miguel