`bayesnec`

There are a range of models available in `bayesnec`

and
the working `bnec`

function supports individual model
fitting, as well as multi-model fitting with Bayesian model
averaging.

The argument `model`

in a `bayesnecformula`

is
a character string indicating the name(s) of the desired model (see
`?models`

for more details, and the list of models
available). If a recognised model name is provided, a single model of
the specified type is fit, and `bnec`

returns a model object
of class `bayesnecfit`

. If a vector of two or more of the
available models are supplied, `bnec`

returns a model object
of class `bayesmanecfit`

containing Bayesian model averaged
predictions for the supplied models, providing they were successfully
fitted.

Model averaging is achieved through a weighted sample of each fitted
models’ posterior predictions, with weights derived using the
`loo_model_weights`

function from `loo`

(Vehtari et al. 2020; Vehtari, Gelman, and Gabry
2017). Individual `brms`

model fits can be extracted
from the `mod_fits`

element and can be examined
individually.

The `model`

may also be one of `"all"`

, meaning
all of the available models will be fit; `"ecx"`

meaning only
models excluding the \(\eta =
\text{NEC}\) step parameter will be fit; `"nec"`

meaning only models with a specific \(\eta =
\text{NEC}\) step parameter will be fit; `"bot_free"`

meaning only models without a `"bot"`

parameter (without a
bottom plateau) will be fit; `"zero_bounded"`

are models that
are bounded to be zero; or `"decline"`

excludes all hormesis
models, i.e., only allows a strict decline in response across the whole
predictor range (see below **Parameter definitions**).
There are a range of other pre-defined model groups available. The full
list of currently implemented model groups can be seen using:

```
library(bayesnec)
models()
#> $nec
#> [1] "nec3param" "nec4param" "nechorme" "nechorme4" "necsigm" "neclin"
#> [7] "neclinhorme" "nechormepwr" "nechorme4pwr" "nechormepwr01"
#>
#> $ecx
#> [1] "ecx4param" "ecxlin" "ecxexp" "ecxsigm" "ecxwb1" "ecxwb2"
#> [7] "ecxwb1p3" "ecxwb2p3" "ecxll5" "ecxll4" "ecxll3" "ecxhormebc4"
#> [13] "ecxhormebc5"
#>
#> $all
#> [1] "nec3param" "nec4param" "nechorme" "nechorme4" "necsigm" "neclin"
#> [7] "neclinhorme" "nechormepwr" "nechorme4pwr" "nechormepwr01" "ecxlin" "ecxexp"
#> [13] "ecxsigm" "ecx4param" "ecxwb1" "ecxwb2" "ecxwb1p3" "ecxwb2p3"
#> [19] "ecxll5" "ecxll4" "ecxll3" "ecxhormebc4" "ecxhormebc5"
#>
#> $bot_free
#> [1] "nec3param" "nechorme" "necsigm" "neclin" "neclinhorme" "nechormepwr"
#> [7] "ecxlin" "ecxexp" "ecxsigm" "ecxwb1p3" "ecxwb2p3" "ecxll3"
#> [13] "ecxhormebc4" "nechormepwr01"
#>
#> $zero_bounded
#> [1] "nec3param" "nechorme" "necsigm" "nechormepwr" "nechormepwr01" "ecxexp"
#> [7] "ecxsigm" "ecxwb1p3" "ecxwb2p3" "ecxll3" "ecxhormebc4"
#>
#> $decline
#> [1] "nec3param" "nec4param" "neclin" "ecxlin" "ecxexp" "ecxsigm" "ecx4param" "ecxwb1"
#> [9] "ecxwb2" "ecxwb1p3" "ecxwb2p3" "ecxll5" "ecxll4" "ecxll3"
#>
#> $hormesis
#> [1] "nechorme" "nechorme4" "neclinhorme" "nechormepwr" "nechorme4pwr" "nechormepwr01"
#> [7] "ecxhormebc4" "ecxhormebc5"
```

Where possible we have aimed for consistency in the interpretable
meaning of the individual parameters across models. Across the currently
implemented model set, models contain from two (basic linear or
exponential decay, see **ecxlin** or
**ecxexp**) to five possible parameters
(**nechorme4**), including:

\(\tau = \text{top}\), usually interpretable as either the y-intercept or the upper plateau representing the mean concentration of the response at zero concentration;

