# Posterior\_predict for non-linear models with truncated response family

**URL:** <https://discourse.mc-stan.org/t/posterior-predict-for-non-linear-models-with-truncated-response-family/19689>\
**Category:** Modeling\
**Tags:** posterior-predictive, brms\
**Created:** [December 8, 2020, 5:34am UTC](https://discourse.mc-stan.org/t/posterior-predict-for-non-linear-models-with-truncated-response-family/19689 "2020-12-08T05:34:22Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![jscamac](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jscamac/32/756_2.png) [@jscamac](https://discourse.mc-stan.org/u/jscamac)\
**Post date:** [December 8, 2020, 5:34am UTC](https://discourse.mc-stan.org/t/posterior-predict-for-non-linear-models-with-truncated-response-family/19689/1 "2020-12-08T05:34:22Z")

</div>

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

 ![Screen Shot 2020-12-08 at 4.11.52 pm](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/7/75bc2c6f82bfefc6c46258fd7352a3ef775cdb1e.jpeg)

---

<div class="post-metadata">

**Author:** ![jscamac](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jscamac/32/756_2.png) [@jscamac](https://discourse.mc-stan.org/u/jscamac)\
**Post date:** [December 8, 2020, 11:53am UTC](https://discourse.mc-stan.org/t/posterior-predict-for-non-linear-models-with-truncated-response-family/19689/2 "2020-12-08T11:53:11Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![martinmodrak](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/martinmodrak/32/133_2.png) [@martinmodrak](https://discourse.mc-stan.org/u/martinmodrak)\
**Post date:** [December 14, 2020, 9:32pm UTC](https://discourse.mc-stan.org/t/posterior-predict-for-non-linear-models-with-truncated-response-family/19689/3 "2020-12-14T21:32:53Z")

</div>

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!