I’m using a Gaussian distributional model to describe the distribution of a certain trait response in individual plant within species. I kept the structure and data quantity to a minimum for the reproducible example. Data consists in individual measures of the trait (called response in the repex dataset), associated with the species name (namely genus and species). Hypotheses of the model in the reproducible example :
- individual response within species distribution is Gaussian;
- we assume a nested random effect of species within genus on distribution mean ie
response; - we assume a random effect of species on distribution standard deviation
sigma.
Inference worked well. Then I predicted new data response based on species to get a sample of individual values for each unique species of the inference dataset, using the predictfunction with a simpler model formula (because in the actual analysis I need to predict distributional parameters without some of the random effects). From those predictions I computed the variances of the individual responses within species. I got a mismatch between those “predicted variances” and the ones that I directly computed from model estimates (namely computing exp(sigma_Intercept + ranef(fit)$’species’) with Estimatevalue).
It appears that using the option re_formula in predict or fitted functions alters the application of random effects on sigma: they are not taken into account when retrieving sigmas. Thus sigmas for all species are equal instead of being deviated from the common intercept according to the estimated random effect value. This behavior shows only on sigma (not on the response/mean). It shows even if the change in formula only affects the response/mean. The different cases are shown in the repex. I’m wondering if it is an actual bug in brms, and if I should open an issue on the GitHub repo ?
Reproductive example is attached: repex_bug_predict.R (4.9 KB)
Additional system/versions info:
- Operating System: Ubuntu 22.04.5 LTS x64
- R version: 4.5.2
- brms Version: 2.23.0
Note: The actual analysis is performed on way more data (more individual measurements, and more individuals belonging to other species/genus), it includes another hierarchical level (family) and additional sources of variability (related to sources & sites). Hence, I’m aware that the model proposed in the repex is not ecologically relevant, which doesn’t matter here.