\(\eta = \text{NEC}\), the
No-Effect-Concentration value (the x concentration value where the
breakpoint in the regression is estimated at, see **Model types
for NEC and EC_{x} estimation** and
(Fox 2010) for more details on parameter
based

\(\beta = \text{beta}\), generally
the exponential decay rate of response, either from 0 concentration or
from the estimated \(\eta\) value, with
the exception of the **neclinhorme** model where it
represents a linear decay from \(\eta\)
because slope (\(\alpha\)) is required
for the linear increase;

\(\delta = \text{bottom}\), representing the lower plateau for the response at infinite concentration;

\(\alpha = \text{slope}\), the
linear decay rate in the models **neclin** and
**ecxlin**, or the linear increase rate prior to \(\eta\) for all hormesis models;

\(\omega\) = \(\text{EC\textsubscript{50}}\) notionally the 50% effect concentration but may be influenced by scaling and should therefore not be strictly interpreted, and

\(\epsilon = \text{d}\), the
exponent in the **ecxsigm** and **necisgm**
models.

\(\zeta = \text{f}\) A scaling
exponent exclusive to model **ecxll5**.

In addition to the model parameters, all
**nec**-containing models have a step function used to
define the breakpoint in the regression, which can be defined as

\[ f(x_i, \eta) = \begin{cases} 0, & x_i - \eta < 0 \\ 1, & x_i - \eta \geq 0 \\ \end{cases} \]

All models provide an estimate for the No-Effect-Concentration
(*NEC*). For model types with **nec** as a prefix,
the *NEC* is directly estimated as parameter \(\eta = \text{NEC}\) in the model, as per
(Fox 2010). Models with
**ecx** as a prefix are continuous curve models, typically
used for extracting *EC _{x}* values from
concentration-response data. In this instance the

`sig_val`

) percentage certainty (based on the Bayesian
posterior estimate) that the response falls below the estimated value of
the upper asymptote (\(\tau =
\text{top}\)) of the response (i.e. the response value is
significantly lower than that expected in the case of no exposure). The
default value for `sig_val`

is 0.01, which corresponds to an
alpha value (Type-I error rate) of 0.01 for a one-sided test of
significance. See `?nsec`

for more details. We currently
recommend only using the `"nec"`

model set for estimation of
*EC _{x}* estimates can be equally obtained from both

`"nec"`

and `"ecx"`

models.
`"ecx"`

models fitted to the same data as
`"nec"`

models (see the Comparing
posterior predictions) vignette for an example. However, we
recommend using `"all"`

models where `"nec"`

models can fit some
datasets better than `"ecx"`

models and the model averaging
approach will place the greatest weight for the outcome that best fits
the supplied data. This approach will yield There is ambiguity in the definition of *EC _{x}*
estimates from hormesis models—these allow an initial increase in the
response (see Mattson 2008) and include
models with the character string

`horme`

