With `bayesnec`

we have included a function that allows bootstrapped comparisons of posterior predictions. The main focus here is to showcase how the user can fit several different `bnec`

model fits and can compare differences in the posterior predictions across these fits for individual endpoint estimates (e.g. *NEC*, *NSEC* or *ECx*) or across a range of predictor values. Below we demonstrate usage of `compare_posterior`

for objects of class `bayesnecfit`

and `bayesmanecfit`

. In this example we compare different types of models and model sets using a single dataset. However, the intent of this function is to allow comparison across different datasets that might represent, for example, different levels of a fixed factor covariate. At this time `bnec`

does not allow inclusion of an interaction with a fixed factor. Including an interaction term within each of the non-linear models implemented in `bayesnec`

is relatively straightforward, and may be introduced in future releases. However, in many cases the functional form of the response may change with different levels of a given factor. The substantial complexity of defining all possible non-linear model combinations at each factor level means it unlikely this could be feasibly implemented in `bayesnec`

in the short term. In the meantime the greatest flexibility in the functional form of individual model fits can be readily obtained using models fitted independently to data within each factor level.

To run this vignette, we will also need some additional packages

```
library(brms)
library(bayesnec)
data(nec_data)
make_ecx_data <- function(top, bot, ec50, beta, x) {
top + (bot - top) / (1 + exp((ec50 - x) * exp(beta)))
}
x <- seq(0, 10, length = 12)
y <- make_ecx_data(x = x, bot = 0, top = 1, beta = 0.5, ec50 = 5)
set.seed(10)
dat <- data.frame(x = rep(x, 15), y = rnorm(15 * length(x), y, 0.2))
# Fit a set of models
exmp <- bnec(y ~ crf(x, model = "decline"), data = dat, iter = 2e3, refresh = 0)
```

```
class(exmp)
#> [1] "bayesmanecfit" "bnecfit"
```

This call fits all models that are suitable for modelling Gaussian-distributed response data, excluding all of the hormesis models, which we are not considering here. We have reduced the iterations back to the `brms`

default to make the example lighter, but we do not recommend this in practice. Now that we have our example fit, we can use this to make different model sets, purely for the purposes of illustrating `compare_posterior`

in `bayesnec`

and highlighting the rich information that is contained in the Bayesian posterior draws.

We can pull out the **nec** models and the **ecx** models separately, to create two alternative model fits of this data that we can compare to each other, as well as the original **all** model fit.

Now we have three different averaged model fits, all of class `bayemanec`

in this case (because they all contain multiple fits). We can compare their posterior estimates of the ecx10 values using `compare_posterior`

.

```
post_comp <- compare_posterior(list("all" = exmp, "ecx" = exmp_ecx,
"nec" = exmp_nec),
comparison = "ecx", ecx_val = 10)
names(post_comp)
#> [1] "posterior_list" "posterior_data" "diff_list"
#> [4] "diff_data" "prob_diff"
```

The `compare_posterior`

function outputs several elements in a named list. This includes the **posterior_data** for each model in the comparison as a `data.frame`

which we can use to plot a `geom_density`

plot of the posterior estimates, so they can be compared visually.

```
ggplot(data = post_comp$posterior_data, mapping = aes(x = value)) +
geom_density(mapping = aes(group = model, colour = model, fill = model),
alpha = 0.3) +
theme_classic()
```

From this you can see that the *EC _{10}* estimates are very similar for the

`wi`

in `exmp$mod_stats`

. The The `data.frame`

`diff_data`

can be used to make a similar plot, but specifically for the differences among models.

```
ggplot(data = post_comp$diff_data, mapping = aes(x = diff)) +
geom_density(mapping = aes(group = comparison, colour = comparison,
fill = comparison), alpha = 0.3) +
theme_classic()
```

This shows the differences among the three estimates. There is no difference in the **ecx** and **all** estimates (the probability density overlaps zero - red shaded curve). This is because for this example the simulated data are from a smooth **ecx** type curve and thus those models have high weight. The **nec** based *EC _{10}* estimates are greater than the

We can formally test the probability that the endpoint estimate for one model set is greater than the other using posterior differencing. This is contained in the compare posterior output as `prob_diff`

. Here you can see there is 48% chance that **all** is greater than **ecx**. There is only a 33% chance that **ecx** is greater than **nec**, and a 67% chance that the **nec** is greater than **ecx**.

```
post_comp$prob_diff
#> comparison prob
#> 1 all-ecx 0.483125
#> 2 all-nec 0.332500
#> 3 ecx-nec 0.329375
```

The user can also compare posterior fitted values across the full range of x values, using `comparison = "fitted"`

.

```
post_comp_fitted <- compare_posterior(list("all" = exmp, "ecx" = exmp_ecx,
"nec" = exmp_nec),
comparison = "fitted")
```

