Posterior_predict for non-linear models with truncated response family

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).

    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!

[1] ‘2.14.0’


NOTE using posterior_epred produces similar predictions to

Ok I’ve progressed a little further here. I now realise that when I run posterior_linpred with transform = TRUE I effectively get the same as posterior_epred. I guess I was under the assumption that because I was using a gaussian I didn’t need to transform… but I guess the truncation does some form of transformation in the background. Still… the untransformed posterior_linpred estimates appear closer to the initial leaf length estimate… so I’m still somewhat confused.

Sorry for taking quite long to get to your question. It is relevant and well written.

posterior_epred gets the estimated mean while posterior_linpred gets the linear predictor. Now for non-truncated gaussian without any link the linear predictor is equal to the mean. But when truncation enters the picture this is no longer the case. This is simple to show in an extreme case: if the linear predictor was -100, the posterior mean of N(-100, \sigma) truncated at 0 has to be > 0 (and thus far from the linear predictor). In fact, the posterior mean will depend on both the linear predictor and \sigma.

On a higher level, I think what you observe is that the truncated normal is problematic for your model. Maybe a gamma distribution would be suitable (as this doesn’t make the variance grow so quickly with the mean as the lognormal does). There are likely other problems as you seem to have needed large adapt_delta and max_treedepth, but if that lets you fit without divergences / other diagnostics signalling a problem, it is likely OK).

Does that make sense?

And best of luck with your model!