Hi! I’m a PhD student trying to do a phylogenetic variance decomposition for the first time - it was meant to be a quick extra analysis for a paper but, predictably, is becoming a bit of a headache.
This is an analysis of simulated data, where we simulated removal processes on social networks, and we’re interested in the effect of things we’ve simulated (removal in this case) vs things inherent to the original data we used (e.g species). The variables I’ll mention are
-
removal method: three different ways of simulating removal
-
species: the species of the network
-
identity: the specific network
-
cov_matrix: a phylogenetic covariance matrix (1s on the diagonal)
We are attempting to model multiplication factor (MF), a positive value which shows the amount of change in our metric from the original network until after removal (1 = no change in the metric). We are using log(MF) in the models, in order to de-emphasise the scale difference between 0 to 1 and 1+. This also helps it approximate normality.
Each individual network has three measurements of MF, one for each removal method. Each network ID belongs to one species, with some species having multiple network IDs and some having just one. These variable groupings may be part of what is causing me trouble!
Here is a toy example (there is more data than this). The only exception is where, for various reasons, the simulation failed and therefore some network IDs do not have all three removal methods (60 of 263).
This is the most simple form, with other variables of interest that I may want to add in later. For now, the model formula looks like this:
clust_formula1 <- brmsformula(log(MF) ~ 1 + (1|removal_method) +
(1|gr(phylo, cov = cov_matrix)) +
(1|species/identity))
clust_model1 <- brm(clust_formula1,
data=filter(MF_clust1, !is.infinite(MF)),
data2 = list(cov_matrix = cov_matrix),
family = student(),
chains=1,
iter = 4000,
prior = c(
prior(normal(0, 2), "Intercept"),
prior(student_t(3, 0, 2.5), "sd"),
prior(student_t(3, 0, 2.5), "sigma")
),
control = list("adapt_delta" = 0.9,
"max_treedepth" = 12))
Just using 1 chain for now while I try to diagnose issues, as using more chains takes a long time while I’m messing with things. While my changes to the control parameters have successfully reduced the number of divergent transitions, there are still a small number. The mcmc_pairs output for this model run with 2 chains looks like this:
and pp_check like this:
The student distribution does a much better job than Gaussian at approximating the peakiness of the data, but does lead to some very extreme outliers. These are not biologically impossible, but very unlikely.
While trying to diagnose this model, I trimmed it down to just
clust_formula0 <- brmsformula(log(MF) ~ 1 + (1|removal_method))
and faced the same issues with divergent transitions. Run with the same parameters as model 1 above, here are the mcmc_pairs plots:
Hard to say as much with such relatively small numbers of divergent transitions, but they seem to pop up at the 0 boundaries of the sd parameters. I’ve tried changing control parameters, using different model families, and have removed some variables from the model that had correlation with other terms.
Would love to hear people’s thoughts, as I’m very new to this form of modelling, and I might have made some pretty obvious mistakes with grouping terms etc…