In the case of `comparison = "fitted"`

most of the elements returned by `compare_posterior`

are class `data.frame`

, with summary values for the posteriors, difference values and probabilities returned for each value of the predictor, for each model or model comparison.

```
head(post_comp_fitted$posterior_data)
#> model x Estimate Q2.5 Q97.5
#> 1 all 0.09090909 0.9560272 0.9079987 1.005893
#> 2 all 0.29313544 0.9559913 0.9079906 1.005071
#> 3 all 0.49536178 0.9559812 0.9079895 1.004756
#> 4 all 0.69758813 0.9559264 0.9079880 1.004655
#> 5 all 0.89981447 0.9558704 0.9079857 1.004118
#> 6 all 1.10204082 0.9557849 0.9079824 1.003546
head(post_comp_fitted$diff_data)
#> comparison diff.Estimate diff.Q2.5 diff.Q97.5
#> 1 all-ecx 0.0005514382 -0.06892429 0.06687044
#> 2 all-ecx 0.0004998985 -0.06877963 0.06680648
#> 3 all-ecx 0.0004271837 -0.06881629 0.06668127
#> 4 all-ecx 0.0004301296 -0.06887058 0.06666006
#> 5 all-ecx 0.0002705520 -0.06895098 0.06624540
#> 6 all-ecx 0.0002613768 -0.06843418 0.06591832
#> x
#> 1 0.09090909
#> 2 0.29313544
#> 3 0.49536178
#> 4 0.69758813
#> 5 0.89981447
#> 6 1.10204082
```

Using the collated posterior_data we can plot the predicted curves with confidence bounds for each of the input models. This shows clearly that the **ecx** model set begins to decline earlier than the **nec** set, which is flat prior to the **nec** step point, and then declines more rapidly.

```
ggplot(data = post_comp_fitted$posterior_data) +
geom_line(mapping = aes(x = x, y = Estimate, color = model), size = 0.5) +
geom_ribbon(mapping = aes(x = x, ymin = Q2.5, ymax = Q97.5, fill = model),
alpha = 0.3)+
labs(x = "Concentration", y = "Posterior estimate") +
theme_classic()
```

We can plot the differences between pairs of models in the list by plotting `diff.Estimate`

from `diff_data`

and using colours for the different comparisons. This plot highlights where the differences among these model sets are the greatest and explains the differences we see above in the *EC _{10}* endpoint estimates. The

```
ggplot(data = post_comp_fitted$diff_data) +
geom_line(mapping = aes(x = x, y = diff.Estimate, color = comparison),
size = 0.5) +
geom_ribbon(mapping = aes(x = x, ymin = diff.Q2.5, ymax = diff.Q97.5,
fill = comparison), alpha = 0.3) +
labs(x = "Concentration", y = "Posterior difference") +
theme_classic()
```

And finally we can plot the probability that one model is greater than the other by plotting `prob`

from `diff_data`

. The pattern of this plot is identical to the plot of differences, but the y axis now shows the probability of these differences. The red line hovers around 0.5 clearly indicating the lack of significant difference in the **ecx** and **all** model sets at any point of the predictor axis. The green and blue curves pass through 0.5 at several points, meaning there are parts of the curve where there is no significant difference between the **nec** and the **ecx** or **all** model set predictions. The greatest probability of difference among these curves is between values of ~3 and ~4 of the predictor, where the probability of difference trends towards low values (**ecx** is higher than **nec**), and at around 5, where the probability difference is very high (**nec** is higher than **ecx**). These are the two points where the **ecx** set deviates most from the **nec** set.

```
ggplot(data = post_comp_fitted$prob_diff) +
geom_line(mapping = aes(x = x, y = prob, color = comparison), size = 0.5) +
labs(x = "Concentration", y = "Probability of difference") +
theme_classic()
```

Here we use the difference in the **ecx** and **nec** models sets purely to showcase the functionality of `compare_posterior`

and highlight some of the strengths of using a Bayesian approach. These differences are expected in this case, because the simulated data we used are from a smooth **exc** function. The **nec** models are fitting a break-point model that has a single mean value up to the **nec**, followed by a relatively sharp decline. When fitted to these smooth curves the **nec** models do not really fit the data well and produce inflated (less conservative) estimates of the *EC _{10}* value. Fortunately, the model averaging approach solves this issue by yielding a very low weight to the

The `compare_posterior`

function was actually developed with the intention of comparing model sets based on different data to assess statistically the difference in either the endpoint estimates, or across the entire fitted curve, amongst factors of interest. Examples would include assessing for changes in toxicity across different WET tests through time, or comparing different levels of a factor treatment variable (such as climate change scenarios), to look at how these interact with concentration to impact toxicity.