Convergence issue with phylogenetic correlation structure model

[Renato Morais]

Hello all (Paul),

I am having convergence issues trying to fit a model in brms with a phylogenetic correlation structure. This was created by the cor_fixed function from my correlation matrix and inserted in the autocor argument. The data has multiple categorical and numeric traits, as well as numeric environmental variables and a method variable. It also has multiple observations per species (intraspecific observations are included in the phylogeny with 0 distance). Numeric or categorical traits vary only between species, but environmental variables can vary both between and within species. Please, see the attached html from Markdown for a better glimpse into data structure and workflow. I ran 8 chains with 5000 steps each, 2500 of which were warmup, and with a thinning of 2. 

fmod <- brmsformula (log10 (Kmax) ~ csstmean + clog2MaxSizeTL + cFormFactor + clog2pelnpp + 
                            Diet + Movimentation + Method)

#mod <- brm (fmod, data = db, family = gaussian (),
#                   autocor = cor_fixed (corBrown),
#                   prior = c (prior (normal (0, 3), class = 'b'),
#                              prior (normal (0, 3), class = 'Intercept')),
#                   chains = 8, iter = 5000, thin = 2,  
#                   control = list(adapt_delta = 0.85),
#                   refresh = 50)

The final model converged nicely for the completely species-unrelated variables (environmental and method) but did pretty bad relative to the trait (species-related) variables. 

## Warning: The model has not converged (some Rhats are > 1.1). Do not analyse the results! 
## We recommend running more iterations and/or setting stronger priors.

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: log10(Kmax) ~ csstmean + clog2MaxSizeTL + cFormFactor + clog2pelnpp + Diet + Movimentation + Method 
##    Data: db (Number of observations: 1922) 
## Samples: 8 chains, each with iter = 5000; warmup = 2500; thin = 2; 
##          total post-warmup samples = 10000
##     ICs: LOO = NA; WAIC = NA; R2 = NA
##  
## Population-Level Effects: 
##                     Estimate Est.Error l-95% CI u-95% CI Eff.Sample  Rhat
## Intercept              -0.85      0.59    -1.80     0.04          4 23.24
## csstmean                0.00      0.00     0.00     0.00       9555  1.00
## clog2MaxSizeTL         -0.38      0.01    -0.41    -0.36          8  1.70
## cFormFactor            -0.54      0.48    -1.62     0.10          4 46.91
## clog2pelnpp             0.01      0.00     0.01     0.01       9961  1.00
## DietHerDet             -0.03      0.10    -0.29     0.13          8  2.30
## DietHerMac              0.18      0.25    -0.61     0.63          7  2.73
## DietInvMob             -0.06      0.05    -0.19     0.03          6  2.51
## DietInvSes              0.07      0.11    -0.16     0.32          5  3.28
## DietOmnivr              0.04      0.10    -0.20     0.20          7  2.36
## DietPlktiv              0.07      0.11    -0.20     0.24          6  2.80
## MovimentationBnthDw     0.08      0.07    -0.04     0.25          5  3.38
## MovimentationBtPlAs    -0.03      0.19    -0.21     0.48          4  5.09
## MovimentationBtPlDw     0.14      0.11    -0.01     0.43          5  4.19
## MovimentationPelgAs     0.15      0.17    -0.07     0.60          5  4.12
## MovimentationPelgDw     0.03      0.24    -0.36     0.67          5  4.62
## MethodMarkRc           -0.10      0.00    -0.10    -0.10      10000  1.00
## MethodOthRin           -0.20      0.00    -0.20    -0.20      10000  1.00
## MethodOtolth           -0.15      0.00    -0.15    -0.15      10000  1.00
## MethodScalRi           -0.14      0.00    -0.14    -0.14      10000  1.00
## MethodUnknow           -0.12      0.00    -0.12    -0.12       9800  1.00
## 
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, and Rhat is the potential 
## scale reduction factor on split chains (at convergence, Rhat = 1).

I could rerun the model with less chains or more iterations but this model took me almost 60 hours to run with the 8 chains in parallel and I am unsure whether running longer chains will achieve convergence or if the problem is something completely unrelated to that. Some weird things are happening in the pairs plot (attached) but I guess with so little effectives it's probably not telling much?

I would really appreciate any insights on how to achieve convergence.

All the best, thank you

Renato Morais

[Paul Buerkner]

Is there a specific reason you chose to use cor_fixed() instead of the 'cov_ranef' argument of brm?

If you have not already done so, I recommend vignette("brms_phylogenetics") as a starting point for how to (usually) fit these models in brms.

--
You received this message because you are subscribed to the Google Groups "brms-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to brms-users+unsubscribe@googlegroups.com.
To post to this group, send email to brms-...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/brms-users/586f7304-b2d8-4019-82f7-54800ce7d06f%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Renato Morais]

Hi Paul, thanks for your answer.

I started with the vignette but actually yes, there is a reason for using cor_fixed.

This was done to avoid including Species as a group-level effect. Including this group effect 'kills' the model because there is no variability left in traits variables (i.e. it makes no sense to look at the relationship between the response and some of these variables per species). So I was using the cor_fixed function to include the correlation structure without having to relate it to a group-level effect.

Is this not a valid procedure? Am I missing something important?

--

Renato Morais

PhD candidate in Marine Science

College of Science and Engineering

ARC Centre of Excellence for Coral Reef Studies

James Cook University

Townsville, Australia

http://thebellwoodreeffishlab.com

[Paul Buerkner]

I might be mistaken but I don't think your reason to be valid. Is this a common procedure in your field?

Both the covariates (constant within species as I understand) as well as the group-level intercept can coexist within one model, with the latter taking the excees variation not explained by your covariates. Unless there is something very special going on with you model, I don't think using cor_fixed() and opposed to cov_ranef is the most reasonable approach.

[Renato Morais]

It's not uncommon and it works alright with a frequentist GLS.

https://lukejharmon.github.io/ilhabela/instruction/2015/07/03/PGLS/

I don't know if it's the best implementation in brms and maybe that is causing the problem though...

[Paul Buerkner]

In general, if their are multiple observations per species, one should apply some grouping structure to account for this dependency.
At the most basic level this is a varying intercept (1 | species). If one then goes ahead and wants to model (phylogenetic) correlations between

species, one models the varying intercepts as correlated according to the pre-defined correlation matrix via cov_ranef.

Using cor_fixed for this has two drawbacks:

(1) Its slower because the correlation matrix is larger (nobs > nspecies)

(2) One may run into numerical problems due to 1s in the off-diagnoal for two observations of the same species.

Thus, at least for use in brms, I recommend going with the multilevel approach.

I am not an expert in this field, but I feel this is the favorable approach in general, as it naturally generalizes from the appropriate model if we did not have any correlation between species

[Renato Morais]

Thank you for your clarifications. 

I understand that the grouping structure approach is more efficient numerically and, thus, computationally. But wouldn't employ a correlation matrix with all observations be able to account for the non-independence just as well as the grouping structure?

Anyway, I appreciate your thoughts on cor_fixed and the most adequate ways of doing. Thank you!

[Paul Buerkner]

I think it would, but the models are not equivalent, as the cor_fixed model is not a multilevel model. By setting the correlation of the two obs of the same specieis to 1, you are creating a lot of problems, which I suggest you should avoid.