Hi all,
I’m playing around with fitting various non-linear growth models… and broadly everything is going fine. However, I’ve noted some issues when I use to model to make posterior predictions to new data. I’m almost certain these issues are more to do with my lack of understanding of what posterior_predict does under the hood.
Details of the model and issue below:
I’m attempting to model the growth of leaf_length over time as a function of the following non-linear growth curve called a von Bertalanffy:
\alpha \times (1 + ((\frac{\text{init leaf length}}{\alpha})^{0.33}-1) \times \exp(-R \times t)
In the above equation, init leaf length is a known quantity (leaf length at the start of the experiment). \alpha is the asymptotic leaf length, R is effectively the growth rate term and t is time.
Anyway, my aim is to estimate both \alpha and R, the latter as a function of various covariates and random effects (i.e. individual_id and site_id). In order to fit this function, I’ve needed to do a few hacks. The first is that I’ve needed to model R on the log scale to ensure leaf length, such that I’m enforcing growth to be either zero or positive. The other thing I’ve had to do is use a half-normal distribution for the response (hence the truncates(lb=0). I’ve tried other distributions (e.g. lognormal(link=“identity”)) where I shouldn’t need to truncate… but they produce dubious fits. The half-normal model fits fine (everything converges, no divergent iterations), and the posterior predictive checks appear to be OK (see below).
brms::brm(
bf(
leaf_length | trunc(lb=0) ~ Asym * (1 + ((init_leaf_length/Asym)^(0.33)-1) * exp(-exp(lograte) * julian_date))^3,
lograte ~ treatment_id * elevation_cs + treatment_id * aspect_id + elevation_cs * aspect_id +
(1|site_id) + (1|ind_id),
Asym ~ 1, nl = TRUE),
data = growth_model_data,
control = list(adapt_delta = 0.99, max_treedepth = 15),
cores = 4,
prior =
prior(normal(0,5),class = "b", nlpar ="lograte") +
prior(uniform(5,30), lb=5, ub=30, nlpar="Asym")
)
Anyway, when I take this model and predict to a new dataset:
new_data <- data.frame(key = rep(c("Control", "Gap"), each =100),
julian_date = rep(seq(0.001, 4, length.out = 100), times = 2),
init_leaf_length = rep(1,200),
treatment_id = rep(c(0,1), each =100),
aspect_id = rep(0,200),
elevation_cs = rep(0,200))
apply(brms::posterior_predict(model, newdata = new_data, re.form = NA), 2, mean)[1]
[1] 3.466569
The mean (or median) posterior estimates when t= 0.001 appear to be significantly higher than expected for a plant with an initial leaf length of 1 cm.
If I run the predictions with posterior_linpred I get the estimate I expect
apply(brms::posterior_linpred(model, newdata = new_data, re.form = NA), 2, mean)[1]
[1] 1.032739
If I rerun the predictions with a higher initial_leaf_length (e.g. from 1 to 10), the magnitude of the difference between the two is much reduced.
So to my question, is this purely an artefact of the model struggling to estimate/make predictions close to the distribution boundary (in this case zero)? I know the predicted uncertainty should be different between posterior_predict vs posterior_linpred as the former accounts for the residual error. But should the means be substantially different?
Any advice welcomed!
packageVersion("brms")
[1] ‘2.14.0’
James
NOTE using posterior_epred produces similar predictions to
