Nothing
Priors"
In bayesnec: A Bayesian No-Effect- Concentration (NEC) Algorithm
Overview
The main focus here is to explain the default priors used in bayesnec and to showcase how the user can interrogate the priors used in a bayesec model and alternatively specify their own priors, should they wish to. This might be needed depending on the model and the data because bayesnec tries to find reasonable and yet only weakly informative priors for the model parameters by default. First we describe the default priors used in bayesnec and then follow up with a demonstration of how the user can specify priors in multiple ways for objects of class bayesnecfit and bayesmanecfit.
Default priors in bayesnec
The default priors used in bayesnec can generally be considered "weakly informative". They are constructed for each parameter of each model being fitted based on the characteristics of either the input response and predictor data, depending on which is relevant to the specific parameter scaling. In the case of parameters that scale with the response, priors are constructed based on the relevant link scaling, whether that be identity or the default (or user specific) link function for a specific distribution. The priors are constructed by bnec internally based on the chosen model, the distribution (including the relevant link function), the predictor and the response.
Priors for response-scaled parameters
Only the parameters top and bot scale specifically with the response. For Gaussian-distributed response data (or any response data for which the link ensures valid values of the response can take from +Inf to -Inf, including log and logit) priors are normal with a standard deviation of 2.5 times the standard deivation of the response data, and a mean set at the 90^th^ and 10^th^ quantiles for top and bot respectively. For Poisson-, Negative Binomial- and Gamma-distributed data the response is bounded by 0 and thus priors are Gamma, with a mean scaled to correspond to the 75^th^ and 25^th^ quantiles for top and bot respectively. The mean is linked mathematically to the shape (s) and rate parameters (r) by the equation $mean = s \times (1/r)$ with the Gamma shape parameter set at 2. For the Binomial, Beta, and Beta Binomial families estimates for top and bot must necessarily be constrained between 0 and 1 when modelled on the "identity" link (the default in bayesnec). Because of this constraint there is no need to adjust scaling based on the response. bayesnec uses beta(5, 2) and beta(2, 5) priors to provide a broad density centred across the upper and lower 0 to 1 range for the top and bot parameters, respectively.
Priors for predictor-scaled parameters
The parameters nec and ec50 scale with respect to the predictor because both of these are parameters in concentration-response curves the are estimated in units of concentration. To stabilise model fitting the nec and ec50 parameters are bounded to the upper and lower observed range in the predictor, under the assumption that the range of concentrations in the experiment were sufficient to cover the full range of the response outcomes. The priors used reflect the characteristics of the observed data that are used to guess the appropriate distribution. If the predictor data are bounded to 0 and >1, a Gamma prior is used, with maximum density (mean, see above) at the median value of the predictor, and a shape parameter of 5. If the predictor data are bounded to 0 and 1 a beta(2, 2) prior is used. For predictor data ranging from +Inf to -Inf a normal prior is used, with a mean set at the median of the predictor values and a standard deviation of 10 times the standard deviation of the predictor values.
Priors for other parameters
For the parameters beta, slope, d and f, we first ensured any relevant transformations in the model formula such that theoretical values of -Inf to +Inf are allowable, and a normal(0, 5) prior is used. For example, in the nec3param model, beta is an exponential decay parameter which must by definition be bounded to 0 and +Inf. Calling exp(beta) in the model formula ensures the exponent meets these requirements. Note also that a mean of 0 and standard deviation of 5 represents a relatively broad prior on the exponential scaling. See the Model details vignette or ?model("all") for more information on all the models available in bayesnec and their specific formulation.
User specified priors
There may be situations were the default bayesnec priors to not behave as desired, or the user wants to provide informative priors. For example the default priors may be too informative, yielding unreasonably tight confidence bands, although this is only likely where there are few data. Conversely, priors may be too vague, leading to poor model convergence, or an inability of bayesnec to find appropriate starting values. Alternatively, as indicated in the example below, the default priors may be of the wrong statistical distribution if there was insufficient information in the provided data for bayesnec to guess correctly the appropriate one to use.
The priors used in the default model fit can be extracted using pull_prior, and a sample or plot of prior values can be obtained from the individual brms model fits through the function sample_priors which samples directly from the prior element in the brm model fit. We can also use the function check_prior (based on the hypothesis function of brms) to assess how the posterior probability density for each parameter differs from that of the prior.
bayesnecfit
To set specified priors, it is simplest to start by letting bnec find the priors on its own, i.e. by not specifying the brm argument prior at all.
library(brms)
library(bayesnec)
data(nec_data)
# a single model
set.seed(333)
exmp_a <- bnec(y ~ crf(x, model = "nec3param"), data = nec_data,
family = Beta(link = "identity"),
iter = 1e4, control = list(adapt_delta = 0.99),
open_progress = FALSE)
class(exmp_a)
#> [1] "bayesnecfit" "bnecfit"
We can view the prior and posterior probability densities of all the parameters in the model using the function check_prior, based on the hypothesis function of brms. This can be useful to assess if priors are suitably vague, and/or if they might be having an undesirable influence on the posterior.
check_priors(exmp_a)
plot of chunk exmp3-checkpriors
In this case the priors seem reasonably vague, however there will be times when it is necessary to modify these priors. The user can take advantage of the function pull_prior to inspect what bnec came up with on its own, and decide how best to modify those priors to be more desirable.
pull_prior(exmp_a)
#> [[1]]
#> prior class coef group resp dpar nlpar lb ub source
#> normal(0, 5) b beta user
#> normal(0, 5) b Intercept beta (vectorized)
#> gamma(5, 2.28313180499098) b nec 0.03234801324009 3.22051966293556 user
#> gamma(5, 2.28313180499098) b Intercept nec 0.03234801324009 3.22051966293556 (vectorized)
#> beta(5, 2) b top 0 1 user
#> beta(5, 2) b Intercept top 0 1 (vectorized)
#> gamma(0.01, 0.01) phi 0 default
bnec chose a gamma prior on the NEC parameter of nec3param because the predictor nec_data$x is non-zero positive. However, imagine that in theory the predictor could have had negative values, it just happened to not have in this particular dataset. So let's go ahead and specify something else, say a normal with larger variance.
set.seed(333)
my_prior <- c(prior_string("beta(5, 1)", nlpar = "top"),
prior_string("normal(1.3, 2.7)", nlpar = "nec"),
prior_string("gamma(0.5, 2)", nlpar = "beta"))
exmp_b <- bnec(y ~ crf(x, model = "nec3param"), data = nec_data,
family = Beta(link = "identity"), prior = my_prior,
iter = 1e4, control = list(adapt_delta = 0.99),
open_progress = FALSE)
Two things are of note. If the user is specifying their own priors, bnec requires them to specify priors for all parameters. The pull_prior function shows the priors after the model was fitted, but suppose the user does not know what parameters were comprised in a particular model. In those instances, the user can call the function show_params(model = "all") to inspect the parameters of each function, or some targeted function in particular.
show_params(model = "nec3param")
#> [[1]]
#> y ~ top * exp(-exp(beta) * (x - nec) * step(x - nec))
#> top ~ 1
#> beta ~ 1
#> nec ~ 1
bayesmanecfit
The user can also specify a named list of priors when one or more models are being fitted to the same dataset.
set.seed(333)
my_priors <- list(nec3param = c(prior_string("beta(5, 1)", nlpar = "top"),
prior_string("normal(1.3, 2.7)", nlpar = "nec"),
prior_string("gamma(0.5, 2)", nlpar = "beta")),
nec4param = c(prior_string("beta(5, 1)", nlpar = "top"),
prior_string("normal(1.3, 2.7)", nlpar = "nec"),
prior_string("gamma(0.5, 2)", nlpar = "beta"),
prior_string("beta(1, 5)", nlpar = "bot")))
exmp_c <- bnec(y ~ crf(x, model = c("nec3param", "nec4param")), data = nec_data,
family = Beta(link = "identity"), prior = my_priors,
iter = 1e4, control = list(adapt_delta = 0.99),
open_progress = FALSE)
pull_prior also works for an object of class bayesmanecfit, as does check_priors which allows an option of passing a filename to save the prior and posterior probability density plots to a pdf.
pull_prior(exmp_c)
check_priors(exmp_c, filename = "Check_priors")
#> $nec3param
#> prior class coef group resp dpar nlpar lb ub source
#> gamma(0.5, 2) b beta user
#> gamma(0.5, 2) b Intercept beta (vectorized)
#> normal(1.3, 2.7) b nec user
#> normal(1.3, 2.7) b Intercept nec (vectorized)
#> beta(5, 1) b top user
#> beta(5, 1) b Intercept top (vectorized)
#> gamma(0.01, 0.01) phi 0 default
#>
#> $nec4param
#> prior class coef group resp dpar nlpar lb ub source
#> gamma(0.5, 2) b beta user
#> gamma(0.5, 2) b Intercept beta (vectorized)
#> beta(1, 5) b bot user
#> beta(1, 5) b Intercept bot (vectorized)
#> normal(1.3, 2.7) b nec user
#> normal(1.3, 2.7) b Intercept nec (vectorized)
#> beta(5, 1) b top user
#> beta(5, 1) b Intercept top (vectorized)
#> gamma(0.01, 0.01) phi 0 default
The user can also specify priors for one model only out of the entire set, bayesnec will return a message stating that it is searching for priors on its own when they are either ill-formed (e.g. incomplete or have a typo), or the user simply decided not to specify priors for a particular model, e.g.
set.seed(333)
my_priors <- list(nec3param = c(prior_string("beta(5, 1)", nlpar = "top"),
prior_string("normal(1.3, 2.7)", nlpar = "nec"),
prior_string("gamma(0.5, 2)", nlpar = "beta")),
nec4param = c(prior_string("beta(5, 1)", nlpar = "top"),
prior_string("normal(1.3, 2.7)", nlpar = "nec"),
prior_string("gamma(0.5, 2)", nlpar = "beta"),
prior_string("beta(1, 5)", nlpar = "bot")))
exmp_d <- bnec(y ~ crf(x, model = c("nec3param", "nec4param")), data = nec_data,
family = Beta(link = "identity"), prior = my_priors[1],
iter = 1e4, control = list(adapt_delta = 0.99),
open_progress = FALSE)
prior = my_priors[[1]] would also have worked because the argument priors can either take a brmsprior object directly, or a named list containing model-specific brmsprior objects.
Finally the user can also extend an existing bayesmanecfit object with the function amend, also by specifying custom-built priors.
set.seed(333)
necsigm_priors <- c(prior_string("beta(5, 1)", nlpar = "top"),
prior_string("gamma(2, 6.5)", nlpar = "beta"),
prior_string("normal(1.3, 2.7)", nlpar = "nec"),
prior_string("normal(0, 2)", nlpar = "d"))
exmp_e <- amend(exmp_d, add = "necsigm", priors = necsigm_priors)
Try the bayesnec package in your browser
Any scripts or data that you put into this service are public.
bayesnec documentation built on Sept. 11, 2024, 6:15 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Overview
The main focus here is to explain the default priors used in bayesnec and to showcase how the user can interrogate the priors used in a bayesec model and alternatively specify their own priors, should they wish to. This might be needed depending on the model and the data because bayesnec tries to find reasonable and yet only weakly informative priors for the model parameters by default. First we describe the default priors used in bayesnec and then follow up with a demonstration of how the user can specify priors in multiple ways for objects of class bayesnecfit and bayesmanecfit.
Default priors in bayesnec
The default priors used in bayesnec can generally be considered "weakly informative". They are constructed for each parameter of each model being fitted based on the characteristics of either the input response and predictor data, depending on which is relevant to the specific parameter scaling. In the case of parameters that scale with the response, priors are constructed based on the relevant link scaling, whether that be identity or the default (or user specific) link function for a specific distribution. The priors are constructed by bnec internally based on the chosen model, the distribution (including the relevant link function), the predictor and the response.
Priors for response-scaled parameters
Only the parameters top and bot scale specifically with the response. For Gaussian-distributed response data (or any response data for which the link ensures valid values of the response can take from +Inf to -Inf, including log and logit) priors are normal with a standard deviation of 2.5 times the standard deivation of the response data, and a mean set at the 90^th^ and 10^th^ quantiles for top and bot respectively. For Poisson-, Negative Binomial- and Gamma-distributed data the response is bounded by 0 and thus priors are Gamma, with a mean scaled to correspond to the 75^th^ and 25^th^ quantiles for top and bot respectively. The mean is linked mathematically to the shape (s) and rate parameters (r) by the equation $mean = s \times (1/r)$ with the Gamma shape parameter set at 2. For the Binomial, Beta, and Beta Binomial families estimates for top and bot must necessarily be constrained between 0 and 1 when modelled on the "identity" link (the default in bayesnec). Because of this constraint there is no need to adjust scaling based on the response. bayesnec uses beta(5, 2) and beta(2, 5) priors to provide a broad density centred across the upper and lower 0 to 1 range for the top and bot parameters, respectively.
Priors for predictor-scaled parameters
The parameters nec and ec50 scale with respect to the predictor because both of these are parameters in concentration-response curves the are estimated in units of concentration. To stabilise model fitting the nec and ec50 parameters are bounded to the upper and lower observed range in the predictor, under the assumption that the range of concentrations in the experiment were sufficient to cover the full range of the response outcomes. The priors used reflect the characteristics of the observed data that are used to guess the appropriate distribution. If the predictor data are bounded to 0 and >1, a Gamma prior is used, with maximum density (mean, see above) at the median value of the predictor, and a shape parameter of 5. If the predictor data are bounded to 0 and 1 a beta(2, 2) prior is used. For predictor data ranging from +Inf to -Inf a normal prior is used, with a mean set at the median of the predictor values and a standard deviation of 10 times the standard deviation of the predictor values.
Priors for other parameters
For the parameters beta, slope, d and f, we first ensured any relevant transformations in the model formula such that theoretical values of -Inf to +Inf are allowable, and a normal(0, 5) prior is used. For example, in the nec3param model, beta is an exponential decay parameter which must by definition be bounded to 0 and +Inf. Calling exp(beta) in the model formula ensures the exponent meets these requirements. Note also that a mean of 0 and standard deviation of 5 represents a relatively broad prior on the exponential scaling. See the Model details vignette or ?model("all") for more information on all the models available in bayesnec and their specific formulation.
User specified priors
There may be situations were the default bayesnec priors to not behave as desired, or the user wants to provide informative priors. For example the default priors may be too informative, yielding unreasonably tight confidence bands, although this is only likely where there are few data. Conversely, priors may be too vague, leading to poor model convergence, or an inability of bayesnec to find appropriate starting values. Alternatively, as indicated in the example below, the default priors may be of the wrong statistical distribution if there was insufficient information in the provided data for bayesnec to guess correctly the appropriate one to use.
The priors used in the default model fit can be extracted using pull_prior, and a sample or plot of prior values can be obtained from the individual brms model fits through the function sample_priors which samples directly from the prior element in the brm model fit. We can also use the function check_prior (based on the hypothesis function of brms) to assess how the posterior probability density for each parameter differs from that of the prior.
bayesnecfit
To set specified priors, it is simplest to start by letting bnec find the priors on its own, i.e. by not specifying the brm argument prior at all.
library(brms) library(bayesnec) data(nec_data) # a single model set.seed(333) exmp_a <- bnec(y ~ crf(x, model = "nec3param"), data = nec_data, family = Beta(link = "identity"), iter = 1e4, control = list(adapt_delta = 0.99), open_progress = FALSE)
class(exmp_a) #> [1] "bayesnecfit" "bnecfit"
We can view the prior and posterior probability densities of all the parameters in the model using the function check_prior, based on the hypothesis function of brms. This can be useful to assess if priors are suitably vague, and/or if they might be having an undesirable influence on the posterior.
check_priors(exmp_a)
plot of chunk exmp3-checkpriors
In this case the priors seem reasonably vague, however there will be times when it is necessary to modify these priors. The user can take advantage of the function pull_prior to inspect what bnec came up with on its own, and decide how best to modify those priors to be more desirable.
pull_prior(exmp_a) #> [[1]] #> prior class coef group resp dpar nlpar lb ub source #> normal(0, 5) b beta user #> normal(0, 5) b Intercept beta (vectorized) #> gamma(5, 2.28313180499098) b nec 0.03234801324009 3.22051966293556 user #> gamma(5, 2.28313180499098) b Intercept nec 0.03234801324009 3.22051966293556 (vectorized) #> beta(5, 2) b top 0 1 user #> beta(5, 2) b Intercept top 0 1 (vectorized) #> gamma(0.01, 0.01) phi 0 default
bnec chose a gamma prior on the NEC parameter of nec3param because the predictor nec_data$x is non-zero positive. However, imagine that in theory the predictor could have had negative values, it just happened to not have in this particular dataset. So let's go ahead and specify something else, say a normal with larger variance.
set.seed(333) my_prior <- c(prior_string("beta(5, 1)", nlpar = "top"), prior_string("normal(1.3, 2.7)", nlpar = "nec"), prior_string("gamma(0.5, 2)", nlpar = "beta")) exmp_b <- bnec(y ~ crf(x, model = "nec3param"), data = nec_data, family = Beta(link = "identity"), prior = my_prior, iter = 1e4, control = list(adapt_delta = 0.99), open_progress = FALSE)
Two things are of note. If the user is specifying their own priors, bnec requires them to specify priors for all parameters. The pull_prior function shows the priors after the model was fitted, but suppose the user does not know what parameters were comprised in a particular model. In those instances, the user can call the function show_params(model = "all") to inspect the parameters of each function, or some targeted function in particular.
show_params(model = "nec3param") #> [[1]] #> y ~ top * exp(-exp(beta) * (x - nec) * step(x - nec)) #> top ~ 1 #> beta ~ 1 #> nec ~ 1
bayesmanecfit
The user can also specify a named list of priors when one or more models are being fitted to the same dataset.
set.seed(333) my_priors <- list(nec3param = c(prior_string("beta(5, 1)", nlpar = "top"), prior_string("normal(1.3, 2.7)", nlpar = "nec"), prior_string("gamma(0.5, 2)", nlpar = "beta")), nec4param = c(prior_string("beta(5, 1)", nlpar = "top"), prior_string("normal(1.3, 2.7)", nlpar = "nec"), prior_string("gamma(0.5, 2)", nlpar = "beta"), prior_string("beta(1, 5)", nlpar = "bot"))) exmp_c <- bnec(y ~ crf(x, model = c("nec3param", "nec4param")), data = nec_data, family = Beta(link = "identity"), prior = my_priors, iter = 1e4, control = list(adapt_delta = 0.99), open_progress = FALSE)
pull_prior also works for an object of class bayesmanecfit, as does check_priors which allows an option of passing a filename to save the prior and posterior probability density plots to a pdf.
pull_prior(exmp_c) check_priors(exmp_c, filename = "Check_priors")
#> $nec3param #> prior class coef group resp dpar nlpar lb ub source #> gamma(0.5, 2) b beta user #> gamma(0.5, 2) b Intercept beta (vectorized) #> normal(1.3, 2.7) b nec user #> normal(1.3, 2.7) b Intercept nec (vectorized) #> beta(5, 1) b top user #> beta(5, 1) b Intercept top (vectorized) #> gamma(0.01, 0.01) phi 0 default #> #> $nec4param #> prior class coef group resp dpar nlpar lb ub source #> gamma(0.5, 2) b beta user #> gamma(0.5, 2) b Intercept beta (vectorized) #> beta(1, 5) b bot user #> beta(1, 5) b Intercept bot (vectorized) #> normal(1.3, 2.7) b nec user #> normal(1.3, 2.7) b Intercept nec (vectorized) #> beta(5, 1) b top user #> beta(5, 1) b Intercept top (vectorized) #> gamma(0.01, 0.01) phi 0 default
The user can also specify priors for one model only out of the entire set, bayesnec will return a message stating that it is searching for priors on its own when they are either ill-formed (e.g. incomplete or have a typo), or the user simply decided not to specify priors for a particular model, e.g.
set.seed(333) my_priors <- list(nec3param = c(prior_string("beta(5, 1)", nlpar = "top"), prior_string("normal(1.3, 2.7)", nlpar = "nec"), prior_string("gamma(0.5, 2)", nlpar = "beta")), nec4param = c(prior_string("beta(5, 1)", nlpar = "top"), prior_string("normal(1.3, 2.7)", nlpar = "nec"), prior_string("gamma(0.5, 2)", nlpar = "beta"), prior_string("beta(1, 5)", nlpar = "bot"))) exmp_d <- bnec(y ~ crf(x, model = c("nec3param", "nec4param")), data = nec_data, family = Beta(link = "identity"), prior = my_priors[1], iter = 1e4, control = list(adapt_delta = 0.99), open_progress = FALSE)
prior = my_priors[[1]] would also have worked because the argument priors can either take a brmsprior object directly, or a named list containing model-specific brmsprior objects.
Finally the user can also extend an existing bayesmanecfit object with the function amend, also by specifying custom-built priors.
set.seed(333) necsigm_priors <- c(prior_string("beta(5, 1)", nlpar = "top"), prior_string("gamma(2, 6.5)", nlpar = "beta"), prior_string("normal(1.3, 2.7)", nlpar = "nec"), prior_string("normal(0, 2)", nlpar = "d")) exmp_e <- amend(exmp_d, add = "necsigm", priors = necsigm_priors)
Try the bayesnec package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.