Hi,
I am working with a dataset where I have measured traits on multiple plant species and populations. For most of the traits I measured, I was able to fit appropriate models and posterior predictive checks corroborated this - chains mixed well and converged (there were still a few divergent transition warnings, but I am not too worried about those).
However, there is one trait where I have a problem figuring out the right approach. It is a continuous positive number, technically bounded (as it is time until midpoint of growth). Unlike other traits I measure however, this turned out to be multimodal. This is all fine and good, as there are biological explanations for this.
My first instinct was to try fitting a distributional model, as upon splitting the data into species, the density plots looked unimodal enough.
prior_norm = c(
prior(normal(0, 1), "b"),
prior(student_t(3, 3.4, 2.5), "Intercept"),
prior(student_t(3, 0, 2.5), "sd"),
prior(student_t(3, 0, 2.5), class="Intercept", dpar="sigma"),
prior(student_t(3, 0, 2.5), class="sd", dpar="sigma"),
prior(normal(0, 1), class="b", dpar="sigma")
)
f1 <- brmsformula(xmid ~ elev_Org_m * treatment * delta_elev +
(1 | gr(phylo, cov=A)) + (1 | species:region),
sigma ~ elev_Org_m * treatment * delta_elev +
(1 | gr(phylo, cov=A)) + (1 | species:region),
family = lognormal)
#elev_Org_m and delta_elev describe the elevational origin of the species and plant
m <- brm(
f1,
data = d,
data2=list(A=A),
prior = prior_norm,
chains = 4, cores = 4,
iter = 4000, warmup = 1000
)
However, this is what the pp_check() looks like:
While the posterior does display bimodality, it doesn’t recover the actual structure of the data all that well. I have tried making sigma simpler (just the species level random effects) but there isn’t much of a difference.
So after some searching I found that mixture models could be an alternative approach, which could also make sense as populations can have slow and fast growers. As I cannot explicitly model this in a distreg framework, I thought of letting mixture models do the hard work for me
f2 <- brmsformula(xmid ~ elev_Org_m * treatment * delta_elev +
(1 | gr(phylo, cov=A)) + (1 | species:region), family = mixture(gaussian, gaussian))
#you will notice i use gaussian here: only did it as a starting point as maybe it makes things easier
#in general i use lognormal dist. to model my traits as all are restricted to be positive
prior_mix = c(
prior(normal(0, 5), "b", dpar="mu1"),
prior(normal(0, 5), "b", dpar="mu2"),
prior(normal(18, .05), "Intercept", dpar="mu1", lb = (5), ub = (70)),
prior(normal(40, .05), "Intercept", dpar="mu2", lb = (5), ub = (70)),
prior(student_t(3, 0, 14), "sd", dpar = "mu1"),
prior(student_t(3, 0, 14), "sd", dpar = "mu2"),
prior(student_t(3, 0, 14), class="sigma1"),
prior(student_t(3, 0, 14), class="sigma2"),
prior(dirichlet(10), class="theta")
)
m <- brm(
f2,
data = d,
data2=list(A=A),
prior = prior_mix,
chains = 4, cores = 4,
iter = 4000, warmup = 1000
)
I read some of the other posts on convergence issues on mixture models, and that’s why my priors for the intercept have such low standard deviations. At least at this stage, parts of the model have good convergence, but other parts don’t. I have very few divergent transitions, high Rhat, and low ESS values.
Looks like I’m a bit closer here. I still think lognormal is the appropriate distribution, but since the log scale is smaller, further reducing SDs for the intercept priors make estimation unnecessarily hard…
I guess essentially this question boils down to mixture models and appropriate model (and prior specification) for them, so maybe I’ve provided too much needless info here (sorry!). So my final question is where to go from here and how to specify the mixture model a bit better and get my models to work?
Thank you in advance, I appreciate any help.
- Operating System: MacOS (Apple Silicon)
- brms Version: 2.20.4
