Posterior inference for dose-response models

Fits the Bayesian hierarchical dose-response model of Miller et al. (2014) to multiple species using rjags (Plummer 2019) , and estimates posterior model probabilities using a Gibbs Variable Selection (GVS) approach (O'Hara and Sillanpää 2009) .

gibbs(
  dat,
  random.effects = FALSE,
  pseudo.n = 10000,
  mcmc.n = 1000,
  burnin = 1000,
  n.chains = 1,
  thin = 1,
  epsilon.upper = 30
)

Arguments

dat	Input data. Must be an object of class `rjtrace` or `brsdata`.
random.effects	Logical. When `TRUE`, uses a random effect model formulation.
pseudo.n	Number of iterations for the pseudo-priors.
mcmc.n	Number of posterior samples.
burnin	Number of iterations to discard as burn-in.
n.chains	Number of MCMC chains.
thin	Thinning interval.
epsilon.upper	Upper bound on the ε parameter used in the random effect model formulation.

Value

A list object of class gvs.

Details

Adapted from original code developed by Dina Sadykova as part of the Mocha project. The function can accommodate species/species groups either as a fixed or a random effect.

References

Miller PJO, Antunes RN, Wensveen PJ, Samarra FIP, Alves AC, Tyack PL, Kvadsheim PH, Kleivane L, Lam FA, Ainslie MA, Thomas L (2014). “Dose-response relationships for the onset of avoidance of sonar by free-ranging killer whales.” J. Acoust. Soc. Am., 135(2), 975. doi: 10.1121/1.4861346 .

O'Hara RB, Sillanpää MJ (2009). “A review of Bayesian variable selection methods: What, how and which.” Bayesian Analysis, 4(1), 85--117.

Plummer M (2019). rjags: Bayesian Graphical Models using MCMC. R package version 4-10, https://CRAN.R-project.org/package=rjags.

Author

Phil J. Bouchet

Examples

if (FALSE) {
library(espresso)

# Simulate data for two species
mydat <- simulate_data(n.species = 2, 
                       n.whales = 16, 
                       max.trials = 3, 
                       covariates = list(exposed = c(0, 5), range = 0.5),
                       mu = c(101, 158), 
                       phi = 20, 
                       sigma = 20, 
                       Rc = c(210, 211), 
                       seed = 58697)
summary(mydat)
     
# Model selection by GVS                        
gvs <- gibbs(dat = mydat, 
             random.effects = FALSE, 
             mcmc.n = 1000, 
             burnin = 500)
}