in their name—as
well as those that have no natural lower bound on the scale of the
response (models with the string `ecx`

function has arguments `hormesis_def`

and
`type`

, both character vectors indicating the desired
behaviour. For `hormesis_def = "max"`

,
`hormesis_def = "control"`

(the
default) indicates that `type = "relative"`

`type = "absolute"`

(the default)
`type = "direct"`

, a direct interpolation of the response on
the predictor is obtained.Models that have an exponential decay (most models with parameter
\(\beta = \text{beta}\)) with no \(\delta = \text{bottom}\) parameter are
0-bounded and are not suitable for the Gaussian family, or any family
modelled using a `"logit"`

or `"log"`

link because
they cannot generate predictions of negative response values.
Conversely, models with a linear decay (containing the string
**lin** in their name) are not suitable for modelling
families that are 0-bounded (Gamma, Poisson, Negative Binomial, Beta,
Binomial, Beta-Binomial) using an `"identity"`

link. These
restrictions do not need to be controlled by the user, as a call to
`bnec`

with `models = "all"`

in the formula will
simply exclude inappropriate models, albeit with a message.

Strictly speaking, models with a linear hormesis increase are not
suitable for modelling responses that are 0, 1-bounded (Binomial-, Beta-
and Beta-Binomial-distributed), however they are currently allowed in
`bayesnec`

, with a reasonable fit achieved through a
combination of the appropriate distribution being applied to the
response, and `bayesnec`

’s `make_inits`

function
which ensures initial values passed to `brms`

yield response
values within the range of the user-supplied response data.

The **ecxlin** model is a basic linear decay model,
given by the equation: \[y_i = \tau -
e^{\alpha} x_i\] with the respective `brmsformula`

being

```
#> y ~ top - exp(slope) * x
#> top ~ 1
#> slope ~ 1
```

Because the model contains linear predictors it is not suitable for
0, 1-bounded data (i.e. Binomial and Beta families with an
`"identity"`

link function). As the model includes a linear
decline with concentration, it is also not suitable for 0,
`Inf`

bounded data (Gamma, Poisson, Negative Binomial with an
`"identity"`

link).

The **ecxexp** model is a basic exponential decay model,
given by the equation: \[y_i = \tau
e^{-e^{\beta} x_i}\] with the respective `brmsformula`

being

```
#> y ~ top * exp(-exp(beta) * x)
#> top ~ 1
#> beta ~ 1
```

The model is 0-bounded, thus not suitable for Gaussian response data
or the use of a `"logit"`

or `"log"`

link
function.

The **ecxsigm** model is a simple sigmoidal decay model,
given by the equation: \[y_i = \tau
e^{-e^{\beta} x_i^{e^\epsilon}}\] with the respective
`brmsformula`

being

```
#> y ~ top * exp(-exp(beta) * x^exp(d))
#> d ~ 1
#> top ~ 1
#> beta ~ 1
```

The model is 0-bounded, thus not suitable for Gaussian response data
or the use of a `"logit"`

or `"log"`

link
function.

The **ecx4param** model is a 4-parameter sigmoidal decay
model, given by the equation: \[y_i = \tau +
(\delta - \tau)/(1 + e^{e^{\beta} (\omega - x_i)})\] with the
respective `brmsformula`

being

```
#> y ~ top + (bot - top)/(1 + exp((ec50 - x) * exp(beta)))
#> bot ~ 1
#> ec50 ~ 1
#> top ~ 1
#> beta ~ 1
```

The **ecxwb1** model is a 4-parameter sigmoidal decay
model which is a slight reformulation of the Weibull1 model of Ritz et al. (2016), given by the equation: \[y_i = \delta + (\tau - \delta) e^{-e^{e^{\beta}
(x_i - \omega)}}\] with the respective `brmsformula`

being

```
#> y ~ bot + (top - bot) * exp(-exp(exp(beta) * (x - ec50)))
#> bot ~ 1
#> ec50 ~ 1
#> top ~ 1
#> beta ~ 1
```

The **ecxwb1p3** model is a 3-parameter sigmoidal decay
model which is a slight reformulation of the Weibull1 model of Ritz et al. (2016), given by the equation: \[y_i = {0} + (\tau - {0}) e^{-e^{e^{\beta} (x_i -
\omega)}}\] with the respective `brmsformula`

being

```
#> y ~ 0 + (top - 0) * exp(-exp(exp(beta) * (x - ec50)))
#> ec50 ~ 1
#> top ~ 1
#> beta ~ 1
```

The model is 0-bounded, thus not suitable for Gaussian response data
or the use of a `"logit"`

or `"log"`

link
function.

The **ecxwb2** model is a 4-parameter sigmoidal decay
model which is a slight reformulation of the Weibull2 model of Ritz et al. (2016), given by the equation: \[y_i = \delta + (\tau - \delta) (1 -
e^{-e^{e^{\beta} (x_i - \omega)}})\] with the respective
`brmsformula`

being

```
#> y ~ bot + (top - bot) * (1 - exp(-exp(-exp(beta) * (x - ec50))))
#> bot ~ 1
#> ec50 ~ 1
#> top ~ 1
#> beta ~ 1
```

While very similar to the **ecxwb1** (according to Ritz et al. 2016), fitted
**ecxwb1** and **ecxwb2** models can differ
slightly.

The **ecxwb2p3** model is a 3-parameter sigmoidal decay
model, which is a slight reformulation of the Weibull2 model of Ritz et al. (2016), given by the equation: \[y_i = {0} + (\tau -{0}) (1 - e^{-e^{e^{\beta}
(x_i - \omega)}})\] with the respective `brmsformula`

being

```
#> y ~ 0 + (top - 0) * (1 - exp(-exp(-exp(beta) * (x - ec50))))
#> ec50 ~ 1
#> top ~ 1
#> beta ~ 1
```

While very similar to the **ecxwb1p3** (according to Ritz et al. 2016), fitted
**ecxwb1p3** and **ecxwb2p3** models can
differ slightly. The model is 0-bounded, thus not suitable for Gaussian
response data or the use of a logit or log link function.

The **ecxll5** model is a 5-parameter sigmoidal
log-logistic decay model, which is a slight reformulation of the LL.5
model of Ritz et al. (2016), given by the
equation: \[y_i = \delta + (\tau - \delta) /
(1 + e^{-e^{\beta} (x_i - \omega)})^{e^\zeta}\] with the
respective `brmsformula`

being

```
#> y ~ bot + (top - bot)/(1 + exp(exp(beta) * (x - ec50)))^exp(f)
#> bot ~ 1
#> top ~ 1
#> beta ~ 1
#> ec50 ~ 1
#> f ~ 1
```

The **ecxll4** model is a 4-parameter sigmoidal
log-logistic decay model which is a slight reformulation of the LL.4
model of Ritz et al. (2016), given by the
equation: \[y_i = \delta + (\tau - \delta)/
(1 + e^{e^{\beta} (x_i - \omega)})\] with the respective
`brmsformula`

being

```
#> y ~ bot + (top - bot)/(1 + exp(exp(beta) * (x - ec50)))
#> bot ~ 1
#> top ~ 1
#> beta ~ 1
#> ec50 ~ 1
```

The **ecxll3** model is a 3-parameter sigmoidal
log-logistic decay model, which is a slight reformulation of the LL.3
model of Ritz et al. (2016), given by the
equation: \[y_i = 0 + (\tau - 0)/ (1 +
e^{e^{\beta} (x_i - \omega)})\] with the respective
`brmsformula`

being

```
#> y ~ 0 + (top - 0)/(1 + exp(exp(beta) * (x - ec50)))
#> top ~ 1
#> beta ~ 1
#> ec50 ~ 1
```

`"logit"`

or `"log"`

link
function.

The **ecxhormebc5** model is a 5 parameter log-logistic
model modified to accommodate a non-linear hormesis at low
concentrations. It has been modified from to the “Brain-Cousens” (BC.5)
model of Ritz et al. (2016), given by the
equation: \[y_i = \delta + (\tau - \delta +
e^{\alpha} x)/ (1 + e^{e^{\beta} (x_i - \omega)})\] with the
respective `brmsformula`

being

```
#> y ~ bot + (top - bot + exp(slope) * x)/(1 + exp(exp(beta) * (x - ec50)))
#> bot ~ 1
#> top ~ 1
#> beta ~ 1
#> ec50 ~ 1
#> slope ~ 1
```

The **ecxhormebc4** model is a 5-parameter log-logistic
model similar to the **exchormebc5** model but with a lower
bound of 0, given by the equation: \[y_i = 0
+ (\tau - 0 + e^{\alpha} x)/ (1 + e^{e^{\beta} (x_i - \omega)})\]
with the respective `brmsformula`

being

```
#> y ~ 0 + (top - 0 + exp(slope) * x)/(1 + exp(exp(beta) * (x - ec50)))
#> top ~ 1
#> beta ~ 1
#> ec50 ~ 1
#> slope ~ 1
```

`"logit"`

or `"log"`

link
function.

The **neclin** model is a basic linear decay model
equivalent to **ecxlin** with the addition of the
*NEC* step function, given by the equation: \[y_i = \tau - e^{\alpha} \left(x_i - \eta \right)
f(x_i, \eta)\] with the respective `brmsformula`

being

```
#> y ~ top - exp(slope) * (x - nec) * step(x - nec)
#> top ~ 1
#> slope ~ 1
#> nec ~ 1
```

Because the model contains linear predictors it is not suitable for
0, 1-bounded data (Binomial and Beta distributions with
`"identity"`

link). As the model includes a linear decline
with concentration, it is also not suitable for 0, `Inf`

bounded data (Gamma, Poisson, Negative Binomial with
`"identity"`

link).

The **nec3param** model is a basic exponential decay
model equivalent to **ecxexp** with the addition of the
*NEC* step function, given by the equation: \[y_i = \tau e^{-e^{\beta} \left(x_i - \eta \right)
f(x_i, \eta)}\] with the respective `brmsformula`

being

```
#> y ~ top * exp(-exp(beta) * (x - nec) * step(x - nec))
#> top ~ 1
#> beta ~ 1
#> nec ~ 1
```

For Binomial-distributed response data in the case of
`"identity"`

link this model is equivalent to that in Fox (2010). The model is 0-bounded, thus not
suitable for Gaussian response data or the use of a `"logit"`

or `"log"`

link function.

The **nec4param** model is a 3-parameter decay model
with the addition of the *NEC* step function, given by the
equation: \[y_i = \delta + (\tau - \delta)
e^{-e^{\beta} \left(x_i - \eta \right) f(x_i, \eta)}\] with the
respective `brmsformula`

being

```
#> y ~ bot + (top - bot) * exp(-exp(beta) * (x - nec) * step(x - nec))
#> bot ~ 1
#> top ~ 1
#> beta ~ 1
#> nec ~ 1
```

The **nechorme** model is a basic exponential decay
model with an *NEC* step function equivalent to
**nec3param**, with the addition of a linear increase prior
to \(\eta\), given by the equation
\[y_i = (\tau + e^{\alpha} x_i) e^{-e^{\beta}
\left(x_i - \eta \right) f(x_i, \eta)}\] with the respective
`brmsformula`

being

```
#> y ~ (top + exp(slope) * x) * exp(-exp(beta) * (x - nec) * step(x - nec))
#> top ~ 1
#> beta ~ 1
#> nec ~ 1
#> slope ~ 1
```

The **nechorme** model is a *hormesis* model
(Mattson 2008), allowing an initial
increase in the response variable at concentrations below \(\eta\). The model is 0-bounded, thus not
suitable for Gaussian response data or the use of a `"logit"`

or `"log"`

link function. In this case the linear version
(**neclinhorme**) should be used.

The **nechormepwr** model is a basic exponential decay
model with an *NEC* step function equivalent to
**nec3param**, with the addition of a power increase prior
to \(\eta\), given by the equation:
\[y_i = (\tau + x_i^{1/(1+e^{\alpha})})
e^{-e^{\beta} \left(x_i - \eta \right) f(x_i, \eta)}\] with the
respective `brmsformula`

being

```
#> y ~ (top + x^(1/(1 + exp(slope)))) * exp(-exp(beta) * (x - nec) * step(x - nec))
#> top ~ 1
#> beta ~ 1
#> nec ~ 1
#> slope ~ 1
```

The **nechormepwr** model is a *hormesis* model
(Mattson 2008), allowing an initial
increase in the response variable at concentrations below \(\eta\). The model is 0-bounded, thus not
suitable for Gaussian response data or the use of a `"logit"`

or `"log"`

link function. Because the model can generate
predictions > 1 it should not be used for Binomial and Beta
distributions with `"identity"`

link. In this case the
**nechromepwr01** model should be used.

The **neclinhorme** model is a basic linear decay model
with an *NEC* step function equivalent to
**neclin**, with the addition of a linear increase prior to
\(\eta\), given by the equation: \[y_i = (\tau + e^{\alpha} x_i) - e^{\beta}
\left(x_i - \eta \right) f(x_i, \eta)\] with the respective
`brmsformula`

being.

```
#> y ~ (top + exp(slope) * x) - exp(beta) * (x - nec) * step(x - nec)
#> top ~ 1
#> beta ~ 1
#> nec ~ 1
#> slope ~ 1
```

The **neclinhorme** model is a *hormesis* model
(Mattson 2008), allowing an initial
increase in the response variable at concentrations below \(\eta\). This model contains linear
predictors and is not suitable for 0, 1-bounded data (Binomial and Beta
distributions with `"identity"`

link). As the model includes
a linear decline with concentration, it is also not suitable for 0,
`Inf`

bounded data (Gamma, Poisson, Negative Binomial with
`"identity"`

link).

The **nechorme4** model is 4 parameter decay model with
an *NEC* step function equivalent to **nec4param**
with the addition of a linear increase prior to \(\eta\), given by the equation: \[y_i = \delta + ((\tau + e^{\alpha} x_i) - \delta
) e^{-e^{\beta} \left(x_i - \eta \right) f(x_i, \eta)}\] with the
respective `brmsformula`

being

```
#> y ~ bot + ((top + exp(slope) * x) - bot) * exp(-exp(beta) * (x - nec) * step(x - nec))
#> bot ~ 1
#> top ~ 1
#> beta ~ 1
#> nec ~ 1
#> slope ~ 1
```

The **nechorme4** model is a *hormesis* model
(Mattson 2008), allowing an initial
increase in the response variable at concentrations below \(\eta\).

The **nechorme4pwr** model is 4 parameter decay model
with an *NEC* step function equivalent to
**nec4param** with the addition of a power increase prior
to \(\eta\), given by the equation:
\[y_i = \delta + ((\tau +
x_i^{1/(1+e^{\alpha})}) - \delta) e^{-e^{\beta} \left(x_i - \eta \right)
f(x_i, \eta)}\] with the respective `brmsformula`

being

```
#> y ~ bot + ((top + x^(1/(1 + exp(slope)))) - bot) * exp(-exp(beta) * (x - nec) * step(x - nec))
#> bot ~ 1
#> top ~ 1
#> beta ~ 1
#> nec ~ 1
#> slope ~ 1
```

The **nechorme4pwr** model is a *hormesis* model
(Mattson 2008), allowing an initial power
increase in the response variable at concentrations below \(\eta\). Because the model can generate
predictions > 1 it should not be used for Binomial and Beta
distributions with `"identity"`

link. In this case the
**nechromepwr01** model should be used.

The **nechormepwr01** model is a basic exponential decay
model with an *NEC* step function equivalent to
**nec3param**, with the addition of a power increase prior
to \(\eta\), given by the equation:
\[y_i = \left(\frac{1}{(1 +
((1/\tau)-1) e^{-e^{\alpha}x_i}}\right) e^{-e^{\beta} \left(x_i - \eta
\right) f(x_i, \eta)}\] with the respective
`brmsformula`

being

```
#> y ~ (1/(1 + ((1/top) - 1) * exp(-exp(slope) * x))) * exp(-exp(beta) * (x - nec) * step(x - nec))
#> top ~ 1
#> beta ~ 1
#> nec ~ 1
#> slope ~ 1
```

The **nechormepwr01** model is a *hormesis* model
(Mattson 2008), allowing an initial
increase in the response variable at concentrations below \(\eta\). The model is 0-bounded, thus not
suitable for Gaussian response data or the use of a `"logit"`

or `"log"`

link function. In this case the linear version
(**neclinhorme**) should be used.

The **necsigm** model is a basic exponential decay model
equivalent to **ecxlin** with the addition of the
*NEC* step function, given by the equation: \[y_i = \tau e^{-e^{\beta} ((x_i - \eta) f(x_i,
\eta))^{e^\epsilon}f(x_i, \eta)}\] with the respective
`brmsformula`

being

```
#> y ~ top * exp(-exp(beta) * (step(x - nec) * (x - nec))^exp(d) * step(x - nec))
#> top ~ 1
#> beta ~ 1
#> nec ~ 1
#> d ~ 1
```

The model is 0-bounded, thus not suitable for Gaussian response data
or the use of a `"logit"`

or `"log"`

link
function. In addition, there may be theoretical issues with combining a
sigmoidal decay model with an *NEC* step function because where
there is an upper plateau in the data the location of \(\eta\) may become ambiguous. Estimation of
No-Effect-Concentrations using this model are not currently recommended
without further testing.

Fox, David R. 2010. “A Bayesian approach for
determining the no effect concentration and hazardous concentration in
ecotoxicology.” *Ecotoxicology and Environmental
Safety* 73 (2): 123–31.

Mattson, Mark P. 2008. “Hormesis Defined.” *Ageing
Research Reviews* 7 (1): 1–7. https://doi.org/https://doi.org/10.1016/j.arr.2007.08.007.

Ritz, Christian, Florent Baty, Jens C Streibig, and Daniel Gerhard.
2016. “Dose-Response Analysis Using R.”
*PLoS ONE* 10 (12): e0146021. https://doi.org/10.1371/journal.pone.0146021.

Vehtari, Aki, Jonah Gabry, Mans Magnusson, Yuling Yao, Paul-Christian
Bürkner, Topi Paananen, and Andrew Gelman. 2020. *Loo: Efficient
Leave-One-Out Cross-Validation and WAIC for Bayesian Models*. https://mc-stan.org/loo/.

Vehtari, Aki, Andrew Gelman, and Jonah Gabry. 2017. “Practical
Bayesian Model Evaluation Using Leave-One-Out Cross-Validation and
WAIC.” *Statistics and Computing* 27: 1413–32. https://doi.org/10.1007/s11222-016-9696-4.