Comparing posterior endpoint values

set.seed(333)
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(333)
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, 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. 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.

exmp_nec <- pull_out(exmp, model = "nec")
exmp_ecx <- pull_out(exmp, model = "ecx")

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"      "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₁₀ estimates are very similar for the ecx and all model sets. This is because the ecx model types dominate the model weights in this all fit, see wi in exmp$mod_stats. The EC₁₀ estimate is slightly lower (more conservative) for the ecx based models.

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₁₀ estimates are greater than the ecx based estimates in this case.

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 55% chance that all is greater than ecx. There is only a 22% chance that ecx is greater than NEC, and a 78% chance that the NEC is greater than ecx.

post_comp$prob_diff
#>   comparison      prob
#> 1    all-ecx 0.5545693
#> 2    all-nec 0.2944118
#> 3    ecx-nec 0.2221528

Comparing posterior fitted values

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 1.019064 0.9661288 1.081557
#> 2   all 0.29313544 1.018926 0.9661288 1.079875
#> 3   all 0.49536178 1.018739 0.9661288 1.078295
#> 4   all 0.69758813 1.018546 0.9661170 1.076359
#> 5   all 0.89981447 1.018291 0.9660999 1.074599
#> 6   all 1.10204082 1.017883 0.9660999 1.072197
head(post_comp_fitted$diff_data)
#>   comparison diff.Estimate   diff.Q2.5 diff.Q97.5          x
#> 1    all-ecx  -0.004558534 -0.08537598 0.07724148 0.09090909
#> 2    all-ecx  -0.004414021 -0.08404266 0.07642943 0.29313544
#> 3    all-ecx  -0.004368055 -0.08321144 0.07519140 0.49536178
#> 4    all-ecx  -0.004269343 -0.08075704 0.07393649 0.69758813
#> 5    all-ecx  -0.004246644 -0.07927299 0.07251096 0.89981447
#> 6    all-ecx  -0.004167774 -0.07736581 0.07094054 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),
            linewidth = 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₁₀ endpoint estimates. The ecx and all model sets are relatively similar across the entire range of concentrations (red band overlaps zero). The green band is the difference between NEC and all (specifically all-NEC) and the blue band is the difference between NEC and ecx (specifically ecx-NEC). You can see from the green and blue bands that the NEC set has slightly lower estimates than all and ecx at low to moderate predictor values; the NEC based curve then crosses the other two, predicting higher values of the response at around a concentration of 4-5; then drops rapidly, producing lower predictions than the other curves at around 5.

ggplot(data = post_comp_fitted$diff_data) +
  geom_line(mapping = aes(x = x, y = diff.Estimate, color = comparison),
            linewidth = 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),
            linewidth = 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₁₀ value. Fortunately, the model averaging approach solves this issue by yielding a very low weight to the NEC based models, and the all model set is indistinguishable from the ecx set, which contains the underlying true generating model in this case. As such, this is not really an interesting example—other than illustrating the behaviour of NEC model types in this setting.

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.