I’m working with some models of binary outcome data where I want the probability parameter of a Bernoulli distribution to depend on several Gamma survival functions assessed at a specified threshold value. To achieve this, I’m using brms’ non-linear functionality, so that the model formula looks something like this (simplified):

```
bf(response ~ 1-gamma_cdf(exp(threshold), exp(alpha), exp(beta))),
alpha ~ condition,
beta ~ condition,
threshold ~ 1,
nl = TRUE,
family = bernoulli(link="identity"))
```

The `gamma_cdf`

function comes from Stan, since as far as I’m aware Stan has no `pgamma`

function that corresponds to the one in R.

The models are fitting the data more or less effectively, and there are clear differences between the data I’ve (manually) simulated from the posteriors of most of them that allow me to rule out some models via eye test. However, there are a couple of models that look to be fitting the data roughly equally well, so I’d like to use `loo`

to get a quantitative estimate of their relative performance—if only to confirm my judgement that the performance is fairly similar. Unfortunately when I try to do this (and also when I try to use `posterior_predict`

and `conditional_effects`

) I get an error message telling me:

```
Error in gamma_cdf(exp(threshold), exp(alpha), exp(beta)) :
could not find function "gamma_cdf"
Most likely this is because you used a Stan function in the non-linear model formula that is not defined in R. If this is a user-defined function, please run 'expose_functions(., vectorize = TRUE)' on your fitted model and try again.
```

I’m just wondering if anyone who’s more experienced with this sort of thing might have a suggestion regarding how I can circumvent this problem so I can perform the comparison I’m interested in. If so, I’d love to hear it!